using Dates
using Pumas
using PumasUtilities
using DataFramesMeta
using PharmaDatasets
using CairoMakie
using AlgebraOfGraphics
using Random
Why are non-Gaussian random effects relevant?
1 Motivation - PK model
Why using a non-Gaussian distribution as the underlying distribution for the random effects? There are a couple of arguments.
First, the Gaussian distribution has unbounded support, i.e. it take any value in \((-\infty, \infty)\). While phamacokinetic parameters typically are (semi) bounded, e.g.:
- clearance and volumes, \((0, \infty)\)
- bioavailability, \([0, 1]\)
Additionally, in order for a Gaussian distribution to work as the underlying distribution, often we need to transform them (e.g. exponentiation and logistic transformation). But these transformations in some settings, when the random effects do not have a great impact, i.e. they do not have large values, may shift the mean of the typical values (\(\theta\)) so that the expectation of the typical values (\(\operatorname{E}\)) are not equal to the mean. For example, the following code block is a traditional 1-compartment PK model with a Gaussian random effect that needs to be constrained to positive values, \((0, \infty)\):
@random begin
ηCL ~ Normal(0.0, ωCL)
ηVc ~ Normal(0.0, ωVc)
end
@pre begin
CL = θCL * exp(ηCL)
Vc = θVc * exp(ηVc)
end
@dynamics begin
Central' = -CL / Vc * Central
endIf we recover the formula for the expectation of the log-normal distribution, we have that:
\[\operatorname{E}[CL] = \exp \left\{ \log(\theta_{CL}) + \frac{\omega^2_{CL}}{2} \right\} \approx \theta_{CL}\]
This approximation only holds for small \(\omega_{CL}\).
Hence, \(\theta_{CL}\) is only the typical value when \(\omega_{CL}\) is small.
Here is a small tabulation for \(\operatorname{E}[CL]\) when \(\theta_{CL} = 0.5\):
ωs = [0.1, 0.2, 0.4, 0.8, 1.6]
DataFrame(; ω_CL = ωs, E_CL = (ω -> exp(log(0.5) + ω^2 / 2)).(ωs))5×2 DataFrame
| Row | ω_CL | E_CL |
|---|---|---|
| Float64 | Float64 | |
| 1 | 0.1 | 0.502506 |
| 2 | 0.2 | 0.510101 |
| 3 | 0.4 | 0.541644 |
| 4 | 0.8 | 0.688564 |
| 5 | 1.6 | 1.79832 |
As you can see, the larger the \(\omega_{CL}\) the more \(\operatorname{E}[CL]\) deviates from \(\theta_{CL}\).
1.1 Gamma distribution for the rescue
We can use the gamma distribution which has the following parametrization:
\[\text{Gamma}(k, \theta)\]
where \(k\) is a shape parameter and \(\theta\) is a scale parameter.
Note
Shape parameters generally control the shape of the distribution rather than shifting it (as a location parameter) of stretching/shrinking it (as a scale parameter)
We can use an alternative parametrization where the mean-value appears directly a parameter:
\[\text{Gamma}(\mu, \sigma)\]
where:
- \(\mu = \theta k\)
- \(\sigma = k^{-\frac{1}{2}}\)
Tip
The \(\sigma\) parameter is the coefficient of variation, i.e.
\[\sigma = \frac{\operatorname{Var} X}{\operatorname{E} X},\]
because that mimics the role of \(\sigma\) in the LogNormal(log(μ), σ) where for small values of \(\sigma\)
\[\sigma \approx \\frac{\operatorname{Var} X}{\operatorname{E} X}.\]
So, our previous PK model now becomes:
@random begin
_CL ~ Gamma(1 / ωCL^2, θCL * ωCL^2)
_Vc ~ Gamma(1 / ωVc^2, θVc * ωVc^2)
end
@pre begin
CL = _CL
Vc = _Vc
end
@dynamics begin
Central' = -CL / Vc * Central
endAs you can see the mean from the gamma distribution becomes:
\[\operatorname{E}[CL] = \theta k = \frac{1}{\omega^2_{CL}} \theta_{CL} \omega^2_{CL} = \theta_{CL}\]
It does not dependent on the between-subject variability \(\omega\)!
Note
We are avoiding η notation here since we are modeling the subject-specific parameter directly.
1.2 Gamma versus Log-Nogmal Numerical Simulations
Before we dive into our PK examples, let us showcase the case for gamma versus log-normal with some numerical simulations.
First, let’s define a mean μ_PK value for a typical value along with an array of possible standard deviations σ values:
μ_PK = 1.0
σ = [0.1, 0.2, 0.5, 1.0, 1.5, 2.0]These will serve as the mean and standard deviations for our gamma and log-normal distributions.
Now let’s compare the coefficient of variation (CV) as a function of σ for LogNormal and Gamma:
num_df_gamma = DataFrame(;
μ = μ_PK,
σ = σ,
meanLogNormal = mean.(LogNormal.(log.(μ_PK), σ)),
cvLogNormal = std.(LogNormal.(log.(μ_PK), σ)) ./ mean.(LogNormal.(log.(μ_PK), σ)),
meanGamma = mean.(Gamma.(1 ./ σ .^ 2, μ_PK .* σ .^ 2)),
cvGamma = std.(Gamma.(1 ./ σ .^ 2, μ_PK .* σ .^ 2)) ./
mean.(Gamma.(1 ./ σ .^ 2, μ_PK .* σ .^ 2)),
)6×6 DataFrame
| Row | μ | σ | meanLogNormal | cvLogNormal | meanGamma | cvGamma |
|---|---|---|---|---|---|---|
| Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | |
| 1 | 1.0 | 0.1 | 1.00501 | 0.100251 | 1.0 | 0.1 |
| 2 | 1.0 | 0.2 | 1.0202 | 0.202017 | 1.0 | 0.2 |
| 3 | 1.0 | 0.5 | 1.13315 | 0.53294 | 1.0 | 0.5 |
| 4 | 1.0 | 1.0 | 1.64872 | 1.31083 | 1.0 | 1.0 |
| 5 | 1.0 | 1.5 | 3.08022 | 2.91337 | 1.0 | 1.5 |
| 6 | 1.0 | 2.0 | 7.38906 | 7.32108 | 1.0 | 2.0 |
f, ax, plotobj = lines(
num_df_gamma.σ,
num_df_gamma.meanLogNormal;
label = "μ - LogNormal",
linewidth = 3,
axis = (; xlabel = L"\sigma", ylabel = "value"),
)
lines!(
num_df_gamma.σ,
num_df_gamma.meanGamma;
label = "μ - Gamma",
linewidth = 3,
linestyle = :dashdot,
)
lines!(num_df_gamma.σ, num_df_gamma.cvLogNormal; label = "CV - LogNormal", linewidth = 3)
lines!(
num_df_gamma.σ,
num_df_gamma.cvGamma;
label = "CV - Gamma",
linewidth = 3,
linestyle = :dashdot,
)
axislegend(ax; position = :lt)
f
In the graph above, the dashed lines correspond to the mean and CV for the gamma distribution, whereas the solid lines correspond to the log-normal distribution.
There is clearly a bias in both the log-normal’s mean and CV that we don’t see in the gamma distribution.
2 Motivation - Bioavailability
Here is a very common model that can benefit from a non-Gaussian random effects distribution.
The model has one-compartment elimination and oral absorption with modeled bioavailability based on a crossover design.
The following code is a traditional PK model with a Gaussian random effect that needs to be constrained to the unit interval, \([0, 1]\):
@param begin
θF ∈ RealDomain(lower = 0.0, upper = 1.0)
ωF ∈ RealDomain(lower = 0.0)
end
@random begin
ηF ~ Normal(0.0, ωF)
end
@dosecontrol begin
bioav = (Depot = logistic(logit(θF) + ηF),)
endThe expectation \(\operatorname{E}[F]\) doesn’t have closed form and is generally different from \(\theta_F\). However, we have that:
\[\operatorname{E}[F] \approx \theta_F\]
when \(ωF\) is small. I.e. \(\theta_F\) is only the typical value when \(ωF\) is small.
2.1 Beta versus Logit-Normal Numerical Simulations
Let’s perform the same type of simulations we did before, but now we will be using the numerical integrator quadgk from the QuadGK.jl package. This is because we don’t have a closed form solution for \(\operatorname{E}[F]\) in the logit-normal parameterization.
using QuadGK: quadgkμ_bioav = 0.7We’ll also reuse the same σ values for the CVs.
num_df_beta = DataFrame(;
μ = μ_bioav,
σ = σ,
meanLogitNormal = map(
σ -> quadgk(
t -> logistic(t) * pdf(Normal(logit(μ_bioav), σ), t),
-100 * σ,
100 * σ,
)[1],
σ,
),
meanBeta = mean.(Beta.(μ_bioav ./ σ, (1 - μ_bioav) ./ σ)),
)6×4 DataFrame
| Row | μ | σ | meanLogitNormal | meanBeta |
|---|---|---|---|---|
| Float64 | Float64 | Float64 | Float64 | |
| 1 | 0.7 | 0.1 | 0.699582 | 0.7 |
| 2 | 0.7 | 0.2 | 0.698345 | 0.7 |
| 3 | 0.7 | 0.5 | 0.690393 | 0.7 |
| 4 | 0.7 | 1.0 | 0.668971 | 0.7 |
| 5 | 0.7 | 1.5 | 0.646064 | 0.7 |
| 6 | 0.7 | 2.0 | 0.626038 | 0.7 |
f, ax, plotobj = lines(
num_df_beta.σ,
num_df_beta.meanLogitNormal;
label = "μ - LogitNormal",
linewidth = 3,
axis = (; xlabel = L"\sigma", ylabel = "value"),
)
lines!(
num_df_beta.σ,
num_df_beta.meanBeta;
label = "μ - Beta",
linewidth = 3,
linestyle = :dashdot,
)
axislegend(ax; position = :lb)
f
In the graph above, the dashed lines correspond to the mean for the beta distribution, whereas the solid lines correspond to the logit-normal distribution.
As before, there is clearly a bias in the logit-normal’s mean that we don’t see in the beta distribution.
3 Warfarin data
We’ll demonstrate those intuitions using the Warfarin dataset.
pop = read_pumas(dataset("pumas/warfarin"))Population
Subjects: 32
Observations: dv4 Models and Simulations
Here we will provide a Gaussian and a non-Gaussian approach for:
- PK IV 1-compartment model fit for the Warfarin dataset
- Bioavaliability parallel absorption model simulation
4.1 Warfarin Gaussian and non-Gaussian PK model
The first model is a simple 1-compartment PK IV model with proportional error. This is for the Gaussian versus gamma random effects:
model_lognormal = @model begin
@param begin
θCL ∈ RealDomain(; lower = 0)
θVc ∈ RealDomain(; lower = 0)
ωCL ∈ RealDomain(; lower = 0)
ωVc ∈ RealDomain(; lower = 0)
σ ∈ RealDomain(; lower = 0)
end
@random begin
_CL ~ LogNormal(log(θCL), ωCL)
_Vc ~ LogNormal(log(θVc), ωVc)
end
# This is equivalent to defining
# CL = θCL*exp(ηCL)
# with
# ηCL = Normal(0, ωCL)
@pre begin
CL = _CL
Vc = _Vc
end
@dynamics Central1
@derived begin
μ := @. Central / Vc
dv ~ @. Normal(μ, μ * σ)
end
endPumasModel
Parameters: θCL, θVc, ωCL, ωVc, σ
Random effects: _CL, _Vc
Covariates:
Dynamical system variables: Central
Dynamical system type: Closed form
Derived: dv
Observed: dvmodel_gamma = @model begin
@param begin
θCL ∈ RealDomain(; lower = 0)
θVc ∈ RealDomain(; lower = 0)
ωCL ∈ RealDomain(; lower = 0)
ωVc ∈ RealDomain(; lower = 0)
σ ∈ RealDomain(; lower = 0)
end
@random begin
_CL ~ Gamma(1 / ωCL^2, θCL * ωCL^2)
_Vc ~ Gamma(1 / ωVc^2, θVc * ωVc^2)
end
@pre begin
CL = _CL
Vc = _Vc
end
@dynamics begin
Central' = -CL / Vc * Central
end
@derived begin
μ := @. Central / Vc
dv ~ @. Normal(μ, μ * σ)
end
endPumasModel
Parameters: θCL, θVc, ωCL, ωVc, σ
Random effects: _CL, _Vc
Covariates:
Dynamical system variables: Central
Dynamical system type: Matrix exponential
Derived: dv
Observed: dvWe also need some initial values for the fitting:
iparams_pk = (; θCL = 1.0, θVc = 5.0, ωCL = 0.1, ωVc = 0.1, σ = 0.2)(θCL = 1.0,
θVc = 5.0,
ωCL = 0.1,
ωVc = 0.1,
σ = 0.2,)We proceed by fitting both models:
fit_lognormal = fit(model_lognormal, pop, iparams_pk, FOCE())[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter Function value Gradient norm
0 5.770212e+03 7.912060e+03
* time: 0.027621030807495117
1 9.433464e+02 6.079483e+02
* time: 0.5276870727539062
2 8.189627e+02 4.423725e+02
* time: 0.5386300086975098
3 5.917683e+02 1.819248e+02
* time: 0.6032760143280029
4 5.421783e+02 1.121313e+02
* time: 0.6096658706665039
5 5.255651e+02 7.407230e+01
* time: 0.6155638694763184
6 5.208427e+02 8.699271e+01
* time: 0.6208839416503906
7 5.174883e+02 8.974584e+01
* time: 0.6265780925750732
8 5.138523e+02 7.328235e+01
* time: 0.6318020820617676
9 5.109883e+02 4.155805e+01
* time: 0.6370000839233398
10 5.094359e+02 3.170517e+01
* time: 0.641718864440918
11 5.086172e+02 3.327331e+01
* time: 0.646385908126831
12 5.080941e+02 2.942077e+01
* time: 0.6511268615722656
13 5.074009e+02 2.839941e+01
* time: 0.6558349132537842
14 5.059302e+02 3.330093e+01
* time: 0.660283088684082
15 5.036399e+02 3.172884e+01
* time: 0.7048540115356445
16 5.017004e+02 3.160020e+01
* time: 0.7099850177764893
17 5.008553e+02 2.599524e+01
* time: 0.7149169445037842
18 5.005913e+02 2.139314e+01
* time: 0.7196469306945801
19 5.003573e+02 2.134778e+01
* time: 0.7246279716491699
20 4.997249e+02 2.069868e+01
* time: 0.7294900417327881
21 4.984453e+02 1.859010e+01
* time: 0.7344329357147217
22 4.959584e+02 2.156209e+01
* time: 0.7397129535675049
23 4.923347e+02 3.030833e+01
* time: 0.7450759410858154
24 4.906916e+02 1.652278e+01
* time: 0.7504889965057373
25 4.902955e+02 6.360800e+00
* time: 0.7555508613586426
26 4.902870e+02 7.028603e+00
* time: 0.7609109878540039
27 4.902193e+02 1.176895e+00
* time: 0.7662758827209473
28 4.902189e+02 1.170642e+00
* time: 0.7936248779296875
29 4.902186e+02 1.167624e+00
* time: 0.797652006149292
30 4.902145e+02 1.110377e+00
* time: 0.8021199703216553
31 4.902079e+02 1.010507e+00
* time: 0.8066680431365967
32 4.901917e+02 9.619218e-01
* time: 0.8114030361175537
33 4.901683e+02 1.001006e+00
* time: 0.8157398700714111
34 4.901473e+02 6.138233e-01
* time: 0.8204829692840576
35 4.901412e+02 1.754342e-01
* time: 0.8254408836364746
36 4.901406e+02 2.617009e-02
* time: 0.8299710750579834
37 4.901405e+02 4.585882e-03
* time: 0.8338890075683594
38 4.901405e+02 7.668184e-04
* time: 0.8372058868408203FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Closed form
Log-likelihood value: -490.14052
Number of subjects: 32
Number of parameters: Fixed Optimized
0 5
Observation records: Active Missing
dv: 256 0
Total: 256 0
-----------------
Estimate
-----------------
θCL 0.16025
θVc 10.262
ωCL 0.23505
ωVc 0.10449
σ 0.3582
-----------------fit_gamma = fit(model_gamma, pop, iparams_pk, FOCE())[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter Function value Gradient norm
0 4.219519e+03 5.960578e+03
* time: 6.29425048828125e-5
1 9.644493e+02 5.541920e+02
* time: 0.014283895492553711
2 8.376083e+02 3.988904e+02
* time: 0.025664091110229492
3 6.503174e+02 1.495777e+02
* time: 0.03482508659362793
4 6.121241e+02 8.826583e+01
* time: 0.04348397254943848
5 5.977360e+02 9.144086e+01
* time: 0.08541107177734375
6 5.911620e+02 1.000933e+02
* time: 0.09413790702819824
7 5.858372e+02 8.423844e+01
* time: 0.1022040843963623
8 5.821934e+02 5.194402e+01
* time: 0.10989594459533691
9 5.801487e+02 3.461331e+01
* time: 0.1177968978881836
10 5.789092e+02 3.888113e+01
* time: 0.12574505805969238
11 5.780054e+02 3.556605e+01
* time: 0.13375401496887207
12 5.769455e+02 3.624436e+01
* time: 0.1414790153503418
13 5.749747e+02 4.322775e+01
* time: 0.14927911758422852
14 5.721322e+02 3.722515e+01
* time: 0.15687894821166992
15 5.695879e+02 3.401586e+01
* time: 0.16476011276245117
16 5.683277e+02 2.854997e+01
* time: 0.17324590682983398
17 5.678285e+02 2.644560e+01
* time: 0.21099209785461426
18 5.673305e+02 2.744429e+01
* time: 0.21860003471374512
19 5.662430e+02 2.793918e+01
* time: 0.22611594200134277
20 5.641877e+02 2.616169e+01
* time: 0.23360705375671387
21 5.606628e+02 2.257667e+01
* time: 0.24142909049987793
22 5.530616e+02 3.832878e+01
* time: 0.25009703636169434
23 5.528349e+02 5.518159e+01
* time: 0.2609899044036865
24 5.497231e+02 3.042064e+01
* time: 0.27074694633483887
25 5.488355e+02 6.929306e+00
* time: 0.2795851230621338
26 5.486095e+02 1.087865e+00
* time: 0.2891659736633301
27 5.486062e+02 6.456402e-01
* time: 0.2982749938964844
28 5.486061e+02 6.467689e-01
* time: 0.33109498023986816
29 5.486060e+02 6.463480e-01
* time: 0.33884692192077637
30 5.486055e+02 6.408914e-01
* time: 0.3466780185699463
31 5.486045e+02 6.208208e-01
* time: 0.35457706451416016
32 5.486020e+02 1.035462e+00
* time: 0.3626549243927002
33 5.485971e+02 1.452099e+00
* time: 0.37085795402526855
34 5.485897e+02 1.482593e+00
* time: 0.3788919448852539
35 5.485839e+02 8.420646e-01
* time: 0.38687610626220703
36 5.485822e+02 2.023876e-01
* time: 0.3948841094970703
37 5.485821e+02 1.885486e-02
* time: 0.40236902236938477
38 5.485821e+02 2.343932e-03
* time: 0.4097158908843994
39 5.485821e+02 4.417566e-04
* time: 0.41697192192077637FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Matrix exponential
Log-likelihood value: -548.58208
Number of subjects: 32
Number of parameters: Fixed Optimized
0 5
Observation records: Active Missing
dv: 256 0
Total: 256 0
-----------------
Estimate
-----------------
θCL 0.16466
θVc 10.329
ωCL 0.23348
ωVc 0.10661
σ 0.35767
-----------------Finally, let’s compare the estimates:
compare_estimates(; lognormal = fit_lognormal, gamma = fit_gamma)5×3 DataFrame
| Row | parameter | lognormal | gamma |
|---|---|---|---|
| String | Float64? | Float64? | |
| 1 | θCL | 0.160253 | 0.164658 |
| 2 | θVc | 10.2617 | 10.3288 |
| 3 | ωCL | 0.235046 | 0.233484 |
| 4 | ωVc | 0.10449 | 0.106611 |
| 5 | σ | 0.358205 | 0.357667 |
As mention above, the mean of a log-normal is \(\exp \left\{ \mu + \frac{\sigma^2}{2} \right\}\).
So let’s compare that with the gamma typical values:
DataFrame(;
parameter = ["θCL", "θVc"],
θLogNormal = [coef(fit_lognormal).θCL, coef(fit_lognormal).θVc],
ELogNormal = [
exp(log(coef(fit_lognormal).θCL) + coef(fit_lognormal).ωCL^2 / 2),
exp(log(coef(fit_lognormal).θVc) + coef(fit_lognormal).ωVc^2 / 2),
],
θGamma = [coef(fit_gamma).θCL, coef(fit_gamma).θVc],
)2×4 DataFrame
| Row | parameter | θLogNormal | ELogNormal | θGamma |
|---|---|---|---|---|
| String | Float64 | Float64 | Float64 | |
| 1 | θCL | 0.160253 | 0.164741 | 0.164658 |
| 2 | θVc | 10.2617 | 10.3178 | 10.3288 |
As you can see the Gaussian model has a slight bias in the estimation of both θCL and θVc.
Let’s also plot the two probability density functions (PDF) for θCL:
plotdataPK = @chain DataFrame(; x = range(0, 0.5; length = 1_000)) begin
@rtransform begin
:LogNormal =
pdf(LogNormal(log(coef(fit_lognormal).θCL), coef(fit_lognormal).ωCL), :x)
:Gamma = pdf(LogNormal(log(coef(fit_gamma).θCL), coef(fit_gamma).ωCL), :x)
end
end
first(plotdataPK, 5)5×3 DataFrame
| Row | x | LogNormal | Gamma |
|---|---|---|---|
| Float64 | Float64 | Float64 | |
| 1 | 0.0 | 0.0 | 0.0 |
| 2 | 0.000500501 | 5.27258e-128 | 5.2502e-131 |
| 3 | 0.001001 | 9.25648e-99 | 3.24235e-101 |
| 4 | 0.0015015 | 2.10132e-83 | 1.45529e-85 |
| 5 | 0.002002 | 2.71651e-73 | 2.9785e-75 |
data(stack(plotdataPK, [:LogNormal, :Gamma])) *
mapping(:x, :value; color = :variable) *
visual(Lines) |> draw
4.2 Bioavaliability Parallel Absorption Simulation
This is a parallel absorption model with bioavaliabity in both the “fast” as the “slow” depots.
First, the traditional approach with a logistic transformation of a Gaussian random variable. This makes the individual relative bioavailibility logit-normally distributed.
model_logitnormal = @model begin
@param begin
θkaFast ∈ RealDomain(; lower = 0)
θkaSlow ∈ RealDomain(; lower = 0)
θCL ∈ RealDomain(; lower = 0)
θVc ∈ RealDomain(; lower = 0)
θbioav ∈ RealDomain(; lower = 0.0, upper = 1.0)
ωCL ∈ RealDomain(; lower = 0)
ωVc ∈ RealDomain(; lower = 0)
# I call this one ξ to distinguish it from ω since the interpretation is NOT a relative error (coefficient of variation)
ξbioav ∈ RealDomain(; lower = 0, init = 0.1)
σ ∈ RealDomain(; lower = 0)
end
@random begin
_CL ~ Gamma(1 / ωCL^2, θCL * ωCL^2)
_Vc ~ Gamma(1 / ωVc^2, θVc * ωVc^2)
# define the latent Gaussian random effect. Notice the logit transform
ηbioavLogit ~ Normal(logit(θbioav), ξbioav)
end
@pre begin
kaFast = θkaFast
kaSlow = θkaSlow
CL = _CL
Vc = _Vc
end
@dosecontrol begin
# _bioav is LogitNormal distributed
_bioav = logistic(ηbioavLogit)
bioav = (; DepotFast = _bioav, DepotSlow = 1 - _bioav)
end
@dynamics begin
DepotFast' = -kaFast * DepotFast
DepotSlow' = -kaSlow * DepotSlow
Central' = kaFast * DepotFast + kaSlow * DepotSlow - CL / Vc * Central
end
@derived begin
μ := @. Central / Vc
dv ~ @. Normal(μ, abs(μ) * σ)
end
endPumasModel
Parameters: θkaFast, θkaSlow, θCL, θVc, θbioav, ωCL, ωVc, ξbioav, σ
Random effects: _CL, _Vc, ηbioavLogit
Covariates:
Dynamical system variables: DepotFast, DepotSlow, Central
Dynamical system type: Matrix exponential
Derived: dv
Observed: dvNow the same model but with the non-Gaussian random-effects using a beta distribution instead of the logit parameterization of the Gaussian distribution:
model_beta = @model begin
@param begin
θkaFast ∈ RealDomain(; lower = 0)
θkaSlow ∈ RealDomain(; lower = 0)
θCL ∈ RealDomain(; lower = 0)
θVc ∈ RealDomain(; lower = 0)
θbioav ∈ RealDomain(; lower = 0.0, upper = 1.0)
ωCL ∈ RealDomain(; lower = 0)
ωVc ∈ RealDomain(; lower = 0)
# We call this one n since the interpretation is like the length of a Binomial distribution
nbioav ∈ RealDomain(; lower = 0, init = 10)
σ ∈ RealDomain(; lower = 0)
end
@random begin
ηCL ~ Gamma(1 / ωCL^2, θCL * ωCL^2)
ηVc ~ Gamma(1 / ωVc^2, θVc * ωVc^2)
# The makes E(_bioav) = θbioav
# See https://en.wikipedia.org/wiki/Beta_distribution
ηbioav ~ Beta(θbioav * nbioav, (1 - θbioav) * nbioav)
end
@pre begin
kaFast = θkaFast
kaSlow = θkaSlow
CL = ηCL
Vc = ηVc
end
@dosecontrol begin
bioav = (; DepotFast = ηbioav, DepotSlow = 1 - ηbioav)
end
@dynamics begin
DepotFast' = -kaFast * DepotFast
DepotSlow' = -kaSlow * DepotSlow
Central' = kaFast * DepotFast + kaSlow * DepotSlow - CL / Vc * Central
end
@derived begin
μ := @. Central / Vc
dv ~ @. Normal(μ, abs(μ) * σ)
end
endPumasModel
Parameters: θkaFast, θkaSlow, θCL, θVc, θbioav, ωCL, ωVc, nbioav, σ
Random effects: ηCL, ηVc, ηbioav
Covariates:
Dynamical system variables: DepotFast, DepotSlow, Central
Dynamical system type: Matrix exponential
Derived: dv
Observed: dvWe have two types of random effects here.
First, as you are already familiar from the previous example, the clearance (CL), volume of concentration (Vc), and absorption rate (ka) have typical values (i.e. fixed effects) and between-subject variability (i.e. random effects) modelled as a gamma distribution.
Second, bioavailability Bioav is modelled as a beta distribution. Generally the beta distribution is parametrized as:
\[\text{Beta}(\alpha, \beta)\]
where both parameters \(\alpha\) and \(\beta\) are shape parameters.
One nice thing about the beta distribution is that it only takes values between and including 0 and 1, i.e. \([0, 1]\). This makes it the perfect candidate to model bioavailability parameters which are generally bounded in that interval. So, we don’t need to do a logistic transformation.
Another nice thing about the beta distribution is that we can use the alternative \((\mu, n)\)-parametrization with with \(\mu\) serving as a mean-value parameter:
\[\text{Beta}(\mu, n)\]
where in the original beta parametrization:
- \(\alpha = \mu n\)
- \(\beta = (1 - \mu) n\)
Hence, our mean is:
\[\operatorname{E}[F] = \mu = \theta_F\]
which, again, does not depend on any other parameters. The variance is
\[\operatorname{Var}(F) = \frac{\mu(1 - \mu)}{n}\]
so similar to the mean of Bernoulli trials.
Now let’s generate some data for the simulation:
dr = DosageRegimen(
DosageRegimen(100; cmt = :DepotFast),
DosageRegimen(100; cmt = :DepotSlow),
)2×10 DataFrame
| Row | time | cmt | amt | evid | ii | addl | rate | duration | ss | route |
|---|---|---|---|---|---|---|---|---|---|---|
| Float64 | Symbol | Float64 | Int8 | Float64 | Int64 | Float64 | Float64 | Int8 | NCA.Route | |
| 1 | 0.0 | DepotFast | 100.0 | 1 | 0.0 | 0 | 0.0 | 0.0 | 0 | NullRoute |
| 2 | 0.0 | DepotSlow | 100.0 | 1 | 0.0 | 0 | 0.0 | 0.0 | 0 | NullRoute |
simtimes = [0.5, 1.0, 2.0, 4.0, 8.0, 24.0]6-element Vector{Float64}:
0.5
1.0
2.0
4.0
8.0
24.0trueparam = (;
θkaFast = 0.9,
θkaSlow = 0.2,
θCL = 1.1,
θVc = 10.0,
θbioav = 0.7,
ωCL = 0.1,
ωVc = 0.1,
nbioav = 40,
σ = 0.1,
)(θkaFast = 0.9,
θkaSlow = 0.2,
θCL = 1.1,
θVc = 10.0,
θbioav = 0.7,
ωCL = 0.1,
ωVc = 0.1,
nbioav = 40,
σ = 0.1,)For simplicity, we just add 20% to the true values for initial values:
initparamBeta = map(t -> 1.2 * t, trueparam)(θkaFast = 1.08,
θkaSlow = 0.24,
θCL = 1.32,
θVc = 12.0,
θbioav = 0.84,
ωCL = 0.12,
ωVc = 0.12,
nbioav = 48.0,
σ = 0.12,)The initial values for the LogitNormal need to have ξbioav defined instead of nbioav:
initparamLogitNormal =
(Base.structdiff(initparamBeta, NamedTuple{(:nbioav,)})..., ξbioav = 0.1)(θkaFast = 1.08,
θkaSlow = 0.24,
θCL = 1.32,
θVc = 12.0,
θbioav = 0.84,
ωCL = 0.12,
ωVc = 0.12,
σ = 0.12,
ξbioav = 0.1,)Setup empty Subjects with the dose information:
skeletonpop = map(i -> Subject(; id = i, events = dr), 1:40)Population
Subjects: 40
Observations: Next, we simulate the data (while setting the seed for reprocibility):
simpop =
Subject.(
simobs(
model_beta,
skeletonpop,
trueparam;
obstimes = simtimes,
rng = Random.seed!(Random.default_rng(), 128),
)
)Finally let’s fit both models:
fit_logitnormal = fit(model_logitnormal, simpop, initparamLogitNormal, FOCE())[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter Function value Gradient norm
0 3.179853e+02 3.723266e+02
* time: 8.392333984375e-5
1 2.618605e+02 1.558889e+02
* time: 0.18386292457580566
2 2.484481e+02 5.652093e+01
* time: 0.28971099853515625
3 2.466691e+02 3.819960e+01
* time: 0.46693992614746094
4 2.430995e+02 2.281161e+01
* time: 0.5504388809204102
5 2.421492e+02 1.230731e+01
* time: 0.6352620124816895
6 2.412891e+02 1.022578e+01
* time: 0.7214739322662354
7 2.398220e+02 4.018354e+00
* time: 0.8125119209289551
8 2.398069e+02 2.747661e+00
* time: 0.9020669460296631
9 2.398011e+02 1.787833e+00
* time: 1.0361700057983398
10 2.397895e+02 2.923095e+00
* time: 1.1217210292816162
11 2.397511e+02 4.766103e+00
* time: 1.2086310386657715
12 2.397274e+02 2.427078e+00
* time: 1.298508882522583
13 2.396969e+02 2.317178e+00
* time: 1.3921010494232178
14 2.396432e+02 6.390959e+00
* time: 1.4863488674163818
15 2.395437e+02 1.431396e+01
* time: 1.5990240573883057
16 2.394242e+02 1.660862e+01
* time: 1.68772292137146
17 2.392883e+02 7.918201e+00
* time: 1.7752149105072021
18 2.392648e+02 5.359476e+00
* time: 1.883673906326294
19 2.392551e+02 8.504104e-01
* time: 1.9758200645446777
20 2.392536e+02 8.658366e-01
* time: 2.0671210289001465
21 2.392502e+02 1.359491e+00
* time: 2.1751298904418945
22 2.392455e+02 1.908257e+00
* time: 2.2630770206451416
23 2.392352e+02 3.434282e+00
* time: 2.3530189990997314
24 2.392094e+02 4.019057e+00
* time: 2.4478719234466553
25 2.391904e+02 2.360120e+00
* time: 2.5586109161376953
26 2.391639e+02 2.232625e+00
* time: 2.687286853790283
27 2.390985e+02 7.831580e+00
* time: 2.7756969928741455
28 2.390462e+02 1.220982e+01
* time: 2.8644769191741943
29 2.389849e+02 1.372374e+01
* time: 2.9748599529266357
30 2.389295e+02 1.012388e+01
* time: 3.068449020385742
31 2.389050e+02 7.364830e+00
* time: 3.200345993041992
32 2.388675e+02 3.493531e+00
* time: 3.2915329933166504
33 2.388348e+02 3.601972e+00
* time: 3.3806920051574707
34 2.388092e+02 9.896791e+00
* time: 3.4761359691619873
35 2.387816e+02 5.816471e+00
* time: 3.5709638595581055
36 2.387667e+02 3.995305e+00
* time: 3.664850950241089
37 2.387624e+02 2.111619e+00
* time: 3.7786519527435303
38 2.387542e+02 3.340445e+00
* time: 3.866576910018921
39 2.387468e+02 3.742750e+00
* time: 3.9536659717559814
40 2.387148e+02 8.270335e+00
* time: 4.045107841491699
41 2.386920e+02 8.042091e+00
* time: 4.135984897613525
42 2.386311e+02 1.786257e+00
* time: 4.228839874267578
43 2.385982e+02 4.295253e+00
* time: 4.337578058242798
44 2.385769e+02 2.426177e+00
* time: 4.425214052200317
45 2.385642e+02 4.507279e-01
* time: 4.51392388343811
46 2.385598e+02 1.783922e+00
* time: 4.604898929595947
47 2.385584e+02 4.345023e-01
* time: 4.695162057876587
48 2.385582e+02 3.504116e-01
* time: 4.783501863479614
49 2.385581e+02 9.328221e-02
* time: 4.888208866119385
50 2.385580e+02 1.265801e-02
* time: 4.969329833984375
51 2.385580e+02 2.371756e-03
* time: 5.048222064971924
52 2.385580e+02 2.371756e-03
* time: 5.1578710079193115
53 2.385580e+02 2.371756e-03
* time: 5.280717849731445
54 2.385580e+02 2.371756e-03
* time: 5.416619062423706
55 2.385580e+02 2.371756e-03
* time: 5.532060861587524
56 2.385580e+02 2.371756e-03
* time: 5.648100852966309
57 2.385580e+02 2.371756e-03
* time: 5.780947923660278
58 2.385580e+02 2.371756e-03
* time: 5.823861837387085FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Matrix exponential
Log-likelihood value: -238.55805
Number of subjects: 40
Number of parameters: Fixed Optimized
0 9
Observation records: Active Missing
dv: 240 0
Total: 240 0
---------------------
Estimate
---------------------
θkaFast 1.8907
θkaSlow 0.60771
θCL 1.0797
θVc 11.375
θbioav 0.14394
ωCL 0.08151
ωVc 0.10157
ξbioav 0.12513
σ 0.10441
---------------------fit_beta = fit(model_beta, simpop, initparamBeta, FOCE())[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter Function value Gradient norm
0 3.554658e+02 3.718530e+02
* time: 6.890296936035156e-5
1 2.991160e+02 1.342892e+02
* time: 0.3098170757293701
2 2.851580e+02 4.511285e+01
* time: 0.43611598014831543
3 2.835992e+02 4.343268e+01
* time: 0.5589299201965332
4 2.821916e+02 4.373720e+01
* time: 0.713310956954956
5 2.794524e+02 1.253857e+01
* time: 0.8081269264221191
6 2.787960e+02 1.102093e+01
* time: 0.9012460708618164
7 2.770092e+02 1.114240e+01
* time: 1.0011999607086182
8 2.769202e+02 5.089120e+00
* time: 1.103287935256958
9 2.769088e+02 3.041932e+00
* time: 1.2042279243469238
10 2.768987e+02 1.603805e+00
* time: 1.333038091659546
11 2.767766e+02 1.411814e+01
* time: 1.4287359714508057
12 2.766170e+02 2.719464e+01
* time: 1.527318000793457
13 2.764580e+02 3.014173e+01
* time: 1.6300499439239502
14 2.762515e+02 1.669149e+01
* time: 1.7316360473632812
15 2.761434e+02 1.957297e+00
* time: 1.8534669876098633
16 2.761353e+02 3.756442e+00
* time: 1.9489850997924805
17 2.761237e+02 1.341491e+00
* time: 2.0422511100769043
18 2.761152e+02 2.182772e+00
* time: 2.1422600746154785
19 2.761036e+02 3.209881e+00
* time: 2.243079900741577
20 2.760912e+02 3.070991e+00
* time: 2.367645025253296
21 2.760762e+02 1.990495e+00
* time: 2.464064121246338
22 2.760587e+02 1.586627e+00
* time: 2.5609889030456543
23 2.760290e+02 2.471592e+00
* time: 2.662287950515747
24 2.759737e+02 1.040809e+01
* time: 2.7668159008026123
25 2.758948e+02 9.671126e+00
* time: 2.8692409992218018
26 2.758264e+02 7.207553e+00
* time: 2.991995096206665
27 2.757728e+02 3.300455e+00
* time: 3.109786033630371
28 2.757344e+02 4.723183e+00
* time: 3.207331895828247
29 2.756731e+02 6.515451e+00
* time: 3.309117078781128
30 2.755936e+02 4.181865e+00
* time: 3.4124910831451416
31 2.755684e+02 5.163221e+00
* time: 3.558938980102539
32 2.755454e+02 8.132621e+00
* time: 3.656682014465332
33 2.755346e+02 2.188011e+00
* time: 3.7570860385894775
34 2.755279e+02 3.127135e+00
* time: 3.8586840629577637
35 2.755221e+02 3.813163e+00
* time: 3.9574599266052246
36 2.754865e+02 6.610197e+00
* time: 4.078087091445923
37 2.754663e+02 5.265615e+00
* time: 4.174292087554932
38 2.754486e+02 9.998959e-01
* time: 4.273252964019775
39 2.754421e+02 1.112426e+00
* time: 4.373258113861084
40 2.754395e+02 3.014869e-01
* time: 4.473398923873901
41 2.754386e+02 5.351895e-01
* time: 4.5944108963012695
42 2.754384e+02 1.402404e-01
* time: 4.6871631145477295
43 2.754383e+02 1.288533e-01
* time: 4.779293060302734
44 2.754381e+02 3.845371e-01
* time: 4.87625789642334
45 2.754376e+02 9.172184e-01
* time: 4.974968910217285
46 2.754363e+02 1.766701e+00
* time: 5.073789119720459
47 2.754331e+02 3.007560e+00
* time: 5.19330096244812
48 2.754256e+02 4.573254e+00
* time: 5.288755893707275
49 2.754117e+02 5.874317e+00
* time: 5.3864099979400635
50 2.753921e+02 6.633185e+00
* time: 5.48725700378418
51 2.753744e+02 2.205111e+00
* time: 5.588349103927612
52 2.753702e+02 2.363930e+00
* time: 5.711046934127808
53 2.753546e+02 8.713730e-01
* time: 5.807805061340332
54 2.753543e+02 3.795710e+00
* time: 5.903501987457275
55 2.753492e+02 1.784042e+00
* time: 6.0016560554504395
56 2.753441e+02 8.514295e-01
* time: 6.101538896560669
57 2.753339e+02 1.341633e+00
* time: 6.225685119628906
58 2.753266e+02 2.387419e+00
* time: 6.322042942047119
59 2.753233e+02 1.246148e+00
* time: 6.417129993438721
60 2.753208e+02 2.973353e-01
* time: 6.515166997909546
61 2.753202e+02 3.696236e-01
* time: 6.613790988922119
62 2.753193e+02 4.553497e-01
* time: 6.7126030921936035
63 2.753186e+02 2.081281e-01
* time: 6.83256196975708
64 2.753183e+02 8.672527e-02
* time: 6.9260571002960205
65 2.753182e+02 1.015286e-01
* time: 7.017467975616455
66 2.753181e+02 1.003887e-01
* time: 7.113807916641235
67 2.753181e+02 5.060689e-02
* time: 7.2108070850372314
68 2.753180e+02 2.051541e-02
* time: 7.30561900138855
69 2.753180e+02 1.864290e-02
* time: 7.418317079544067
70 2.753180e+02 1.911113e-02
* time: 7.507455110549927
71 2.753180e+02 8.991085e-03
* time: 7.597229957580566
72 2.753180e+02 4.925215e-03
* time: 7.690983057022095
73 2.753180e+02 2.135748e-03
* time: 7.783713102340698
74 2.753180e+02 2.135926e-03
* time: 7.940264940261841FittedPumasModel
Successful minimization: false
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Matrix exponential
Log-likelihood value: -275.31802
Number of subjects: 40
Number of parameters: Fixed Optimized
0 9
Observation records: Active Missing
dv: 240 0
Total: 240 0
----------------------
Estimate
----------------------
θkaFast 1.9882
θkaSlow 0.61333
θCL 1.0797
θVc 11.381
θbioav 0.13274
ωCL 0.081471
ωVc 0.10178
nbioav 1.4823e7
σ 0.10464
----------------------As before, let’s compare the estimates:
compare_estimates(; logitnormal = fit_logitnormal, beta = fit_beta)10×3 DataFrame
| Row | parameter | logitnormal | beta |
|---|---|---|---|
| String | Float64? | Float64? | |
| 1 | θkaFast | 1.89067 | 1.98819 |
| 2 | θkaSlow | 0.607708 | 0.613333 |
| 3 | θCL | 1.07968 | 1.07967 |
| 4 | θVc | 11.3745 | 11.3811 |
| 5 | θbioav | 0.143942 | 0.132743 |
| 6 | ωCL | 0.0815101 | 0.081471 |
| 7 | ωVc | 0.101565 | 0.101784 |
| 8 | σ | 0.10441 | 0.104644 |
| 9 | ξbioav | 0.125128 | missing |
| 10 | nbioav | missing | 1.48228e7 |
Again, we’ll both PDFs from the estimated values:
plotdatabioav = @chain DataFrame(; x = range(0, 1; length = 1_000)) begin
@rtransform begin
:logitnormal =
1 / coef(fit_logitnormal).ξbioav / √(2π) / (:x * (1 - :x)) * exp(
-(logit(:x) - logit(coef(fit_logitnormal).θbioav))^2 /
(2 * coef(fit_logitnormal).ξbioav^2),
)
:beta = pdf(
Beta(
coef(fit_beta).θbioav * coef(fit_beta).nbioav,
(1 - coef(fit_beta).θbioav) * coef(fit_beta).nbioav,
),
:x,
)
end
end
first(plotdatabioav, 5)5×3 DataFrame
| Row | x | logitnormal | beta |
|---|---|---|---|
| Float64 | Float64 | Float64 | |
| 1 | 0.0 | NaN | 0.0 |
| 2 | 0.001001 | 0.0 | 0.0 |
| 3 | 0.002002 | 1.54002e-269 | 0.0 |
| 4 | 0.003003 | 4.49003e-222 | 0.0 |
| 5 | 0.004004 | 3.80322e-191 | 0.0 |
plt_pdf_bioav =
data(stack(plotdatabioav, [:logitnormal, :beta])) *
mapping(:x, :value; color = :variable) *
visual(Lines);
draw(plt_pdf_bioav; axis = (; xticks = 0.1:0.1:1.0))
For this dataset, the two distributions differ significantly with the Beta model producing a distribution much closer to the truth but for other realizations of the simulated data they are closer to each other.