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.03017592430114746
1 9.433464e+02 6.079483e+02
* time: 1.3123559951782227
2 8.189627e+02 4.423725e+02
* time: 1.3510699272155762
3 5.917683e+02 1.819248e+02
* time: 1.3806648254394531
4 5.421783e+02 1.121313e+02
* time: 1.4042999744415283
5 5.255651e+02 7.407230e+01
* time: 1.4252288341522217
6 5.208427e+02 8.699271e+01
* time: 1.4452698230743408
7 5.174883e+02 8.974584e+01
* time: 1.4661200046539307
8 5.138523e+02 7.328235e+01
* time: 1.485672950744629
9 5.109883e+02 4.155805e+01
* time: 1.502781867980957
10 5.094359e+02 3.170517e+01
* time: 1.5186269283294678
11 5.086172e+02 3.327331e+01
* time: 1.5351648330688477
12 5.080941e+02 2.942077e+01
* time: 1.550691843032837
13 5.074009e+02 2.839941e+01
* time: 1.5666368007659912
14 5.059302e+02 3.330093e+01
* time: 1.5827698707580566
15 5.036399e+02 3.172884e+01
* time: 1.5999388694763184
16 5.017004e+02 3.160020e+01
* time: 1.6177518367767334
17 5.008553e+02 2.599524e+01
* time: 1.634917974472046
18 5.005913e+02 2.139314e+01
* time: 1.6508688926696777
19 5.003573e+02 2.134778e+01
* time: 1.6665170192718506
20 4.997249e+02 2.069868e+01
* time: 1.6827068328857422
21 4.984453e+02 1.859010e+01
* time: 1.6993370056152344
22 4.959584e+02 2.156209e+01
* time: 1.7159039974212646
23 4.923347e+02 3.030833e+01
* time: 1.8090839385986328
24 4.906916e+02 1.652278e+01
* time: 1.8293588161468506
25 4.902955e+02 6.360800e+00
* time: 1.8474438190460205
26 4.902870e+02 7.028603e+00
* time: 1.8665649890899658
27 4.902193e+02 1.176895e+00
* time: 1.8847289085388184
28 4.902189e+02 1.170642e+00
* time: 1.8991680145263672
29 4.902186e+02 1.167624e+00
* time: 1.9116218090057373
30 4.902145e+02 1.110377e+00
* time: 1.9273109436035156
31 4.902079e+02 1.010507e+00
* time: 1.9430749416351318
32 4.901917e+02 9.619218e-01
* time: 1.9590568542480469
33 4.901683e+02 1.001006e+00
* time: 1.9728479385375977
34 4.901473e+02 6.138233e-01
* time: 1.9870309829711914
35 4.901412e+02 1.754342e-01
* time: 2.001046895980835
36 4.901406e+02 2.617009e-02
* time: 2.0136919021606445
37 4.901405e+02 4.585882e-03
* time: 2.0244979858398438
38 4.901405e+02 7.668184e-04
* time: 2.03357195854187FittedPumasModel
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.160707e+03 5.960578e+03
* time: 0.00010609626770019531
1 9.056373e+02 5.541920e+02
* time: 0.042034149169921875
2 7.787963e+02 3.988904e+02
* time: 0.07386994361877441
3 5.915054e+02 1.495777e+02
* time: 0.09892892837524414
4 5.533120e+02 8.826583e+01
* time: 0.12057805061340332
5 5.389239e+02 9.144086e+01
* time: 0.140625
6 5.323499e+02 1.000933e+02
* time: 0.16185307502746582
7 5.270252e+02 8.423844e+01
* time: 0.18349695205688477
8 5.233813e+02 5.194402e+01
* time: 0.20279407501220703
9 5.213366e+02 3.461331e+01
* time: 0.22287607192993164
10 5.200972e+02 3.888113e+01
* time: 0.24251699447631836
11 5.191933e+02 3.556605e+01
* time: 0.26099514961242676
12 5.181335e+02 3.624436e+01
* time: 0.2797069549560547
13 5.161626e+02 4.322775e+01
* time: 0.2997879981994629
14 5.133202e+02 3.722515e+01
* time: 0.31940698623657227
15 5.107758e+02 3.401586e+01
* time: 0.33930301666259766
16 5.095157e+02 2.854997e+01
* time: 0.3600029945373535
17 5.090165e+02 2.644560e+01
* time: 0.37935400009155273
18 5.085184e+02 2.744429e+01
* time: 0.40030503273010254
19 5.074309e+02 2.793918e+01
* time: 0.4205441474914551
20 5.053757e+02 2.616169e+01
* time: 0.4399740695953369
21 5.018507e+02 2.257667e+01
* time: 0.4599001407623291
22 4.942495e+02 3.832878e+01
* time: 0.48278093338012695
23 4.940229e+02 5.518159e+01
* time: 0.5145440101623535
24 4.909110e+02 3.042064e+01
* time: 0.5399010181427002
25 4.900234e+02 6.929306e+00
* time: 0.561089038848877
26 4.897974e+02 1.087865e+00
* time: 0.5819461345672607
27 4.897942e+02 6.456402e-01
* time: 0.6017501354217529
28 4.897940e+02 6.467689e-01
* time: 0.6189711093902588
29 4.897939e+02 6.463480e-01
* time: 0.6362340450286865
30 4.897935e+02 6.408914e-01
* time: 0.653188943862915
31 4.897924e+02 6.208208e-01
* time: 0.6707279682159424
32 4.897900e+02 1.035462e+00
* time: 0.6889641284942627
33 4.897850e+02 1.452099e+00
* time: 0.8061599731445312
34 4.897776e+02 1.482593e+00
* time: 0.824099063873291
35 4.897718e+02 8.420646e-01
* time: 0.8426380157470703
36 4.897702e+02 2.023876e-01
* time: 0.8613109588623047
37 4.897700e+02 1.885486e-02
* time: 0.8782949447631836
38 4.897700e+02 2.343932e-03
* time: 0.8933961391448975
39 4.897700e+02 4.417566e-04
* time: 0.9061450958251953FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Matrix exponential
Log-likelihood value: -489.77002
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):
Random.seed!(128)
simpop = Subject.(simobs(model_beta, skeletonpop, trueparam; obstimes = simtimes))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 2.512964e+02 3.650794e+02
* time: 7.104873657226562e-5
1 2.005689e+02 1.252637e+02
* time: 0.3348069190979004
2 1.875728e+02 4.225518e+01
* time: 0.4770019054412842
3 1.863803e+02 3.429785e+01
* time: 1.7753379344940186
4 1.845581e+02 3.694229e+01
* time: 1.8780949115753174
5 1.828416e+02 1.170915e+01
* time: 1.9805359840393066
6 1.823277e+02 1.004017e+01
* time: 2.0814740657806396
7 1.810629e+02 5.326780e+00
* time: 2.1874430179595947
8 1.810479e+02 1.767987e+00
* time: 2.2868549823760986
9 1.810272e+02 3.852037e+00
* time: 2.394360065460205
10 1.809414e+02 1.196625e+01
* time: 2.506303071975708
11 1.806489e+02 1.861440e+01
* time: 2.6196470260620117
12 1.801984e+02 1.361774e+01
* time: 3.015937089920044
13 1.800427e+02 2.177039e+01
* time: 3.196716070175171
14 1.796554e+02 7.079039e+00
* time: 3.300755023956299
15 1.795832e+02 1.581499e+01
* time: 3.4056270122528076
16 1.795220e+02 6.552120e+00
* time: 3.5103049278259277
17 1.794621e+02 6.612645e+00
* time: 3.6206350326538086
18 1.793643e+02 6.603903e+00
* time: 3.7496819496154785
19 1.793502e+02 1.025787e+01
* time: 3.8529539108276367
20 1.793178e+02 1.202199e+01
* time: 3.9573590755462646
21 1.792424e+02 1.625661e+01
* time: 4.063184022903442
22 1.791560e+02 1.128856e+01
* time: 4.169548034667969
23 1.790991e+02 7.625637e+00
* time: 4.280467987060547
24 1.790756e+02 8.557258e+00
* time: 4.403537034988403
25 1.790633e+02 4.848482e+00
* time: 4.506968021392822
26 1.790408e+02 5.859306e+00
* time: 4.608319997787476
27 1.789734e+02 1.106035e+01
* time: 4.712269067764282
28 1.789070e+02 1.143361e+01
* time: 4.8209240436553955
29 1.788321e+02 7.098448e+00
* time: 4.936424970626831
30 1.787896e+02 5.709857e+00
* time: 5.06594705581665
31 1.787809e+02 7.456555e+00
* time: 5.172116041183472
32 1.787671e+02 7.320931e-01
* time: 5.277366876602173
33 1.787660e+02 5.023990e-01
* time: 5.380575895309448
34 1.787656e+02 3.813994e-01
* time: 5.489840030670166
35 1.787639e+02 8.608909e-01
* time: 5.617810010910034
36 1.787607e+02 2.047321e+00
* time: 5.7222900390625
37 1.787525e+02 3.859529e+00
* time: 5.8261260986328125
38 1.787349e+02 5.864920e+00
* time: 5.934006929397583
39 1.787017e+02 6.966045e+00
* time: 6.04525089263916
40 1.786574e+02 5.348939e+00
* time: 6.166522026062012
41 1.786270e+02 1.642025e+00
* time: 6.301836013793945
42 1.786181e+02 5.211823e-01
* time: 6.406646966934204
43 1.786153e+02 1.186061e+00
* time: 6.512150049209595
44 1.786125e+02 1.292005e+00
* time: 6.614824056625366
45 1.786099e+02 7.814598e-01
* time: 6.719577074050903
46 1.786086e+02 1.369456e-01
* time: 6.833713054656982
47 1.786082e+02 1.912170e-01
* time: 6.9619269371032715
48 1.786080e+02 2.670802e-01
* time: 7.0650599002838135
49 1.786078e+02 1.979262e-01
* time: 7.166836977005005
50 1.786077e+02 5.177918e-02
* time: 7.267877101898193
51 1.786076e+02 2.998328e-02
* time: 7.370927095413208
52 1.786076e+02 5.799706e-02
* time: 7.478882074356079
53 1.786076e+02 4.280434e-02
* time: 7.604887008666992
54 1.786076e+02 1.467829e-02
* time: 7.7057459354400635
55 1.786076e+02 6.928178e-03
* time: 7.799515008926392
56 1.786076e+02 1.236015e-02
* time: 7.89811110496521
57 1.786076e+02 9.950826e-03
* time: 7.995542049407959
58 1.786076e+02 2.974370e-03
* time: 8.097429037094116
59 1.786076e+02 1.374576e-03
* time: 8.202800035476685
60 1.786076e+02 2.860078e-03
* time: 8.314227104187012
61 1.786076e+02 2.100127e-03
* time: 8.404349088668823
62 1.786076e+02 2.100127e-03
* time: 8.538606882095337
63 1.786076e+02 6.412834e-03
* time: 8.650691032409668
64 1.786076e+02 6.412834e-03
* time: 8.823524951934814
65 1.786076e+02 6.412834e-03
* time: 9.048969030380249
66 1.786076e+02 6.412834e-03
* time: 9.355406999588013
67 1.786076e+02 6.412834e-03
* time: 9.624810934066772FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Matrix exponential
Log-likelihood value: -178.60759
Number of subjects: 40
Number of parameters: Fixed Optimized
0 9
Observation records: Active Missing
dv: 240 0
Total: 240 0
-----------------------
Estimate
-----------------------
θkaFast 0.91097
θkaSlow 0.13112
θCL 1.0854
θVc 7.1008
θbioav 0.4802
ωCL 0.088113
ωVc 0.12133
ξbioav 1.8429e-5
σ 0.10545
-----------------------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 2.523111e+02 3.644346e+02
* time: 7.987022399902344e-5
1 2.014577e+02 1.265001e+02
* time: 0.4730849266052246
2 1.880885e+02 4.190708e+01
* time: 0.703819990158081
3 1.870317e+02 8.825666e+01
* time: 0.9377110004425049
4 1.846027e+02 4.400156e+01
* time: 1.1289198398590088
5 1.834445e+02 1.906624e+01
* time: 1.287186861038208
6 1.828599e+02 1.113882e+01
* time: 1.4451210498809814
7 1.815719e+02 7.449355e+00
* time: 1.6241929531097412
8 1.815131e+02 2.164678e+00
* time: 1.811180830001831
9 1.814896e+02 2.167319e+00
* time: 2.0310909748077393
10 1.814458e+02 4.615738e+00
* time: 2.193889856338501
11 1.813173e+02 9.576967e+00
* time: 2.357351064682007
12 1.809756e+02 2.052077e+01
* time: 2.52400803565979
13 1.807625e+02 4.553366e+01
* time: 2.690600872039795
14 1.802224e+02 6.550892e+00
* time: 2.865304946899414
15 1.800862e+02 2.865509e+00
* time: 3.046628952026367
16 1.800780e+02 1.164611e+00
* time: 3.217371940612793
17 1.800737e+02 7.952462e-01
* time: 3.3608429431915283
18 1.800352e+02 4.860618e+00
* time: 3.514164924621582
19 1.800089e+02 5.176689e+00
* time: 3.6559040546417236
20 1.799679e+02 4.303892e+00
* time: 3.8076770305633545
21 1.799153e+02 4.612832e+00
* time: 3.963465929031372
22 1.798423e+02 1.209387e+01
* time: 4.146183967590332
23 1.796821e+02 1.712256e+01
* time: 4.286916971206665
24 1.794275e+02 1.435100e+01
* time: 4.428407907485962
25 1.793773e+02 4.137313e+00
* time: 4.566288948059082
26 1.793488e+02 1.846545e+00
* time: 4.75635290145874
27 1.793331e+02 5.502533e+00
* time: 4.908764839172363
28 1.793234e+02 2.894037e+00
* time: 5.075472831726074
29 1.793119e+02 1.453372e+00
* time: 5.210052967071533
30 1.792879e+02 4.109884e+00
* time: 5.345427989959717
31 1.792599e+02 4.613173e+00
* time: 5.477136850357056
32 1.792384e+02 4.254549e+00
* time: 5.616346836090088
33 1.792251e+02 4.415024e+00
* time: 5.75923490524292
34 1.792026e+02 3.229145e+00
* time: 5.936112880706787
35 1.791841e+02 3.422094e+00
* time: 6.079349040985107
36 1.791705e+02 1.621228e+00
* time: 6.22663688659668
37 1.791555e+02 3.462667e+00
* time: 6.374189853668213
38 1.791293e+02 5.586685e+00
* time: 6.5241379737854
39 1.790713e+02 9.927638e+00
* time: 6.6753740310668945
40 1.789891e+02 1.133271e+01
* time: 6.830858945846558
41 1.788585e+02 1.279172e+01
* time: 7.005619049072266
42 1.787407e+02 5.291681e+00
* time: 7.151219844818115
43 1.786633e+02 5.367971e+00
* time: 7.29197096824646
44 1.786405e+02 2.571548e+00
* time: 7.425806999206543
45 1.786339e+02 2.720314e+00
* time: 7.569604873657227
46 1.786252e+02 1.930583e+00
* time: 7.704701900482178
47 1.786182e+02 1.523164e+00
* time: 7.854750871658325
48 1.786126e+02 5.027920e-01
* time: 7.977946043014526
49 1.786100e+02 3.733023e-01
* time: 8.094667911529541
50 1.786089e+02 2.749553e-01
* time: 8.21721887588501
51 1.786083e+02 1.981772e-01
* time: 8.339988946914673
52 1.786079e+02 1.005262e-01
* time: 8.506510019302368
53 1.786078e+02 2.549641e-02
* time: 8.674316883087158
54 1.786077e+02 4.452931e-02
* time: 8.84359884262085
55 1.786076e+02 2.967125e-02
* time: 8.968178033828735
56 1.786076e+02 1.068274e-02
* time: 9.090214967727661
57 1.786076e+02 5.355447e-03
* time: 9.210286855697632
58 1.786076e+02 8.180920e-03
* time: 9.32889199256897
59 1.786076e+02 4.935439e-03
* time: 9.45156192779541
60 1.786076e+02 4.387986e-03
* time: 9.599244832992554
61 1.786076e+02 4.388023e-03
* time: 9.800833940505981
62 1.786076e+02 4.387934e-03
* time: 9.958776950836182
63 1.786076e+02 4.387934e-03
* time: 10.087278842926025FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Matrix exponential
Log-likelihood value: -178.60759
Number of subjects: 40
Number of parameters: Fixed Optimized
0 9
Observation records: Active Missing
dv: 240 0
Total: 240 0
----------------------
Estimate
----------------------
θkaFast 0.91099
θkaSlow 0.13112
θCL 1.0854
θVc 7.1007
θbioav 0.48019
ωCL 0.088114
ωVc 0.12134
nbioav 2.7333e7
σ 0.10545
----------------------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 | 0.910971 | 0.910988 |
| 2 | θkaSlow | 0.131115 | 0.131117 |
| 3 | θCL | 1.0854 | 1.0854 |
| 4 | θVc | 7.10076 | 7.10073 |
| 5 | θbioav | 0.480202 | 0.480191 |
| 6 | ωCL | 0.0881132 | 0.0881137 |
| 7 | ωVc | 0.121334 | 0.121335 |
| 8 | σ | 0.105448 | 0.105449 |
| 9 | ξbioav | 1.84292e-5 | missing |
| 10 | nbioav | missing | 2.73333e7 |
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 | 0.0 | 0.0 |
| 4 | 0.003003 | 0.0 | 0.0 |
| 5 | 0.004004 | 0.0 | 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.