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 variables: Central
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 variables: Central
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.01983499526977539
1 9.433464e+02 6.079483e+02
* time: 0.3522830009460449
2 8.189627e+02 4.423725e+02
* time: 0.3585999011993408
3 5.917683e+02 1.819248e+02
* time: 0.39913392066955566
4 5.421783e+02 1.121313e+02
* time: 0.40285205841064453
5 5.255651e+02 7.407230e+01
* time: 0.40619587898254395
6 5.208427e+02 8.699271e+01
* time: 0.4095308780670166
7 5.174883e+02 8.974584e+01
* time: 0.4128880500793457
8 5.138523e+02 7.328235e+01
* time: 0.4165370464324951
9 5.109883e+02 4.155805e+01
* time: 0.42019009590148926
10 5.094359e+02 3.170517e+01
* time: 0.4236130714416504
11 5.086172e+02 3.327331e+01
* time: 0.4271080493927002
12 5.080941e+02 2.942077e+01
* time: 0.4308750629425049
13 5.074009e+02 2.839941e+01
* time: 0.4348618984222412
14 5.059302e+02 3.330093e+01
* time: 0.46715402603149414
15 5.036399e+02 3.172884e+01
* time: 0.4700050354003906
16 5.017004e+02 3.160020e+01
* time: 0.47279810905456543
17 5.008553e+02 2.599524e+01
* time: 0.47548794746398926
18 5.005913e+02 2.139314e+01
* time: 0.4780139923095703
19 5.003573e+02 2.134778e+01
* time: 0.4804420471191406
20 4.997249e+02 2.069868e+01
* time: 0.4829740524291992
21 4.984453e+02 1.859010e+01
* time: 0.48569703102111816
22 4.959584e+02 2.156209e+01
* time: 0.48854994773864746
23 4.923347e+02 3.030833e+01
* time: 0.4916520118713379
24 4.906916e+02 1.652278e+01
* time: 0.49504995346069336
25 4.902955e+02 6.360800e+00
* time: 0.49878406524658203
26 4.902870e+02 7.028603e+00
* time: 0.522650957107544
27 4.902193e+02 1.176895e+00
* time: 0.525702953338623
28 4.902189e+02 1.170642e+00
* time: 0.5281078815460205
29 4.902186e+02 1.167624e+00
* time: 0.5301039218902588
30 4.902145e+02 1.110377e+00
* time: 0.5324840545654297
31 4.902079e+02 1.010507e+00
* time: 0.5348629951477051
32 4.901917e+02 9.619218e-01
* time: 0.5372788906097412
33 4.901683e+02 1.001006e+00
* time: 0.5395519733428955
34 4.901473e+02 6.138233e-01
* time: 0.5420799255371094
35 4.901412e+02 1.754342e-01
* time: 0.544605016708374
36 4.901406e+02 2.617009e-02
* time: 0.5469300746917725
37 4.901405e+02 4.585882e-03
* time: 0.5489730834960938
38 4.901405e+02 7.668184e-04
* time: 0.5512189865112305FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
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: 7.295608520507812e-5
1 9.644493e+02 5.541920e+02
* time: 0.010722875595092773
2 8.376083e+02 3.988904e+02
* time: 0.01936793327331543
3 6.503174e+02 1.495777e+02
* time: 0.02622199058532715
4 6.121241e+02 8.826583e+01
* time: 0.032891035079956055
5 5.977360e+02 9.144086e+01
* time: 0.03947901725769043
6 5.911620e+02 1.000933e+02
* time: 0.04573488235473633
7 5.858372e+02 8.423844e+01
* time: 0.05179595947265625
8 5.821934e+02 5.194402e+01
* time: 0.05798792839050293
9 5.801487e+02 3.461331e+01
* time: 0.11454200744628906
10 5.789092e+02 3.888113e+01
* time: 0.11999297142028809
11 5.780054e+02 3.556605e+01
* time: 0.12565994262695312
12 5.769455e+02 3.624436e+01
* time: 0.13124489784240723
13 5.749747e+02 4.322775e+01
* time: 0.13687586784362793
14 5.721322e+02 3.722515e+01
* time: 0.14253592491149902
15 5.695879e+02 3.401586e+01
* time: 0.14832091331481934
16 5.683277e+02 2.854997e+01
* time: 0.15430903434753418
17 5.678285e+02 2.644560e+01
* time: 0.15984582901000977
18 5.673305e+02 2.744429e+01
* time: 0.16522002220153809
19 5.662430e+02 2.793918e+01
* time: 0.1705918312072754
20 5.641877e+02 2.616169e+01
* time: 0.20295000076293945
21 5.606628e+02 2.257667e+01
* time: 0.2082669734954834
22 5.530616e+02 3.832878e+01
* time: 0.21451902389526367
23 5.528349e+02 5.518159e+01
* time: 0.2220628261566162
24 5.497231e+02 3.042064e+01
* time: 0.22947001457214355
25 5.488355e+02 6.929306e+00
* time: 0.2362828254699707
26 5.486095e+02 1.087865e+00
* time: 0.24308085441589355
27 5.486062e+02 6.456402e-01
* time: 0.2489919662475586
28 5.486061e+02 6.467689e-01
* time: 0.25444984436035156
29 5.486060e+02 6.463480e-01
* time: 0.259929895401001
30 5.486055e+02 6.408914e-01
* time: 0.2813568115234375
31 5.486045e+02 6.208208e-01
* time: 0.2867159843444824
32 5.486020e+02 1.035462e+00
* time: 0.29250192642211914
33 5.485971e+02 1.452099e+00
* time: 0.29828596115112305
34 5.485897e+02 1.482593e+00
* time: 0.30391383171081543
35 5.485839e+02 8.420646e-01
* time: 0.30998802185058594
36 5.485822e+02 2.023876e-01
* time: 0.31621789932250977
37 5.485821e+02 1.885486e-02
* time: 0.32202601432800293
38 5.485821e+02 2.343932e-03
* time: 0.3272058963775635
39 5.485821e+02 4.417566e-04
* time: 0.3317708969116211FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
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 variables: DepotFast, DepotSlow, Central
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 variables: DepotFast, DepotSlow, Central
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: 6.604194641113281e-5
1 2.618605e+02 1.558889e+02
* time: 0.1091611385345459
2 2.484481e+02 5.652093e+01
* time: 0.17518401145935059
3 2.466691e+02 3.819960e+01
* time: 0.24222707748413086
4 2.430995e+02 2.281161e+01
* time: 0.3675220012664795
5 2.421492e+02 1.230731e+01
* time: 0.41727209091186523
6 2.412891e+02 1.022578e+01
* time: 0.47058796882629395
7 2.398220e+02 4.018353e+00
* time: 0.5258369445800781
8 2.398069e+02 2.747659e+00
* time: 0.580855131149292
9 2.398011e+02 1.787831e+00
* time: 0.6340110301971436
10 2.397895e+02 2.923097e+00
* time: 0.738166093826294
11 2.397511e+02 4.766107e+00
* time: 0.7881531715393066
12 2.397274e+02 2.427088e+00
* time: 0.8417630195617676
13 2.396969e+02 2.317178e+00
* time: 0.896320104598999
14 2.396432e+02 6.390949e+00
* time: 0.9515480995178223
15 2.395437e+02 1.431395e+01
* time: 1.007277011871338
16 2.394242e+02 1.660863e+01
* time: 1.0843911170959473
17 2.392883e+02 7.918244e+00
* time: 1.1349470615386963
18 2.392648e+02 5.359505e+00
* time: 1.1981170177459717
19 2.392551e+02 8.504113e-01
* time: 1.2525131702423096
20 2.392536e+02 8.658360e-01
* time: 1.3063499927520752
21 2.392502e+02 1.359484e+00
* time: 1.3852150440216064
22 2.392455e+02 1.908268e+00
* time: 1.4358949661254883
23 2.392352e+02 3.434293e+00
* time: 1.4874939918518066
24 2.392094e+02 4.019050e+00
* time: 1.5426909923553467
25 2.391904e+02 2.360105e+00
* time: 1.6081969738006592
26 2.391639e+02 2.232579e+00
* time: 1.6986620426177979
27 2.390985e+02 7.831657e+00
* time: 1.7499029636383057
28 2.390462e+02 1.220981e+01
* time: 1.8006000518798828
29 2.389849e+02 1.372376e+01
* time: 1.8652350902557373
30 2.389295e+02 1.012359e+01
* time: 1.9208149909973145
31 2.389050e+02 7.364624e+00
* time: 1.9853839874267578
32 2.388675e+02 3.493348e+00
* time: 2.0625741481781006
33 2.388348e+02 3.602127e+00
* time: 2.112338066101074
34 2.388092e+02 9.896348e+00
* time: 2.166102170944214
35 2.387816e+02 5.815979e+00
* time: 2.2213401794433594
36 2.387667e+02 3.995259e+00
* time: 2.276092052459717
37 2.387624e+02 2.110864e+00
* time: 2.3313679695129395
38 2.387542e+02 3.340382e+00
* time: 2.407546043395996
39 2.387468e+02 3.743004e+00
* time: 2.4566750526428223
40 2.387148e+02 8.270607e+00
* time: 2.5101630687713623
41 2.386920e+02 8.043639e+00
* time: 2.564116954803467
42 2.386311e+02 1.788043e+00
* time: 2.6195180416107178
43 2.385981e+02 4.286869e+00
* time: 2.6741011142730713
44 2.385768e+02 2.412656e+00
* time: 2.7502100467681885
45 2.385642e+02 4.499325e-01
* time: 2.8002090454101562
46 2.385598e+02 1.789469e+00
* time: 2.852627992630005
47 2.385584e+02 4.262494e-01
* time: 2.9064149856567383
48 2.385582e+02 3.514292e-01
* time: 2.9596941471099854
49 2.385581e+02 9.399545e-02
* time: 3.0143351554870605
50 2.385580e+02 1.262250e-02
* time: 3.103199005126953
51 2.385580e+02 2.333105e-03
* time: 3.150290012359619
52 2.385580e+02 2.333105e-03
* time: 3.213879108428955
53 2.385580e+02 2.333105e-03
* time: 3.2365920543670654FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
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: 7.414817810058594e-5
1 2.991160e+02 1.342892e+02
* time: 0.19059014320373535
2 2.851580e+02 4.511285e+01
* time: 0.25725817680358887
3 2.835992e+02 4.343268e+01
* time: 0.3268871307373047
4 2.821916e+02 4.373720e+01
* time: 0.38590407371520996
5 2.794524e+02 1.253857e+01
* time: 0.4928750991821289
6 2.787961e+02 1.102093e+01
* time: 0.5464382171630859
7 2.770092e+02 1.114239e+01
* time: 0.6014261245727539
8 2.769202e+02 5.089120e+00
* time: 0.6556169986724854
9 2.769088e+02 3.041932e+00
* time: 0.711414098739624
10 2.768987e+02 1.603805e+00
* time: 0.7683110237121582
11 2.767766e+02 1.411813e+01
* time: 0.8573930263519287
12 2.766170e+02 2.719462e+01
* time: 0.9122650623321533
13 2.764580e+02 3.014175e+01
* time: 0.9667601585388184
14 2.762515e+02 1.669154e+01
* time: 1.0208611488342285
15 2.761434e+02 1.957266e+00
* time: 1.0766830444335938
16 2.761353e+02 3.756416e+00
* time: 1.152433156967163
17 2.761237e+02 1.341505e+00
* time: 1.2063422203063965
18 2.761152e+02 2.182776e+00
* time: 1.2606310844421387
19 2.761036e+02 3.209868e+00
* time: 1.314601182937622
20 2.760912e+02 3.070994e+00
* time: 1.3681871891021729
21 2.760762e+02 1.990495e+00
* time: 1.4219532012939453
22 2.760587e+02 1.586630e+00
* time: 1.4944031238555908
23 2.760290e+02 2.471581e+00
* time: 1.549814224243164
24 2.759737e+02 1.040810e+01
* time: 1.605112075805664
25 2.758948e+02 9.671069e+00
* time: 1.6590301990509033
26 2.758264e+02 7.208229e+00
* time: 1.7137031555175781
27 2.757728e+02 3.300150e+00
* time: 1.794323205947876
28 2.757344e+02 4.723361e+00
* time: 1.8482511043548584
29 2.756731e+02 6.515164e+00
* time: 1.902961015701294
30 2.755937e+02 4.176876e+00
* time: 1.957705020904541
31 2.755684e+02 5.161233e+00
* time: 2.021742105484009
32 2.755453e+02 8.103754e+00
* time: 2.0933451652526855
33 2.755346e+02 2.183052e+00
* time: 2.1484410762786865
34 2.755279e+02 3.129360e+00
* time: 2.2026331424713135
35 2.755221e+02 3.813180e+00
* time: 2.2554261684417725
36 2.754865e+02 6.625383e+00
* time: 2.3092920780181885
37 2.754662e+02 5.264434e+00
* time: 2.3629980087280273
38 2.754486e+02 9.924709e-01
* time: 2.4353370666503906
39 2.754420e+02 1.102256e+00
* time: 2.4895992279052734
40 2.754395e+02 3.092900e-01
* time: 2.5439600944519043
41 2.754386e+02 5.369253e-01
* time: 2.5979180335998535
42 2.754384e+02 1.423783e-01
* time: 2.6498842239379883
43 2.754383e+02 1.288451e-01
* time: 2.7033281326293945
44 2.754381e+02 3.887239e-01
* time: 2.773451089859009
45 2.754376e+02 9.266434e-01
* time: 2.8265621662139893
46 2.754363e+02 1.786193e+00
* time: 2.880441188812256
47 2.754331e+02 3.042628e+00
* time: 2.9338772296905518
48 2.754256e+02 4.634171e+00
* time: 2.9869232177734375
49 2.754117e+02 5.954940e+00
* time: 3.0589351654052734
50 2.753920e+02 6.732566e+00
* time: 3.113706111907959
51 2.753738e+02 2.208977e+00
* time: 3.1682941913604736
52 2.753698e+02 2.481125e+00
* time: 3.2225260734558105
53 2.753538e+02 9.838345e-01
* time: 3.2760350704193115
54 2.753535e+02 3.815520e+00
* time: 3.331543207168579
55 2.753483e+02 1.735894e+00
* time: 3.4025251865386963
56 2.753436e+02 8.278884e-01
* time: 3.456317186355591
57 2.753339e+02 1.256117e+00
* time: 3.5102531909942627
58 2.753262e+02 2.343600e+00
* time: 3.5647661685943604
59 2.753231e+02 1.270829e+00
* time: 3.6176891326904297
60 2.753207e+02 3.103092e-01
* time: 3.6710050106048584
61 2.753202e+02 3.488786e-01
* time: 3.740567207336426
62 2.753193e+02 4.518882e-01
* time: 3.793715000152588
63 2.753186e+02 2.027309e-01
* time: 3.8470911979675293
64 2.753183e+02 9.363635e-02
* time: 3.8996801376342773
65 2.753182e+02 9.589112e-02
* time: 3.9510531425476074
66 2.753181e+02 9.484170e-02
* time: 4.00370717048645
67 2.753181e+02 4.831472e-02
* time: 4.0727081298828125
68 2.753180e+02 2.328249e-02
* time: 4.1244590282440186
69 2.753180e+02 1.740503e-02
* time: 4.174759149551392
70 2.753180e+02 1.793621e-02
* time: 4.2256669998168945
71 2.753180e+02 8.549994e-03
* time: 4.275990009307861
72 2.753180e+02 5.684889e-03
* time: 4.327735185623169
73 2.753180e+02 4.082116e-03
* time: 4.399227142333984
74 2.753180e+02 7.506601e-03
* time: 4.449720144271851
75 2.753180e+02 6.945870e-03
* time: 4.499985218048096
76 2.753180e+02 5.205184e-03
* time: 4.554447174072266
77 2.753180e+02 2.313568e-03
* time: 4.603983163833618
78 2.753180e+02 2.313568e-03
* time: 4.66821813583374FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Log-likelihood value: -275.31801
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.13275
ωCL 0.081471
ωVc 0.10178
nbioav 1.3974e8
σ 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.9882 |
| 2 | θkaSlow | 0.607708 | 0.613328 |
| 3 | θCL | 1.07968 | 1.07967 |
| 4 | θVc | 11.3745 | 11.381 |
| 5 | θbioav | 0.143942 | 0.132746 |
| 6 | ωCL | 0.0815101 | 0.0814711 |
| 7 | ωVc | 0.101565 | 0.101784 |
| 8 | σ | 0.10441 | 0.104644 |
| 9 | ξbioav | 0.125128 | missing |
| 10 | nbioav | missing | 1.39739e8 |
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.54e-269 | 0.0 |
| 4 | 0.003003 | 4.48998e-222 | 0.0 |
| 5 | 0.004004 | 3.80318e-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.