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.02042984962463379
1 9.433464e+02 6.079483e+02
* time: 0.8397860527038574
2 8.189627e+02 4.423725e+02
* time: 0.8657889366149902
3 5.917683e+02 1.819248e+02
* time: 0.8857529163360596
4 5.421783e+02 1.121313e+02
* time: 0.9018678665161133
5 5.255651e+02 7.407230e+01
* time: 0.9165918827056885
6 5.208427e+02 8.699271e+01
* time: 0.9299800395965576
7 5.174883e+02 8.974584e+01
* time: 0.9430439472198486
8 5.138523e+02 7.328235e+01
* time: 0.9557769298553467
9 5.109883e+02 4.155805e+01
* time: 0.9680249691009521
10 5.094359e+02 3.170517e+01
* time: 1.005375862121582
11 5.086172e+02 3.327331e+01
* time: 1.0168569087982178
12 5.080941e+02 2.942077e+01
* time: 1.0283479690551758
13 5.074009e+02 2.839941e+01
* time: 1.0397119522094727
14 5.059302e+02 3.330093e+01
* time: 1.0505728721618652
15 5.036399e+02 3.172884e+01
* time: 1.0623178482055664
16 5.017004e+02 3.160020e+01
* time: 1.0745270252227783
17 5.008553e+02 2.599524e+01
* time: 1.0867209434509277
18 5.005913e+02 2.139314e+01
* time: 1.098296880722046
19 5.003573e+02 2.134778e+01
* time: 1.1096160411834717
20 4.997249e+02 2.069868e+01
* time: 1.121438980102539
21 4.984453e+02 1.859010e+01
* time: 1.1334278583526611
22 4.959584e+02 2.156209e+01
* time: 1.1456449031829834
23 4.923347e+02 3.030833e+01
* time: 1.158552885055542
24 4.906916e+02 1.652278e+01
* time: 1.1724109649658203
25 4.902955e+02 6.360800e+00
* time: 1.1852660179138184
26 4.902870e+02 7.028603e+00
* time: 1.1991069316864014
27 4.902193e+02 1.176895e+00
* time: 1.2316389083862305
28 4.902189e+02 1.170642e+00
* time: 1.2420358657836914
29 4.902186e+02 1.167624e+00
* time: 1.2509489059448242
30 4.902145e+02 1.110377e+00
* time: 1.2619140148162842
31 4.902079e+02 1.010507e+00
* time: 1.2728948593139648
32 4.901917e+02 9.619218e-01
* time: 1.2841930389404297
33 4.901683e+02 1.001006e+00
* time: 1.2943298816680908
34 4.901473e+02 6.138233e-01
* time: 1.3071980476379395
35 4.901412e+02 1.754342e-01
* time: 1.3189139366149902
36 4.901406e+02 2.617009e-02
* time: 1.329360008239746
37 4.901405e+02 4.585882e-03
* time: 1.3385679721832275
38 4.901405e+02 7.668184e-04
* time: 1.3461809158325195FittedPumasModel
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: 5.602836608886719e-5
1 9.056373e+02 5.541920e+02
* time: 0.02894306182861328
2 7.787963e+02 3.988904e+02
* time: 0.1283550262451172
3 5.915054e+02 1.495777e+02
* time: 0.1464521884918213
4 5.533120e+02 8.826583e+01
* time: 0.16358208656311035
5 5.389239e+02 9.144086e+01
* time: 0.17986607551574707
6 5.323499e+02 1.000933e+02
* time: 0.1958010196685791
7 5.270252e+02 8.423844e+01
* time: 0.2110581398010254
8 5.233813e+02 5.194402e+01
* time: 0.22545504570007324
9 5.213366e+02 3.461331e+01
* time: 0.2401599884033203
10 5.200972e+02 3.888113e+01
* time: 0.2547750473022461
11 5.191933e+02 3.556605e+01
* time: 0.2692751884460449
12 5.181335e+02 3.624436e+01
* time: 0.28419017791748047
13 5.161626e+02 4.322775e+01
* time: 0.29912614822387695
14 5.133202e+02 3.722515e+01
* time: 0.31336212158203125
15 5.107758e+02 3.401586e+01
* time: 0.3279430866241455
16 5.095157e+02 2.854997e+01
* time: 0.34288597106933594
17 5.090165e+02 2.644560e+01
* time: 0.3566451072692871
18 5.085184e+02 2.744429e+01
* time: 0.36995506286621094
19 5.074309e+02 2.793918e+01
* time: 0.3829841613769531
20 5.053757e+02 2.616169e+01
* time: 0.39652299880981445
21 5.018507e+02 2.257667e+01
* time: 0.45429515838623047
22 4.942495e+02 3.832878e+01
* time: 0.4710509777069092
23 4.940229e+02 5.518159e+01
* time: 0.49210619926452637
24 4.909110e+02 3.042064e+01
* time: 0.512160062789917
25 4.900234e+02 6.929306e+00
* time: 0.5299201011657715
26 4.897974e+02 1.087865e+00
* time: 0.5490851402282715
27 4.897942e+02 6.456402e-01
* time: 0.5645411014556885
28 4.897940e+02 6.467689e-01
* time: 0.5785460472106934
29 4.897939e+02 6.463480e-01
* time: 0.5926971435546875
30 4.897935e+02 6.408914e-01
* time: 0.607105016708374
31 4.897924e+02 6.208208e-01
* time: 0.6217780113220215
32 4.897900e+02 1.035462e+00
* time: 0.6368839740753174
33 4.897850e+02 1.452099e+00
* time: 0.6521401405334473
34 4.897776e+02 1.482593e+00
* time: 0.6705591678619385
35 4.897718e+02 8.420646e-01
* time: 0.6862890720367432
36 4.897702e+02 2.023876e-01
* time: 0.7013580799102783
37 4.897700e+02 1.885486e-02
* time: 0.7151050567626953
38 4.897700e+02 2.343932e-03
* time: 0.7274811267852783
39 4.897700e+02 4.417566e-04
* time: 0.7507190704345703FittedPumasModel
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):
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 2.512964e+02 3.650794e+02
* time: 6.580352783203125e-5
1 2.005689e+02 1.252637e+02
* time: 0.27806782722473145
2 1.875728e+02 4.225518e+01
* time: 0.40436387062072754
3 1.863803e+02 3.429785e+01
* time: 0.5180678367614746
4 1.845581e+02 3.694229e+01
* time: 0.6123218536376953
5 1.828416e+02 1.170915e+01
* time: 0.7706289291381836
6 1.823277e+02 1.004017e+01
* time: 0.8625869750976562
7 1.810629e+02 5.326780e+00
* time: 0.9559628963470459
8 1.810479e+02 1.767987e+00
* time: 1.0455639362335205
9 1.810272e+02 3.852037e+00
* time: 1.140915870666504
10 1.809414e+02 1.196625e+01
* time: 1.2429289817810059
11 1.806489e+02 1.861440e+01
* time: 1.348660945892334
12 1.801984e+02 1.361774e+01
* time: 1.764153003692627
13 1.800427e+02 2.177039e+01
* time: 1.8610119819641113
14 1.796554e+02 7.079039e+00
* time: 1.9564948081970215
15 1.795832e+02 1.581499e+01
* time: 2.060091972351074
16 1.795220e+02 6.552120e+00
* time: 2.164411783218384
17 1.794621e+02 6.612645e+00
* time: 2.2678539752960205
18 1.793643e+02 6.603903e+00
* time: 2.378734827041626
19 1.793502e+02 1.025787e+01
* time: 2.513934850692749
20 1.793178e+02 1.202199e+01
* time: 2.6095688343048096
21 1.792424e+02 1.625661e+01
* time: 2.7060768604278564
22 1.791560e+02 1.128856e+01
* time: 2.8015809059143066
23 1.790991e+02 7.625637e+00
* time: 2.9023327827453613
24 1.790756e+02 8.557258e+00
* time: 3.00312876701355
25 1.790633e+02 4.848482e+00
* time: 3.1059207916259766
26 1.790408e+02 5.859306e+00
* time: 3.207359790802002
27 1.789734e+02 1.106035e+01
* time: 3.3109359741210938
28 1.789070e+02 1.143361e+01
* time: 3.434027910232544
29 1.788321e+02 7.098448e+00
* time: 3.529693841934204
30 1.787896e+02 5.709857e+00
* time: 3.624307870864868
31 1.787809e+02 7.456555e+00
* time: 3.7208139896392822
32 1.787671e+02 7.320931e-01
* time: 3.82366681098938
33 1.787660e+02 5.023990e-01
* time: 3.926435947418213
34 1.787656e+02 3.813994e-01
* time: 4.0245139598846436
35 1.787639e+02 8.608909e-01
* time: 4.123638868331909
36 1.787607e+02 2.047321e+00
* time: 4.248727798461914
37 1.787525e+02 3.859529e+00
* time: 4.338005781173706
38 1.787349e+02 5.864920e+00
* time: 4.429638862609863
39 1.787017e+02 6.966045e+00
* time: 4.5223047733306885
40 1.786574e+02 5.348939e+00
* time: 4.622135877609253
41 1.786270e+02 1.642025e+00
* time: 4.722083806991577
42 1.786181e+02 5.211823e-01
* time: 4.820537805557251
43 1.786153e+02 1.186061e+00
* time: 4.919969797134399
44 1.786125e+02 1.292005e+00
* time: 5.015107870101929
45 1.786099e+02 7.814598e-01
* time: 5.131313800811768
46 1.786086e+02 1.369456e-01
* time: 5.221433877944946
47 1.786082e+02 1.912170e-01
* time: 5.3089189529418945
48 1.786080e+02 2.670802e-01
* time: 5.3984479904174805
49 1.786078e+02 1.979262e-01
* time: 5.493316888809204
50 1.786077e+02 5.177918e-02
* time: 5.58785080909729
51 1.786076e+02 2.998328e-02
* time: 5.683206796646118
52 1.786076e+02 5.799706e-02
* time: 5.7755959033966064
53 1.786076e+02 4.280434e-02
* time: 5.869516849517822
54 1.786076e+02 1.467829e-02
* time: 5.986361980438232
55 1.786076e+02 6.928178e-03
* time: 6.066188812255859
56 1.786076e+02 1.236015e-02
* time: 6.150690793991089
57 1.786076e+02 9.950826e-03
* time: 6.234149932861328
58 1.786076e+02 2.974370e-03
* time: 6.322774887084961
59 1.786076e+02 1.374576e-03
* time: 6.414756774902344
60 1.786076e+02 2.860078e-03
* time: 6.4990479946136475
61 1.786076e+02 2.100127e-03
* time: 6.58353590965271
62 1.786076e+02 2.100127e-03
* time: 6.7107977867126465
63 1.786076e+02 6.412834e-03
* time: 6.834594964981079
64 1.786076e+02 6.412834e-03
* time: 6.9713218212127686
65 1.786076e+02 6.412834e-03
* time: 7.156803846359253
66 1.786076e+02 6.412834e-03
* time: 7.444414854049683
67 1.786076e+02 6.412834e-03
* time: 7.665327787399292FittedPumasModel
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: 5.316734313964844e-5
1 2.014577e+02 1.265001e+02
* time: 0.2949249744415283
2 1.880885e+02 4.190708e+01
* time: 0.44132208824157715
3 1.870317e+02 8.825666e+01
* time: 0.5463001728057861
4 1.846027e+02 4.400156e+01
* time: 0.6741721630096436
5 1.834445e+02 1.906624e+01
* time: 0.8503191471099854
6 1.828599e+02 1.113882e+01
* time: 0.9424631595611572
7 1.815719e+02 7.449355e+00
* time: 1.0367910861968994
8 1.815131e+02 2.164678e+00
* time: 1.1308059692382812
9 1.814896e+02 2.167319e+00
* time: 1.2327191829681396
10 1.814458e+02 4.615738e+00
* time: 1.3339321613311768
11 1.813173e+02 9.576967e+00
* time: 1.4360170364379883
12 1.809756e+02 2.052077e+01
* time: 1.5425829887390137
13 1.807625e+02 4.553366e+01
* time: 1.6894891262054443
14 1.802224e+02 6.550892e+00
* time: 1.7825870513916016
15 1.800862e+02 2.865509e+00
* time: 1.8777351379394531
16 1.800780e+02 1.164611e+00
* time: 1.9698760509490967
17 1.800737e+02 7.952462e-01
* time: 2.065891981124878
18 1.800352e+02 4.860618e+00
* time: 2.171504020690918
19 1.800089e+02 5.176689e+00
* time: 2.2761080265045166
20 1.799679e+02 4.303892e+00
* time: 2.381277084350586
21 1.799153e+02 4.612832e+00
* time: 2.487298011779785
22 1.798423e+02 1.209387e+01
* time: 2.612706184387207
23 1.796821e+02 1.712256e+01
* time: 2.709073066711426
24 1.794275e+02 1.435100e+01
* time: 2.8036751747131348
25 1.793773e+02 4.137313e+00
* time: 2.8986880779266357
26 1.793488e+02 1.846545e+00
* time: 3.0291590690612793
27 1.793331e+02 5.502533e+00
* time: 3.1314611434936523
28 1.793234e+02 2.894037e+00
* time: 3.2332980632781982
29 1.793119e+02 1.453372e+00
* time: 3.3321101665496826
30 1.792879e+02 4.109884e+00
* time: 3.4589531421661377
31 1.792599e+02 4.613173e+00
* time: 3.551457166671753
32 1.792384e+02 4.254549e+00
* time: 3.643778085708618
33 1.792251e+02 4.415024e+00
* time: 3.7352311611175537
34 1.792026e+02 3.229145e+00
* time: 3.8345530033111572
35 1.791841e+02 3.422094e+00
* time: 3.9347541332244873
36 1.791705e+02 1.621228e+00
* time: 4.037853002548218
37 1.791555e+02 3.462667e+00
* time: 4.135825157165527
38 1.791293e+02 5.586685e+00
* time: 4.236175060272217
39 1.790713e+02 9.927638e+00
* time: 4.3542890548706055
40 1.789891e+02 1.133271e+01
* time: 4.447052001953125
41 1.788585e+02 1.279172e+01
* time: 4.536435127258301
42 1.787407e+02 5.291681e+00
* time: 4.62983512878418
43 1.786633e+02 5.367971e+00
* time: 4.733007192611694
44 1.786405e+02 2.571548e+00
* time: 4.830978155136108
45 1.786339e+02 2.720314e+00
* time: 4.93030309677124
46 1.786252e+02 1.930583e+00
* time: 5.029952049255371
47 1.786182e+02 1.523164e+00
* time: 5.155141115188599
48 1.786126e+02 5.027920e-01
* time: 5.245054006576538
49 1.786100e+02 3.733023e-01
* time: 5.332525014877319
50 1.786089e+02 2.749553e-01
* time: 5.42055606842041
51 1.786083e+02 1.981772e-01
* time: 5.508686065673828
52 1.786079e+02 1.005262e-01
* time: 5.6030189990997314
53 1.786078e+02 2.549641e-02
* time: 5.698633193969727
54 1.786077e+02 4.452931e-02
* time: 5.790536165237427
55 1.786076e+02 2.967125e-02
* time: 5.883052110671997
56 1.786076e+02 1.068274e-02
* time: 5.975605010986328
57 1.786076e+02 5.355447e-03
* time: 6.0869550704956055
58 1.786076e+02 8.180920e-03
* time: 6.170389175415039
59 1.786076e+02 4.935439e-03
* time: 6.251726150512695
60 1.786076e+02 4.387986e-03
* time: 6.351366996765137
61 1.786076e+02 4.388023e-03
* time: 6.496145009994507
62 1.786076e+02 4.387934e-03
* time: 6.625404119491577
63 1.786076e+02 4.387934e-03
* time: 6.733436107635498FittedPumasModel
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.