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.016993999481201172
1 9.433464e+02 6.079483e+02
* time: 0.3506028652191162
2 8.189627e+02 4.423725e+02
* time: 0.35567378997802734
3 5.917683e+02 1.819248e+02
* time: 0.35980987548828125
4 5.421783e+02 1.121313e+02
* time: 0.3629779815673828
5 5.255651e+02 7.407230e+01
* time: 0.3657708168029785
6 5.208427e+02 8.699271e+01
* time: 0.36940884590148926
7 5.174883e+02 8.974584e+01
* time: 0.40048789978027344
8 5.138523e+02 7.328235e+01
* time: 0.4030458927154541
9 5.109883e+02 4.155805e+01
* time: 0.4054219722747803
10 5.094359e+02 3.170517e+01
* time: 0.4075469970703125
11 5.086172e+02 3.327331e+01
* time: 0.4096338748931885
12 5.080941e+02 2.942077e+01
* time: 0.41170597076416016
13 5.074009e+02 2.839941e+01
* time: 0.4138059616088867
14 5.059302e+02 3.330093e+01
* time: 0.41605591773986816
15 5.036399e+02 3.172884e+01
* time: 0.41826796531677246
16 5.017004e+02 3.160020e+01
* time: 0.42084193229675293
17 5.008553e+02 2.599524e+01
* time: 0.42333197593688965
18 5.005913e+02 2.139314e+01
* time: 0.42587900161743164
19 5.003573e+02 2.134778e+01
* time: 0.4282388687133789
20 4.997249e+02 2.069868e+01
* time: 0.4306519031524658
21 4.984453e+02 1.859010e+01
* time: 0.45484089851379395
22 4.959584e+02 2.156209e+01
* time: 0.457291841506958
23 4.923347e+02 3.030833e+01
* time: 0.45973801612854004
24 4.906916e+02 1.652278e+01
* time: 0.4623849391937256
25 4.902955e+02 6.360800e+00
* time: 0.46474480628967285
26 4.902870e+02 7.028603e+00
* time: 0.46723294258117676
27 4.902193e+02 1.176895e+00
* time: 0.46964001655578613
28 4.902189e+02 1.170642e+00
* time: 0.4716348648071289
29 4.902186e+02 1.167624e+00
* time: 0.4734938144683838
30 4.902145e+02 1.110377e+00
* time: 0.4755527973175049
31 4.902079e+02 1.010507e+00
* time: 0.4777078628540039
32 4.901917e+02 9.619218e-01
* time: 0.4797499179840088
33 4.901683e+02 1.001006e+00
* time: 0.48164796829223633
34 4.901473e+02 6.138233e-01
* time: 0.4996058940887451
35 4.901412e+02 1.754342e-01
* time: 0.5018389225006104
36 4.901406e+02 2.617009e-02
* time: 0.503835916519165
37 4.901405e+02 4.585882e-03
* time: 0.5055568218231201
38 4.901405e+02 7.668184e-04
* time: 0.5070269107818604FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Closed form
Log-likelihood value: -490.14052
Number of subjects: 32
Number of parameters: Fixed Optimized
0 5
Observation records: Active Missing
dv: 256 0
Total: 256 0
-----------------
Estimate
-----------------
θCL 0.16025
θVc 10.262
ωCL 0.23505
ωVc 0.10449
σ 0.3582
-----------------fit_gamma = fit(model_gamma, pop, iparams_pk, FOCE())[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter Function value Gradient norm
0 4.219519e+03 5.960578e+03
* time: 5.3882598876953125e-5
1 9.644493e+02 5.541920e+02
* time: 0.041838884353637695
2 8.376083e+02 3.988904e+02
* time: 0.04776191711425781
3 6.503174e+02 1.495777e+02
* time: 0.052551984786987305
4 6.121241e+02 8.826583e+01
* time: 0.057118892669677734
5 5.977360e+02 9.144086e+01
* time: 0.061666011810302734
6 5.911620e+02 1.000933e+02
* time: 0.06609702110290527
7 5.858372e+02 8.423844e+01
* time: 0.07022905349731445
8 5.821934e+02 5.194402e+01
* time: 0.07438087463378906
9 5.801487e+02 3.461331e+01
* time: 0.07834887504577637
10 5.789092e+02 3.888113e+01
* time: 0.10428500175476074
11 5.780054e+02 3.556605e+01
* time: 0.10871601104736328
12 5.769455e+02 3.624436e+01
* time: 0.1129000186920166
13 5.749747e+02 4.322775e+01
* time: 0.11698794364929199
14 5.721322e+02 3.722515e+01
* time: 0.12088894844055176
15 5.695879e+02 3.401586e+01
* time: 0.12487602233886719
16 5.683277e+02 2.854997e+01
* time: 0.12891387939453125
17 5.678285e+02 2.644560e+01
* time: 0.13291096687316895
18 5.673305e+02 2.744429e+01
* time: 0.1367950439453125
19 5.662430e+02 2.793918e+01
* time: 0.14050602912902832
20 5.641877e+02 2.616169e+01
* time: 0.1445300579071045
21 5.606628e+02 2.257667e+01
* time: 0.14914894104003906
22 5.530616e+02 3.832878e+01
* time: 0.1724720001220703
23 5.528349e+02 5.518159e+01
* time: 0.17825603485107422
24 5.497231e+02 3.042064e+01
* time: 0.18340086936950684
25 5.488355e+02 6.929306e+00
* time: 0.18809103965759277
26 5.486095e+02 1.087865e+00
* time: 0.19270801544189453
27 5.486062e+02 6.456402e-01
* time: 0.19695806503295898
28 5.486061e+02 6.467689e-01
* time: 0.20110487937927246
29 5.486060e+02 6.463480e-01
* time: 0.20499086380004883
30 5.486055e+02 6.408914e-01
* time: 0.20890498161315918
31 5.486045e+02 6.208208e-01
* time: 0.21326184272766113
32 5.486020e+02 1.035462e+00
* time: 0.21806597709655762
33 5.485971e+02 1.452099e+00
* time: 0.2319018840789795
34 5.485897e+02 1.482593e+00
* time: 0.23620390892028809
35 5.485839e+02 8.420646e-01
* time: 0.24057984352111816
36 5.485822e+02 2.023876e-01
* time: 0.2447960376739502
37 5.485821e+02 1.885486e-02
* time: 0.24873900413513184
38 5.485821e+02 2.343932e-03
* time: 0.2523689270019531
39 5.485821e+02 4.417566e-04
* time: 0.25571703910827637FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Matrix exponential
Log-likelihood value: -548.58208
Number of subjects: 32
Number of parameters: Fixed Optimized
0 5
Observation records: Active Missing
dv: 256 0
Total: 256 0
-----------------
Estimate
-----------------
θCL 0.16466
θVc 10.329
ωCL 0.23348
ωVc 0.10661
σ 0.35767
-----------------Finally, let’s compare the estimates:
compare_estimates(; lognormal = fit_lognormal, gamma = fit_gamma)5×3 DataFrame
| Row | parameter | lognormal | gamma |
|---|---|---|---|
| String | Float64? | Float64? | |
| 1 | θCL | 0.160253 | 0.164658 |
| 2 | θVc | 10.2617 | 10.3288 |
| 3 | ωCL | 0.235046 | 0.233484 |
| 4 | ωVc | 0.10449 | 0.106611 |
| 5 | σ | 0.358205 | 0.357667 |
As mention above, the mean of a log-normal is \(\exp \left\{ \mu + \frac{\sigma^2}{2} \right\}\).
So let’s compare that with the gamma typical values:
DataFrame(;
parameter = ["θCL", "θVc"],
θLogNormal = [coef(fit_lognormal).θCL, coef(fit_lognormal).θVc],
ELogNormal = [
exp(log(coef(fit_lognormal).θCL) + coef(fit_lognormal).ωCL^2 / 2),
exp(log(coef(fit_lognormal).θVc) + coef(fit_lognormal).ωVc^2 / 2),
],
θGamma = [coef(fit_gamma).θCL, coef(fit_gamma).θVc],
)2×4 DataFrame
| Row | parameter | θLogNormal | ELogNormal | θGamma |
|---|---|---|---|---|
| String | Float64 | Float64 | Float64 | |
| 1 | θCL | 0.160253 | 0.164741 | 0.164658 |
| 2 | θVc | 10.2617 | 10.3178 | 10.3288 |
As you can see the Gaussian model has a slight bias in the estimation of both θCL and θVc.
Let’s also plot the two probability density functions (PDF) for θCL:
plotdataPK = @chain DataFrame(; x = range(0, 0.5; length = 1_000)) begin
@rtransform begin
:LogNormal =
pdf(LogNormal(log(coef(fit_lognormal).θCL), coef(fit_lognormal).ωCL), :x)
:Gamma = pdf(LogNormal(log(coef(fit_gamma).θCL), coef(fit_gamma).ωCL), :x)
end
end
first(plotdataPK, 5)5×3 DataFrame
| Row | x | LogNormal | Gamma |
|---|---|---|---|
| Float64 | Float64 | Float64 | |
| 1 | 0.0 | 0.0 | 0.0 |
| 2 | 0.000500501 | 5.27258e-128 | 5.2502e-131 |
| 3 | 0.001001 | 9.25648e-99 | 3.24235e-101 |
| 4 | 0.0015015 | 2.10132e-83 | 1.45529e-85 |
| 5 | 0.002002 | 2.71651e-73 | 2.9785e-75 |
data(stack(plotdataPK, [:LogNormal, :Gamma])) *
mapping(:x, :value; color = :variable) *
visual(Lines) |> draw
4.2 Bioavaliability Parallel Absorption Simulation
This is a parallel absorption model with bioavaliabity in both the “fast” as the “slow” depots.
First, the traditional approach with a logistic transformation of a Gaussian random variable. This makes the individual relative bioavailibility logit-normally distributed.
model_logitnormal = @model begin
@param begin
θkaFast ∈ RealDomain(; lower = 0)
θkaSlow ∈ RealDomain(; lower = 0)
θCL ∈ RealDomain(; lower = 0)
θVc ∈ RealDomain(; lower = 0)
θbioav ∈ RealDomain(; lower = 0.0, upper = 1.0)
ωCL ∈ RealDomain(; lower = 0)
ωVc ∈ RealDomain(; lower = 0)
# I call this one ξ to distinguish it from ω since the interpretation is NOT a relative error (coefficient of variation)
ξbioav ∈ RealDomain(; lower = 0, init = 0.1)
σ ∈ RealDomain(; lower = 0)
end
@random begin
_CL ~ Gamma(1 / ωCL^2, θCL * ωCL^2)
_Vc ~ Gamma(1 / ωVc^2, θVc * ωVc^2)
# define the latent Gaussian random effect. Notice the logit transform
ηbioavLogit ~ Normal(logit(θbioav), ξbioav)
end
@pre begin
kaFast = θkaFast
kaSlow = θkaSlow
CL = _CL
Vc = _Vc
end
@dosecontrol begin
# _bioav is LogitNormal distributed
_bioav = logistic(ηbioavLogit)
bioav = (; DepotFast = _bioav, DepotSlow = 1 - _bioav)
end
@dynamics begin
DepotFast' = -kaFast * DepotFast
DepotSlow' = -kaSlow * DepotSlow
Central' = kaFast * DepotFast + kaSlow * DepotSlow - CL / Vc * Central
end
@derived begin
μ := @. Central / Vc
dv ~ @. Normal(μ, abs(μ) * σ)
end
endPumasModel
Parameters: θkaFast, θkaSlow, θCL, θVc, θbioav, ωCL, ωVc, ξbioav, σ
Random effects: _CL, _Vc, ηbioavLogit
Covariates:
Dynamical system variables: DepotFast, DepotSlow, Central
Dynamical system type: Matrix exponential
Derived: dv
Observed: dvNow the same model but with the non-Gaussian random-effects using a beta distribution instead of the logit parameterization of the Gaussian distribution:
model_beta = @model begin
@param begin
θkaFast ∈ RealDomain(; lower = 0)
θkaSlow ∈ RealDomain(; lower = 0)
θCL ∈ RealDomain(; lower = 0)
θVc ∈ RealDomain(; lower = 0)
θbioav ∈ RealDomain(; lower = 0.0, upper = 1.0)
ωCL ∈ RealDomain(; lower = 0)
ωVc ∈ RealDomain(; lower = 0)
# We call this one n since the interpretation is like the length of a Binomial distribution
nbioav ∈ RealDomain(; lower = 0, init = 10)
σ ∈ RealDomain(; lower = 0)
end
@random begin
ηCL ~ Gamma(1 / ωCL^2, θCL * ωCL^2)
ηVc ~ Gamma(1 / ωVc^2, θVc * ωVc^2)
# The makes E(_bioav) = θbioav
# See https://en.wikipedia.org/wiki/Beta_distribution
ηbioav ~ Beta(θbioav * nbioav, (1 - θbioav) * nbioav)
end
@pre begin
kaFast = θkaFast
kaSlow = θkaSlow
CL = ηCL
Vc = ηVc
end
@dosecontrol begin
bioav = (; DepotFast = ηbioav, DepotSlow = 1 - ηbioav)
end
@dynamics begin
DepotFast' = -kaFast * DepotFast
DepotSlow' = -kaSlow * DepotSlow
Central' = kaFast * DepotFast + kaSlow * DepotSlow - CL / Vc * Central
end
@derived begin
μ := @. Central / Vc
dv ~ @. Normal(μ, abs(μ) * σ)
end
endPumasModel
Parameters: θkaFast, θkaSlow, θCL, θVc, θbioav, ωCL, ωVc, nbioav, σ
Random effects: ηCL, ηVc, ηbioav
Covariates:
Dynamical system variables: DepotFast, DepotSlow, Central
Dynamical system type: Matrix exponential
Derived: dv
Observed: dvWe have two types of random effects here.
First, as you are already familiar from the previous example, the clearance (CL), volume of concentration (Vc), and absorption rate (ka) have typical values (i.e. fixed effects) and between-subject variability (i.e. random effects) modelled as a gamma distribution.
Second, bioavailability Bioav is modelled as a beta distribution. Generally the beta distribution is parametrized as:
\[\text{Beta}(\alpha, \beta)\]
where both parameters \(\alpha\) and \(\beta\) are shape parameters.
One nice thing about the beta distribution is that it only takes values between and including 0 and 1, i.e. \([0, 1]\). This makes it the perfect candidate to model bioavailability parameters which are generally bounded in that interval. So, we don’t need to do a logistic transformation.
Another nice thing about the beta distribution is that we can use the alternative \((\mu, n)\)-parametrization with with \(\mu\) serving as a mean-value parameter:
\[\text{Beta}(\mu, n)\]
where in the original beta parametrization:
- \(\alpha = \mu n\)
- \(\beta = (1 - \mu) n\)
Hence, our mean is:
\[\operatorname{E}[F] = \mu = \theta_F\]
which, again, does not depend on any other parameters. The variance is
\[\operatorname{Var}(F) = \frac{\mu(1 - \mu)}{n}\]
so similar to the mean of Bernoulli trials.
Now let’s generate some data for the simulation:
dr = DosageRegimen(
DosageRegimen(100; cmt = :DepotFast),
DosageRegimen(100; cmt = :DepotSlow),
)2×10 DataFrame
| Row | time | cmt | amt | evid | ii | addl | rate | duration | ss | route |
|---|---|---|---|---|---|---|---|---|---|---|
| Float64 | Symbol | Float64 | Int8 | Float64 | Int64 | Float64 | Float64 | Int8 | NCA.Route | |
| 1 | 0.0 | DepotFast | 100.0 | 1 | 0.0 | 0 | 0.0 | 0.0 | 0 | NullRoute |
| 2 | 0.0 | DepotSlow | 100.0 | 1 | 0.0 | 0 | 0.0 | 0.0 | 0 | NullRoute |
simtimes = [0.5, 1.0, 2.0, 4.0, 8.0, 24.0]6-element Vector{Float64}:
0.5
1.0
2.0
4.0
8.0
24.0trueparam = (;
θkaFast = 0.9,
θkaSlow = 0.2,
θCL = 1.1,
θVc = 10.0,
θbioav = 0.7,
ωCL = 0.1,
ωVc = 0.1,
nbioav = 40,
σ = 0.1,
)(θkaFast = 0.9,
θkaSlow = 0.2,
θCL = 1.1,
θVc = 10.0,
θbioav = 0.7,
ωCL = 0.1,
ωVc = 0.1,
nbioav = 40,
σ = 0.1,)For simplicity, we just add 20% to the true values for initial values:
initparamBeta = map(t -> 1.2 * t, trueparam)(θkaFast = 1.08,
θkaSlow = 0.24,
θCL = 1.32,
θVc = 12.0,
θbioav = 0.84,
ωCL = 0.12,
ωVc = 0.12,
nbioav = 48.0,
σ = 0.12,)The initial values for the LogitNormal need to have ξbioav defined instead of nbioav:
initparamLogitNormal =
(Base.structdiff(initparamBeta, NamedTuple{(:nbioav,)})..., ξbioav = 0.1)(θkaFast = 1.08,
θkaSlow = 0.24,
θCL = 1.32,
θVc = 12.0,
θbioav = 0.84,
ωCL = 0.12,
ωVc = 0.12,
σ = 0.12,
ξbioav = 0.1,)Setup empty Subjects with the dose information:
skeletonpop = map(i -> Subject(; id = i, events = dr), 1:40)Population
Subjects: 40
Observations: Next, we simulate the data (while setting the seed for reprocibility):
simpop =
Subject.(
simobs(
model_beta,
skeletonpop,
trueparam;
obstimes = simtimes,
rng = Random.seed!(Random.default_rng(), 128),
)
)Finally let’s fit both models:
fit_logitnormal = fit(model_logitnormal, simpop, initparamLogitNormal, FOCE())[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter Function value Gradient norm
0 3.179853e+02 3.723266e+02
* time: 5.698204040527344e-5
1 2.618605e+02 1.558889e+02
* time: 0.10254502296447754
2 2.484481e+02 5.652093e+01
* time: 0.16713285446166992
3 2.466691e+02 3.819960e+01
* time: 0.2308359146118164
4 2.430995e+02 2.281161e+01
* time: 0.2838099002838135
5 2.421492e+02 1.230731e+01
* time: 0.3349289894104004
6 2.412891e+02 1.022578e+01
* time: 0.4566500186920166
7 2.398220e+02 4.018354e+00
* time: 0.5063719749450684
8 2.398069e+02 2.747661e+00
* time: 0.5547029972076416
9 2.398011e+02 1.787833e+00
* time: 0.6033508777618408
10 2.397895e+02 2.923095e+00
* time: 0.6529090404510498
11 2.397511e+02 4.766103e+00
* time: 0.702988862991333
12 2.397274e+02 2.427078e+00
* time: 0.7527618408203125
13 2.396969e+02 2.317178e+00
* time: 0.8027989864349365
14 2.396432e+02 6.390959e+00
* time: 0.8529949188232422
15 2.395437e+02 1.431396e+01
* time: 0.9070849418640137
16 2.394242e+02 1.660862e+01
* time: 0.9612760543823242
17 2.392883e+02 7.918201e+00
* time: 1.060032844543457
18 2.392648e+02 5.359476e+00
* time: 1.116588830947876
19 2.392551e+02 8.504104e-01
* time: 1.1649038791656494
20 2.392536e+02 8.658366e-01
* time: 1.21303391456604
21 2.392502e+02 1.359491e+00
* time: 1.2616398334503174
22 2.392455e+02 1.908257e+00
* time: 1.3114039897918701
23 2.392352e+02 3.434282e+00
* time: 1.365121841430664
24 2.392094e+02 4.019057e+00
* time: 1.4191720485687256
25 2.391904e+02 2.360120e+00
* time: 1.482105016708374
26 2.391639e+02 2.232625e+00
* time: 1.544886827468872
27 2.390985e+02 7.831580e+00
* time: 1.5989069938659668
28 2.390462e+02 1.220982e+01
* time: 1.6910488605499268
29 2.389849e+02 1.372374e+01
* time: 1.7481908798217773
30 2.389295e+02 1.012388e+01
* time: 1.7980568408966064
31 2.389050e+02 7.364830e+00
* time: 1.8541200160980225
32 2.388675e+02 3.493531e+00
* time: 1.9037668704986572
33 2.388348e+02 3.601972e+00
* time: 1.9551060199737549
34 2.388092e+02 9.896791e+00
* time: 2.0091769695281982
35 2.387816e+02 5.816471e+00
* time: 2.0633649826049805
36 2.387667e+02 3.995305e+00
* time: 2.1168789863586426
37 2.387624e+02 2.111619e+00
* time: 2.1708879470825195
38 2.387542e+02 3.340445e+00
* time: 2.2243080139160156
39 2.387468e+02 3.742750e+00
* time: 2.3018908500671387
40 2.387148e+02 8.270335e+00
* time: 2.3513739109039307
41 2.386920e+02 8.042091e+00
* time: 2.399440050125122
42 2.386311e+02 1.786257e+00
* time: 2.448331832885742
43 2.385982e+02 4.295253e+00
* time: 2.4970598220825195
44 2.385769e+02 2.426177e+00
* time: 2.5468239784240723
45 2.385642e+02 4.507279e-01
* time: 2.599139928817749
46 2.385598e+02 1.783922e+00
* time: 2.65162992477417
47 2.385584e+02 4.345023e-01
* time: 2.7043638229370117
48 2.385582e+02 3.504116e-01
* time: 2.7592599391937256
49 2.385581e+02 9.328221e-02
* time: 2.811314821243286
50 2.385580e+02 1.265801e-02
* time: 2.8626649379730225
51 2.385580e+02 2.371756e-03
* time: 2.95894193649292
52 2.385580e+02 2.371756e-03
* time: 3.01590895652771
53 2.385580e+02 2.371756e-03
* time: 3.0803170204162598
54 2.385580e+02 2.371756e-03
* time: 3.1487629413604736
55 2.385580e+02 2.371756e-03
* time: 3.2155539989471436
56 2.385580e+02 2.371756e-03
* time: 3.284559965133667
57 2.385580e+02 2.371756e-03
* time: 3.3542349338531494
58 2.385580e+02 2.371756e-03
* time: 3.380823850631714FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Matrix exponential
Log-likelihood value: -238.55805
Number of subjects: 40
Number of parameters: Fixed Optimized
0 9
Observation records: Active Missing
dv: 240 0
Total: 240 0
---------------------
Estimate
---------------------
θkaFast 1.8907
θkaSlow 0.60771
θCL 1.0797
θVc 11.375
θbioav 0.14394
ωCL 0.08151
ωVc 0.10157
ξbioav 0.12513
σ 0.10441
---------------------fit_beta = fit(model_beta, simpop, initparamBeta, FOCE())[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter Function value Gradient norm
0 3.554658e+02 3.718530e+02
* time: 5.0067901611328125e-5
1 2.991160e+02 1.342892e+02
* time: 0.12029099464416504
2 2.851580e+02 4.511285e+01
* time: 0.2619180679321289
3 2.835992e+02 4.343268e+01
* time: 0.3237738609313965
4 2.821916e+02 4.373720e+01
* time: 0.3754849433898926
5 2.794524e+02 1.253857e+01
* time: 0.4278850555419922
6 2.787960e+02 1.102093e+01
* time: 0.48188304901123047
7 2.770092e+02 1.114240e+01
* time: 0.5395979881286621
8 2.769202e+02 5.089120e+00
* time: 0.5967860221862793
9 2.769088e+02 3.041932e+00
* time: 0.6544899940490723
10 2.768987e+02 1.603805e+00
* time: 0.7104768753051758
11 2.767766e+02 1.411814e+01
* time: 0.7680070400238037
12 2.766170e+02 2.719464e+01
* time: 0.8784968852996826
13 2.764580e+02 3.014173e+01
* time: 0.9316120147705078
14 2.762515e+02 1.669149e+01
* time: 0.9843778610229492
15 2.761434e+02 1.957297e+00
* time: 1.0373950004577637
16 2.761353e+02 3.756442e+00
* time: 1.0897319316864014
17 2.761237e+02 1.341491e+00
* time: 1.1422028541564941
18 2.761152e+02 2.182772e+00
* time: 1.1951918601989746
19 2.761036e+02 3.209881e+00
* time: 1.2488129138946533
20 2.760912e+02 3.070991e+00
* time: 1.302799940109253
21 2.760762e+02 1.990495e+00
* time: 1.3569989204406738
22 2.760587e+02 1.586627e+00
* time: 1.4117820262908936
23 2.760290e+02 2.471592e+00
* time: 1.509998083114624
24 2.759737e+02 1.040809e+01
* time: 1.563704013824463
25 2.758948e+02 9.671126e+00
* time: 1.6165728569030762
26 2.758264e+02 7.207553e+00
* time: 1.6699950695037842
27 2.757728e+02 3.300455e+00
* time: 1.7319090366363525
28 2.757344e+02 4.723183e+00
* time: 1.7877070903778076
29 2.756731e+02 6.515451e+00
* time: 1.8425710201263428
30 2.755936e+02 4.181865e+00
* time: 1.8993630409240723
31 2.755684e+02 5.163221e+00
* time: 1.9684650897979736
32 2.755454e+02 8.132621e+00
* time: 2.054144859313965
33 2.755346e+02 2.188011e+00
* time: 2.107814073562622
34 2.755279e+02 3.127135e+00
* time: 2.1603548526763916
35 2.755221e+02 3.813163e+00
* time: 2.2121899127960205
36 2.754865e+02 6.610197e+00
* time: 2.2647600173950195
37 2.754663e+02 5.265615e+00
* time: 2.317209005355835
38 2.754486e+02 9.998959e-01
* time: 2.3725180625915527
39 2.754421e+02 1.112426e+00
* time: 2.4284889698028564
40 2.754395e+02 3.014869e-01
* time: 2.484915018081665
41 2.754386e+02 5.351895e-01
* time: 2.541581869125366
42 2.754384e+02 1.402404e-01
* time: 2.596606969833374
43 2.754383e+02 1.288533e-01
* time: 2.6790239810943604
44 2.754381e+02 3.845371e-01
* time: 2.730478048324585
45 2.754376e+02 9.172184e-01
* time: 2.782032012939453
46 2.754363e+02 1.766701e+00
* time: 2.833767890930176
47 2.754331e+02 3.007560e+00
* time: 2.885891914367676
48 2.754256e+02 4.573254e+00
* time: 2.938135862350464
49 2.754117e+02 5.874317e+00
* time: 2.993433952331543
50 2.753921e+02 6.633185e+00
* time: 3.0498218536376953
51 2.753744e+02 2.205111e+00
* time: 3.10664701461792
52 2.753702e+02 2.363930e+00
* time: 3.1635239124298096
53 2.753546e+02 8.713730e-01
* time: 3.2202000617980957
54 2.753543e+02 3.795710e+00
* time: 3.302341938018799
55 2.753492e+02 1.784042e+00
* time: 3.353848934173584
56 2.753441e+02 8.514295e-01
* time: 3.404778003692627
57 2.753339e+02 1.341633e+00
* time: 3.4569430351257324
58 2.753266e+02 2.387419e+00
* time: 3.5089609622955322
59 2.753233e+02 1.246148e+00
* time: 3.5605180263519287
60 2.753208e+02 2.973353e-01
* time: 3.615043878555298
61 2.753202e+02 3.696236e-01
* time: 3.669847011566162
62 2.753193e+02 4.553497e-01
* time: 3.7249529361724854
63 2.753186e+02 2.081281e-01
* time: 3.7801358699798584
64 2.753183e+02 8.672527e-02
* time: 3.8351938724517822
65 2.753182e+02 1.015286e-01
* time: 3.9160399436950684
66 2.753181e+02 1.003887e-01
* time: 3.9662959575653076
67 2.753181e+02 5.060689e-02
* time: 4.016520977020264
68 2.753180e+02 2.051541e-02
* time: 4.066173076629639
69 2.753180e+02 1.864290e-02
* time: 4.115004062652588
70 2.753180e+02 1.911113e-02
* time: 4.164100885391235
71 2.753180e+02 8.991085e-03
* time: 4.215550899505615
72 2.753180e+02 4.925215e-03
* time: 4.268024921417236
73 2.753180e+02 2.135748e-03
* time: 4.320307016372681
74 2.753180e+02 2.135926e-03
* time: 4.393136978149414FittedPumasModel
Successful minimization: false
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Matrix exponential
Log-likelihood value: -275.31802
Number of subjects: 40
Number of parameters: Fixed Optimized
0 9
Observation records: Active Missing
dv: 240 0
Total: 240 0
----------------------
Estimate
----------------------
θkaFast 1.9882
θkaSlow 0.61333
θCL 1.0797
θVc 11.381
θbioav 0.13274
ωCL 0.081471
ωVc 0.10178
nbioav 1.4823e7
σ 0.10464
----------------------As before, let’s compare the estimates:
compare_estimates(; logitnormal = fit_logitnormal, beta = fit_beta)10×3 DataFrame
| Row | parameter | logitnormal | beta |
|---|---|---|---|
| String | Float64? | Float64? | |
| 1 | θkaFast | 1.89067 | 1.98819 |
| 2 | θkaSlow | 0.607708 | 0.613333 |
| 3 | θCL | 1.07968 | 1.07967 |
| 4 | θVc | 11.3745 | 11.3811 |
| 5 | θbioav | 0.143942 | 0.132743 |
| 6 | ωCL | 0.0815101 | 0.081471 |
| 7 | ωVc | 0.101565 | 0.101784 |
| 8 | σ | 0.10441 | 0.104644 |
| 9 | ξbioav | 0.125128 | missing |
| 10 | nbioav | missing | 1.48228e7 |
Again, we’ll both PDFs from the estimated values:
plotdatabioav = @chain DataFrame(; x = range(0, 1; length = 1_000)) begin
@rtransform begin
:logitnormal =
1 / coef(fit_logitnormal).ξbioav / √(2π) / (:x * (1 - :x)) * exp(
-(logit(:x) - logit(coef(fit_logitnormal).θbioav))^2 /
(2 * coef(fit_logitnormal).ξbioav^2),
)
:beta = pdf(
Beta(
coef(fit_beta).θbioav * coef(fit_beta).nbioav,
(1 - coef(fit_beta).θbioav) * coef(fit_beta).nbioav,
),
:x,
)
end
end
first(plotdatabioav, 5)5×3 DataFrame
| Row | x | logitnormal | beta |
|---|---|---|---|
| Float64 | Float64 | Float64 | |
| 1 | 0.0 | NaN | 0.0 |
| 2 | 0.001001 | 0.0 | 0.0 |
| 3 | 0.002002 | 1.54002e-269 | 0.0 |
| 4 | 0.003003 | 4.49003e-222 | 0.0 |
| 5 | 0.004004 | 3.80322e-191 | 0.0 |
plt_pdf_bioav =
data(stack(plotdatabioav, [:logitnormal, :beta])) *
mapping(:x, :value; color = :variable) *
visual(Lines);
draw(plt_pdf_bioav; axis = (; xticks = 0.1:0.1:1.0))
For this dataset, the two distributions differ significantly with the Beta model producing a distribution much closer to the truth but for other realizations of the simulated data they are closer to each other.