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.02742290496826172
1 9.433464e+02 6.079483e+02
* time: 0.5898759365081787
2 8.189627e+02 4.423725e+02
* time: 0.6015350818634033
3 5.917683e+02 1.819248e+02
* time: 0.6104691028594971
4 5.421783e+02 1.121313e+02
* time: 0.6170809268951416
5 5.255651e+02 7.407230e+01
* time: 0.6229109764099121
6 5.208427e+02 8.699271e+01
* time: 0.6688261032104492
7 5.174883e+02 8.974584e+01
* time: 0.6745929718017578
8 5.138523e+02 7.328235e+01
* time: 0.6798639297485352
9 5.109883e+02 4.155805e+01
* time: 0.6848280429840088
10 5.094359e+02 3.170517e+01
* time: 0.6899681091308594
11 5.086172e+02 3.327331e+01
* time: 0.6961710453033447
12 5.080941e+02 2.942077e+01
* time: 0.701308012008667
13 5.074009e+02 2.839941e+01
* time: 0.7065548896789551
14 5.059302e+02 3.330093e+01
* time: 0.7118659019470215
15 5.036399e+02 3.172884e+01
* time: 0.7169070243835449
16 5.017004e+02 3.160020e+01
* time: 0.7225089073181152
17 5.008553e+02 2.599524e+01
* time: 0.7283320426940918
18 5.005913e+02 2.139314e+01
* time: 0.7339780330657959
19 5.003573e+02 2.134778e+01
* time: 0.7409060001373291
20 4.997249e+02 2.069868e+01
* time: 0.7912790775299072
21 4.984453e+02 1.859010e+01
* time: 0.7966780662536621
22 4.959584e+02 2.156209e+01
* time: 0.8017818927764893
23 4.923347e+02 3.030833e+01
* time: 0.8071770668029785
24 4.906916e+02 1.652278e+01
* time: 0.8128509521484375
25 4.902955e+02 6.360800e+00
* time: 0.8182449340820312
26 4.902870e+02 7.028603e+00
* time: 0.8240621089935303
27 4.902193e+02 1.176895e+00
* time: 0.8296310901641846
28 4.902189e+02 1.170642e+00
* time: 0.8344950675964355
29 4.902186e+02 1.167624e+00
* time: 0.8392229080200195
30 4.902145e+02 1.110377e+00
* time: 0.8448901176452637
31 4.902079e+02 1.010507e+00
* time: 0.8508849143981934
32 4.901917e+02 9.619218e-01
* time: 0.8568611145019531
33 4.901683e+02 1.001006e+00
* time: 0.8833889961242676
34 4.901473e+02 6.138233e-01
* time: 0.8886830806732178
35 4.901412e+02 1.754342e-01
* time: 0.8935399055480957
36 4.901406e+02 2.617009e-02
* time: 0.8981199264526367
37 4.901405e+02 4.585882e-03
* time: 0.9023289680480957
38 4.901405e+02 7.668184e-04
* time: 0.9059770107269287FittedPumasModel
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: 7.009506225585938e-5
1 9.644493e+02 5.541920e+02
* time: 0.015491962432861328
2 8.376083e+02 3.988904e+02
* time: 0.027552127838134766
3 6.503174e+02 1.495777e+02
* time: 0.03701210021972656
4 6.121241e+02 8.826583e+01
* time: 0.04558992385864258
5 5.977360e+02 9.144086e+01
* time: 0.05426597595214844
6 5.911620e+02 1.000933e+02
* time: 0.06309795379638672
7 5.858372e+02 8.423844e+01
* time: 0.07215595245361328
8 5.821934e+02 5.194402e+01
* time: 0.08104491233825684
9 5.801487e+02 3.461331e+01
* time: 0.09004402160644531
10 5.789092e+02 3.888113e+01
* time: 0.13114500045776367
11 5.780054e+02 3.556605e+01
* time: 0.13995790481567383
12 5.769455e+02 3.624436e+01
* time: 0.14841699600219727
13 5.749747e+02 4.322775e+01
* time: 0.15669894218444824
14 5.721322e+02 3.722515e+01
* time: 0.16459012031555176
15 5.695879e+02 3.401586e+01
* time: 0.17236590385437012
16 5.683277e+02 2.854997e+01
* time: 0.18023300170898438
17 5.678285e+02 2.644560e+01
* time: 0.18794703483581543
18 5.673305e+02 2.744429e+01
* time: 0.19611501693725586
19 5.662430e+02 2.793918e+01
* time: 0.20341706275939941
20 5.641877e+02 2.616169e+01
* time: 0.21138405799865723
21 5.606628e+02 2.257667e+01
* time: 0.21925711631774902
22 5.530616e+02 3.832878e+01
* time: 0.2511141300201416
23 5.528349e+02 5.518159e+01
* time: 0.26285505294799805
24 5.497231e+02 3.042064e+01
* time: 0.27362704277038574
25 5.488355e+02 6.929306e+00
* time: 0.2837510108947754
26 5.486095e+02 1.087865e+00
* time: 0.29372310638427734
27 5.486062e+02 6.456402e-01
* time: 0.3025541305541992
28 5.486061e+02 6.467689e-01
* time: 0.31062912940979004
29 5.486060e+02 6.463480e-01
* time: 0.318418025970459
30 5.486055e+02 6.408914e-01
* time: 0.3264000415802002
31 5.486045e+02 6.208208e-01
* time: 0.3352820873260498
32 5.486020e+02 1.035462e+00
* time: 0.3566420078277588
33 5.485971e+02 1.452099e+00
* time: 0.3658571243286133
34 5.485897e+02 1.482593e+00
* time: 0.3747739791870117
35 5.485839e+02 8.420646e-01
* time: 0.3838520050048828
36 5.485822e+02 2.023876e-01
* time: 0.392880916595459
37 5.485821e+02 1.885486e-02
* time: 0.40131211280822754
38 5.485821e+02 2.343932e-03
* time: 0.409088134765625
39 5.485821e+02 4.417566e-04
* time: 0.4160768985748291FittedPumasModel
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: 9.179115295410156e-5
1 2.618605e+02 1.558889e+02
* time: 0.3922598361968994
2 2.484481e+02 5.652093e+01
* time: 0.5027158260345459
3 2.466691e+02 3.819960e+01
* time: 0.6111328601837158
4 2.430995e+02 2.281161e+01
* time: 0.6994597911834717
5 2.421492e+02 1.230731e+01
* time: 0.8579978942871094
6 2.412891e+02 1.022578e+01
* time: 0.9463238716125488
7 2.398220e+02 4.018354e+00
* time: 1.0364480018615723
8 2.398069e+02 2.747661e+00
* time: 1.124128818511963
9 2.398011e+02 1.787833e+00
* time: 1.2126619815826416
10 2.397895e+02 2.923095e+00
* time: 1.3037970066070557
11 2.397511e+02 4.766103e+00
* time: 1.4302668571472168
12 2.397274e+02 2.427078e+00
* time: 1.5188188552856445
13 2.396969e+02 2.317178e+00
* time: 1.6075279712677002
14 2.396432e+02 6.390959e+00
* time: 1.6990318298339844
15 2.395437e+02 1.431396e+01
* time: 1.7927398681640625
16 2.394242e+02 1.660862e+01
* time: 1.8866698741912842
17 2.392883e+02 7.918201e+00
* time: 2.0048749446868896
18 2.392648e+02 5.359476e+00
* time: 2.1099510192871094
19 2.392551e+02 8.504104e-01
* time: 2.19770884513855
20 2.392536e+02 8.658366e-01
* time: 2.285472869873047
21 2.392502e+02 1.359491e+00
* time: 2.3779349327087402
22 2.392455e+02 1.908257e+00
* time: 2.4712088108062744
23 2.392352e+02 3.434282e+00
* time: 2.577631950378418
24 2.392094e+02 4.019057e+00
* time: 2.668590784072876
25 2.391904e+02 2.360120e+00
* time: 2.774880886077881
26 2.391639e+02 2.232625e+00
* time: 2.882330894470215
27 2.390985e+02 7.831580e+00
* time: 2.979834794998169
28 2.390462e+02 1.220982e+01
* time: 3.0853848457336426
29 2.389849e+02 1.372374e+01
* time: 3.192112922668457
30 2.389295e+02 1.012388e+01
* time: 3.2825980186462402
31 2.389050e+02 7.364830e+00
* time: 3.388415813446045
32 2.388675e+02 3.493531e+00
* time: 3.4817957878112793
33 2.388348e+02 3.601972e+00
* time: 3.586423873901367
34 2.388092e+02 9.896791e+00
* time: 3.677142858505249
35 2.387816e+02 5.816471e+00
* time: 3.767258882522583
36 2.387667e+02 3.995305e+00
* time: 3.856663942337036
37 2.387624e+02 2.111619e+00
* time: 3.9497108459472656
38 2.387542e+02 3.340445e+00
* time: 4.04262900352478
39 2.387468e+02 3.742750e+00
* time: 4.144883871078491
40 2.387148e+02 8.270335e+00
* time: 4.232879877090454
41 2.386920e+02 8.042091e+00
* time: 4.32021689414978
42 2.386311e+02 1.786257e+00
* time: 4.410562992095947
43 2.385982e+02 4.295253e+00
* time: 4.503330945968628
44 2.385769e+02 2.426177e+00
* time: 4.595682859420776
45 2.385642e+02 4.507279e-01
* time: 4.699665784835815
46 2.385598e+02 1.783922e+00
* time: 4.787667989730835
47 2.385584e+02 4.345023e-01
* time: 4.875746011734009
48 2.385582e+02 3.504116e-01
* time: 4.962616920471191
49 2.385581e+02 9.328221e-02
* time: 5.052999973297119
50 2.385580e+02 1.265801e-02
* time: 5.141560792922974
51 2.385580e+02 2.371756e-03
* time: 5.239545822143555
52 2.385580e+02 2.371756e-03
* time: 5.3445048332214355
53 2.385580e+02 2.371756e-03
* time: 5.463953018188477
54 2.385580e+02 2.371756e-03
* time: 5.593372821807861
55 2.385580e+02 2.371756e-03
* time: 5.725511789321899
56 2.385580e+02 2.371756e-03
* time: 5.843072891235352
57 2.385580e+02 2.371756e-03
* time: 5.96366286277771
58 2.385580e+02 2.371756e-03
* time: 6.0137410163879395FittedPumasModel
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: 7.009506225585938e-5
1 2.991160e+02 1.342892e+02
* time: 0.30132603645324707
2 2.851580e+02 4.511285e+01
* time: 0.42742419242858887
3 2.835992e+02 4.343268e+01
* time: 0.5467081069946289
4 2.821916e+02 4.373720e+01
* time: 0.6433651447296143
5 2.794524e+02 1.253857e+01
* time: 0.7446281909942627
6 2.787960e+02 1.102093e+01
* time: 0.885930061340332
7 2.770092e+02 1.114240e+01
* time: 0.984464168548584
8 2.769202e+02 5.089120e+00
* time: 1.0821120738983154
9 2.769088e+02 3.041932e+00
* time: 1.1805870532989502
10 2.768987e+02 1.603805e+00
* time: 1.2825140953063965
11 2.767766e+02 1.411814e+01
* time: 1.432922124862671
12 2.766170e+02 2.719464e+01
* time: 1.5329110622406006
13 2.764580e+02 3.014173e+01
* time: 1.6329901218414307
14 2.762515e+02 1.669149e+01
* time: 1.7325329780578613
15 2.761434e+02 1.957297e+00
* time: 1.836616039276123
16 2.761353e+02 3.756442e+00
* time: 1.9387249946594238
17 2.761237e+02 1.341491e+00
* time: 2.0556721687316895
18 2.761152e+02 2.182772e+00
* time: 2.155661106109619
19 2.761036e+02 3.209881e+00
* time: 2.2545790672302246
20 2.760912e+02 3.070991e+00
* time: 2.3558311462402344
21 2.760762e+02 1.990495e+00
* time: 2.4625790119171143
22 2.760587e+02 1.586627e+00
* time: 2.5862300395965576
23 2.760290e+02 2.471592e+00
* time: 2.6862261295318604
24 2.759737e+02 1.040809e+01
* time: 2.7870140075683594
25 2.758948e+02 9.671126e+00
* time: 2.8889620304107666
26 2.758264e+02 7.207553e+00
* time: 2.99590802192688
27 2.757728e+02 3.300455e+00
* time: 3.1380741596221924
28 2.757344e+02 4.723183e+00
* time: 3.2379989624023438
29 2.756731e+02 6.515451e+00
* time: 3.3388500213623047
30 2.755936e+02 4.181865e+00
* time: 3.440657138824463
31 2.755684e+02 5.163221e+00
* time: 3.5701091289520264
32 2.755454e+02 8.132621e+00
* time: 3.690001964569092
33 2.755346e+02 2.188011e+00
* time: 3.7900431156158447
34 2.755279e+02 3.127135e+00
* time: 3.8890411853790283
35 2.755221e+02 3.813163e+00
* time: 3.9898581504821777
36 2.754865e+02 6.610197e+00
* time: 4.0919740200042725
37 2.754663e+02 5.265615e+00
* time: 4.2114081382751465
38 2.754486e+02 9.998959e-01
* time: 4.311330080032349
39 2.754421e+02 1.112426e+00
* time: 4.408746004104614
40 2.754395e+02 3.014869e-01
* time: 4.507755994796753
41 2.754386e+02 5.351895e-01
* time: 4.609272003173828
42 2.754384e+02 1.402404e-01
* time: 4.707880973815918
43 2.754383e+02 1.288533e-01
* time: 4.8225531578063965
44 2.754381e+02 3.845371e-01
* time: 4.917587995529175
45 2.754376e+02 9.172184e-01
* time: 5.012852191925049
46 2.754363e+02 1.766701e+00
* time: 5.109822034835815
47 2.754331e+02 3.007560e+00
* time: 5.210883140563965
48 2.754256e+02 4.573254e+00
* time: 5.329270124435425
49 2.754117e+02 5.874317e+00
* time: 5.427206993103027
50 2.753921e+02 6.633185e+00
* time: 5.525132179260254
51 2.753744e+02 2.205111e+00
* time: 5.623741149902344
52 2.753702e+02 2.363930e+00
* time: 5.725638151168823
53 2.753546e+02 8.713730e-01
* time: 5.846112966537476
54 2.753543e+02 3.795710e+00
* time: 5.943678140640259
55 2.753492e+02 1.784042e+00
* time: 6.038923978805542
56 2.753441e+02 8.514295e-01
* time: 6.133052110671997
57 2.753339e+02 1.341633e+00
* time: 6.231644153594971
58 2.753266e+02 2.387419e+00
* time: 6.333148002624512
59 2.753233e+02 1.246148e+00
* time: 6.448354005813599
60 2.753208e+02 2.973353e-01
* time: 6.544607162475586
61 2.753202e+02 3.696236e-01
* time: 6.637044191360474
62 2.753193e+02 4.553497e-01
* time: 6.7320170402526855
63 2.753186e+02 2.081281e-01
* time: 6.832717180252075
64 2.753183e+02 8.672527e-02
* time: 6.930605173110962
65 2.753182e+02 1.015286e-01
* time: 7.068118095397949
66 2.753181e+02 1.003887e-01
* time: 7.1599531173706055
67 2.753181e+02 5.060689e-02
* time: 7.252223014831543
68 2.753180e+02 2.051541e-02
* time: 7.3452980518341064
69 2.753180e+02 1.864290e-02
* time: 7.438685178756714
70 2.753180e+02 1.911113e-02
* time: 7.550669193267822
71 2.753180e+02 8.991085e-03
* time: 7.641968011856079
72 2.753180e+02 4.925215e-03
* time: 7.732015132904053
73 2.753180e+02 2.135748e-03
* time: 7.821482181549072
74 2.753180e+02 2.135926e-03
* time: 7.954208135604858FittedPumasModel
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.