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.02231884002685547
1 9.433464e+02 6.079483e+02
* time: 0.9089598655700684
2 8.189627e+02 4.423725e+02
* time: 0.9359259605407715
3 5.917683e+02 1.819248e+02
* time: 0.9564468860626221
4 5.421783e+02 1.121313e+02
* time: 0.973153829574585
5 5.255651e+02 7.407230e+01
* time: 0.987724781036377
6 5.208427e+02 8.699271e+01
* time: 1.0012898445129395
7 5.174883e+02 8.974584e+01
* time: 1.0145249366760254
8 5.138523e+02 7.328235e+01
* time: 1.028228998184204
9 5.109883e+02 4.155805e+01
* time: 1.041111946105957
10 5.094359e+02 3.170517e+01
* time: 1.0528509616851807
11 5.086172e+02 3.327331e+01
* time: 1.0648598670959473
12 5.080941e+02 2.942077e+01
* time: 1.1041529178619385
13 5.074009e+02 2.839941e+01
* time: 1.1163158416748047
14 5.059302e+02 3.330093e+01
* time: 1.127551794052124
15 5.036399e+02 3.172884e+01
* time: 1.1397809982299805
16 5.017004e+02 3.160020e+01
* time: 1.1522009372711182
17 5.008553e+02 2.599524e+01
* time: 1.166186809539795
18 5.005913e+02 2.139314e+01
* time: 1.1781549453735352
19 5.003573e+02 2.134778e+01
* time: 1.1899518966674805
20 4.997249e+02 2.069868e+01
* time: 1.2019009590148926
21 4.984453e+02 1.859010e+01
* time: 1.213824987411499
22 4.959584e+02 2.156209e+01
* time: 1.2262349128723145
23 4.923347e+02 3.030833e+01
* time: 1.2404627799987793
24 4.906916e+02 1.652278e+01
* time: 1.2547249794006348
25 4.902955e+02 6.360800e+00
* time: 1.2680227756500244
26 4.902870e+02 7.028603e+00
* time: 1.282315969467163
27 4.902193e+02 1.176895e+00
* time: 1.295658826828003
28 4.902189e+02 1.170642e+00
* time: 1.3066627979278564
29 4.902186e+02 1.167624e+00
* time: 1.3348939418792725
30 4.902145e+02 1.110377e+00
* time: 1.3460578918457031
31 4.902079e+02 1.010507e+00
* time: 1.3576288223266602
32 4.901917e+02 9.619218e-01
* time: 1.3692197799682617
33 4.901683e+02 1.001006e+00
* time: 1.3800227642059326
34 4.901473e+02 6.138233e-01
* time: 1.3919439315795898
35 4.901412e+02 1.754342e-01
* time: 1.4039578437805176
36 4.901406e+02 2.617009e-02
* time: 1.4149038791656494
37 4.901405e+02 4.585882e-03
* time: 1.4243087768554688
38 4.901405e+02 7.668184e-04
* time: 1.4325029850006104FittedPumasModel
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.91278076171875e-5
1 9.056373e+02 5.541920e+02
* time: 0.032911062240600586
2 7.787963e+02 3.988904e+02
* time: 0.05771303176879883
3 5.915054e+02 1.495777e+02
* time: 0.07742905616760254
4 5.533120e+02 8.826583e+01
* time: 0.14090895652770996
5 5.389239e+02 9.144086e+01
* time: 0.15808606147766113
6 5.323499e+02 1.000933e+02
* time: 0.17420601844787598
7 5.270252e+02 8.423844e+01
* time: 0.19119000434875488
8 5.233813e+02 5.194402e+01
* time: 0.2059791088104248
9 5.213366e+02 3.461331e+01
* time: 0.22099709510803223
10 5.200972e+02 3.888113e+01
* time: 0.23606300354003906
11 5.191933e+02 3.556605e+01
* time: 0.251162052154541
12 5.181335e+02 3.624436e+01
* time: 0.26656603813171387
13 5.161626e+02 4.322775e+01
* time: 0.2820441722869873
14 5.133202e+02 3.722515e+01
* time: 0.2968721389770508
15 5.107758e+02 3.401586e+01
* time: 0.31187915802001953
16 5.095157e+02 2.854997e+01
* time: 0.32714295387268066
17 5.090165e+02 2.644560e+01
* time: 0.3416321277618408
18 5.085184e+02 2.744429e+01
* time: 0.3553340435028076
19 5.074309e+02 2.793918e+01
* time: 0.3684101104736328
20 5.053757e+02 2.616169e+01
* time: 0.3820769786834717
21 5.018507e+02 2.257667e+01
* time: 0.39679408073425293
22 4.942495e+02 3.832878e+01
* time: 0.41304898262023926
23 4.940229e+02 5.518159e+01
* time: 0.43338799476623535
24 4.909110e+02 3.042064e+01
* time: 0.49193811416625977
25 4.900234e+02 6.929306e+00
* time: 0.5105400085449219
26 4.897974e+02 1.087865e+00
* time: 0.5286710262298584
27 4.897942e+02 6.456402e-01
* time: 0.5451710224151611
28 4.897940e+02 6.467689e-01
* time: 0.5603351593017578
29 4.897939e+02 6.463480e-01
* time: 0.5756680965423584
30 4.897935e+02 6.408914e-01
* time: 0.591055154800415
31 4.897924e+02 6.208208e-01
* time: 0.6068849563598633
32 4.897900e+02 1.035462e+00
* time: 0.6227350234985352
33 4.897850e+02 1.452099e+00
* time: 0.638700008392334
34 4.897776e+02 1.482593e+00
* time: 0.6544599533081055
35 4.897718e+02 8.420646e-01
* time: 0.6710870265960693
36 4.897702e+02 2.023876e-01
* time: 0.6890699863433838
37 4.897700e+02 1.885486e-02
* time: 0.704380989074707
38 4.897700e+02 2.343932e-03
* time: 0.7179141044616699
39 4.897700e+02 4.417566e-04
* time: 0.7296750545501709FittedPumasModel
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.605552673339844e-5
1 2.005689e+02 1.252637e+02
* time: 0.32099008560180664
2 1.875728e+02 4.225518e+01
* time: 0.465972900390625
3 1.863803e+02 3.429785e+01
* time: 0.5986950397491455
4 1.845581e+02 3.694229e+01
* time: 0.7933950424194336
5 1.828416e+02 1.170915e+01
* time: 0.8881440162658691
6 1.823277e+02 1.004017e+01
* time: 0.9829208850860596
7 1.810629e+02 5.326780e+00
* time: 1.081355094909668
8 1.810479e+02 1.767987e+00
* time: 1.1842050552368164
9 1.810272e+02 3.852037e+00
* time: 1.2874798774719238
10 1.809414e+02 1.196625e+01
* time: 1.3936851024627686
11 1.806489e+02 1.861440e+01
* time: 1.5028009414672852
12 1.801984e+02 1.361774e+01
* time: 1.9089529514312744
13 1.800427e+02 2.177039e+01
* time: 2.0121569633483887
14 1.796554e+02 7.079039e+00
* time: 2.122056007385254
15 1.795832e+02 1.581499e+01
* time: 2.230936050415039
16 1.795220e+02 6.552120e+00
* time: 2.338670015335083
17 1.794621e+02 6.612645e+00
* time: 2.4450759887695312
18 1.793643e+02 6.603903e+00
* time: 2.581969976425171
19 1.793502e+02 1.025787e+01
* time: 2.6800789833068848
20 1.793178e+02 1.202199e+01
* time: 2.7783150672912598
21 1.792424e+02 1.625661e+01
* time: 2.8785970211029053
22 1.791560e+02 1.128856e+01
* time: 2.9864299297332764
23 1.790991e+02 7.625637e+00
* time: 3.091507911682129
24 1.790756e+02 8.557258e+00
* time: 3.197374105453491
25 1.790633e+02 4.848482e+00
* time: 3.30410099029541
26 1.790408e+02 5.859306e+00
* time: 3.434760093688965
27 1.789734e+02 1.106035e+01
* time: 3.5329620838165283
28 1.789070e+02 1.143361e+01
* time: 3.6306231021881104
29 1.788321e+02 7.098448e+00
* time: 3.7301080226898193
30 1.787896e+02 5.709857e+00
* time: 3.8344199657440186
31 1.787809e+02 7.456555e+00
* time: 3.9404289722442627
32 1.787671e+02 7.320931e-01
* time: 4.04642391204834
33 1.787660e+02 5.023990e-01
* time: 4.151829957962036
34 1.787656e+02 3.813994e-01
* time: 4.254441976547241
35 1.787639e+02 8.608909e-01
* time: 4.380604028701782
36 1.787607e+02 2.047321e+00
* time: 4.472848892211914
37 1.787525e+02 3.859529e+00
* time: 4.563884019851685
38 1.787349e+02 5.864920e+00
* time: 4.658329963684082
39 1.787017e+02 6.966045e+00
* time: 4.763763904571533
40 1.786574e+02 5.348939e+00
* time: 4.867713928222656
41 1.786270e+02 1.642025e+00
* time: 4.970438003540039
42 1.786181e+02 5.211823e-01
* time: 5.0714030265808105
43 1.786153e+02 1.186061e+00
* time: 5.198873996734619
44 1.786125e+02 1.292005e+00
* time: 5.289319038391113
45 1.786099e+02 7.814598e-01
* time: 5.3817009925842285
46 1.786086e+02 1.369456e-01
* time: 5.474068880081177
47 1.786082e+02 1.912170e-01
* time: 5.5681750774383545
48 1.786080e+02 2.670802e-01
* time: 5.6691389083862305
49 1.786078e+02 1.979262e-01
* time: 5.76788592338562
50 1.786077e+02 5.177918e-02
* time: 5.86606502532959
51 1.786076e+02 2.998328e-02
* time: 5.965444087982178
52 1.786076e+02 5.799706e-02
* time: 6.090219020843506
53 1.786076e+02 4.280434e-02
* time: 6.1806581020355225
54 1.786076e+02 1.467829e-02
* time: 6.27051305770874
55 1.786076e+02 6.928178e-03
* time: 6.353360891342163
56 1.786076e+02 1.236015e-02
* time: 6.442904949188232
57 1.786076e+02 9.950826e-03
* time: 6.536591053009033
58 1.786076e+02 2.974370e-03
* time: 6.627985954284668
59 1.786076e+02 1.374576e-03
* time: 6.7227630615234375
60 1.786076e+02 2.860078e-03
* time: 6.810518980026245
61 1.786076e+02 2.100127e-03
* time: 6.8981099128723145
62 1.786076e+02 2.100127e-03
* time: 7.047461986541748
63 1.786076e+02 6.412834e-03
* time: 7.146109104156494
64 1.786076e+02 6.412834e-03
* time: 7.288794040679932
65 1.786076e+02 6.412834e-03
* time: 7.4922709465026855
66 1.786076e+02 6.412834e-03
* time: 7.819143056869507
67 1.786076e+02 6.412834e-03
* time: 8.001981973648071FittedPumasModel
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.104873657226562e-5
1 2.014577e+02 1.265001e+02
* time: 0.32230496406555176
2 1.880885e+02 4.190708e+01
* time: 0.46334004402160645
3 1.870317e+02 8.825666e+01
* time: 0.5654659271240234
4 1.846027e+02 4.400156e+01
* time: 0.7401859760284424
5 1.834445e+02 1.906624e+01
* time: 0.8332438468933105
6 1.828599e+02 1.113882e+01
* time: 0.9253349304199219
7 1.815719e+02 7.449355e+00
* time: 1.021101951599121
8 1.815131e+02 2.164678e+00
* time: 1.1221299171447754
9 1.814896e+02 2.167319e+00
* time: 1.2270328998565674
10 1.814458e+02 4.615738e+00
* time: 1.3274109363555908
11 1.813173e+02 9.576967e+00
* time: 1.4305899143218994
12 1.809756e+02 2.052077e+01
* time: 1.599350929260254
13 1.807625e+02 4.553366e+01
* time: 1.6952738761901855
14 1.802224e+02 6.550892e+00
* time: 1.7876758575439453
15 1.800862e+02 2.865509e+00
* time: 1.8818058967590332
16 1.800780e+02 1.164611e+00
* time: 1.9791738986968994
17 1.800737e+02 7.952462e-01
* time: 2.0778019428253174
18 1.800352e+02 4.860618e+00
* time: 2.181325912475586
19 1.800089e+02 5.176689e+00
* time: 2.2830710411071777
20 1.799679e+02 4.303892e+00
* time: 2.3859329223632812
21 1.799153e+02 4.612832e+00
* time: 2.513260841369629
22 1.798423e+02 1.209387e+01
* time: 2.6092779636383057
23 1.796821e+02 1.712256e+01
* time: 2.7037599086761475
24 1.794275e+02 1.435100e+01
* time: 2.7998239994049072
25 1.793773e+02 4.137313e+00
* time: 2.901576042175293
26 1.793488e+02 1.846545e+00
* time: 3.0334458351135254
27 1.793331e+02 5.502533e+00
* time: 3.136209011077881
28 1.793234e+02 2.894037e+00
* time: 3.236629009246826
29 1.793119e+02 1.453372e+00
* time: 3.3619890213012695
30 1.792879e+02 4.109884e+00
* time: 3.455335855484009
31 1.792599e+02 4.613173e+00
* time: 3.550433874130249
32 1.792384e+02 4.254549e+00
* time: 3.6455299854278564
33 1.792251e+02 4.415024e+00
* time: 3.7456319332122803
34 1.792026e+02 3.229145e+00
* time: 3.850181818008423
35 1.791841e+02 3.422094e+00
* time: 3.954767942428589
36 1.791705e+02 1.621228e+00
* time: 4.061128854751587
37 1.791555e+02 3.462667e+00
* time: 4.164203882217407
38 1.791293e+02 5.586685e+00
* time: 4.289515972137451
39 1.790713e+02 9.927638e+00
* time: 4.3800249099731445
40 1.789891e+02 1.133271e+01
* time: 4.471454858779907
41 1.788585e+02 1.279172e+01
* time: 4.5610198974609375
42 1.787407e+02 5.291681e+00
* time: 4.667816877365112
43 1.786633e+02 5.367971e+00
* time: 4.773988962173462
44 1.786405e+02 2.571548e+00
* time: 4.877553939819336
45 1.786339e+02 2.720314e+00
* time: 4.987169027328491
46 1.786252e+02 1.930583e+00
* time: 5.095970869064331
47 1.786182e+02 1.523164e+00
* time: 5.224249839782715
48 1.786126e+02 5.027920e-01
* time: 5.319468975067139
49 1.786100e+02 3.733023e-01
* time: 5.415748834609985
50 1.786089e+02 2.749553e-01
* time: 5.510723829269409
51 1.786083e+02 1.981772e-01
* time: 5.613173007965088
52 1.786079e+02 1.005262e-01
* time: 5.711282968521118
53 1.786078e+02 2.549641e-02
* time: 5.808081865310669
54 1.786077e+02 4.452931e-02
* time: 5.906101942062378
55 1.786076e+02 2.967125e-02
* time: 6.008654832839966
56 1.786076e+02 1.068274e-02
* time: 6.1342689990997314
57 1.786076e+02 5.355447e-03
* time: 6.226784944534302
58 1.786076e+02 8.180920e-03
* time: 6.311630964279175
59 1.786076e+02 4.935439e-03
* time: 6.396829843521118
60 1.786076e+02 4.387986e-03
* time: 6.51107382774353
61 1.786076e+02 4.388023e-03
* time: 6.671905040740967
62 1.786076e+02 4.387934e-03
* time: 6.8193018436431885
63 1.786076e+02 4.387934e-03
* time: 6.938470840454102FittedPumasModel
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.