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 ~ @. ProportionalNormal(μ, σ)
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 ~ @. ProportionalNormal(μ, σ)
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.04231905937194824 1 9.433464e+02 6.079483e+02 * time: 3.1404380798339844 2 8.189627e+02 4.423725e+02 * time: 3.1579980850219727 3 5.917683e+02 1.819248e+02 * time: 3.1714701652526855 4 5.421783e+02 1.121313e+02 * time: 3.182737112045288 5 5.255651e+02 7.407230e+01 * time: 3.192981004714966 6 5.208427e+02 8.699271e+01 * time: 3.2024052143096924 7 5.174883e+02 8.974584e+01 * time: 3.211433172225952 8 5.138523e+02 7.328235e+01 * time: 3.2203471660614014 9 5.109883e+02 4.155805e+01 * time: 3.2291340827941895 10 5.094359e+02 3.170517e+01 * time: 3.2373621463775635 11 5.086172e+02 3.327331e+01 * time: 3.2455360889434814 12 5.080941e+02 2.942077e+01 * time: 3.253608226776123 13 5.074009e+02 2.839941e+01 * time: 3.261691093444824 14 5.059302e+02 3.330093e+01 * time: 3.2695491313934326 15 5.036399e+02 3.172884e+01 * time: 3.277830123901367 16 5.017004e+02 3.160020e+01 * time: 3.286426067352295 17 5.008553e+02 2.599524e+01 * time: 3.295100212097168 18 5.005913e+02 2.139314e+01 * time: 3.3033151626586914 19 5.003573e+02 2.134778e+01 * time: 3.311368227005005 20 4.997249e+02 2.069868e+01 * time: 3.319634199142456 21 4.984453e+02 1.859010e+01 * time: 3.328266143798828 22 4.959584e+02 2.156209e+01 * time: 3.336933135986328 23 4.923347e+02 3.030833e+01 * time: 3.346140146255493 24 4.906916e+02 1.652278e+01 * time: 3.357802152633667 25 4.902955e+02 6.360800e+00 * time: 3.3688771724700928 26 4.902870e+02 7.028603e+00 * time: 3.3803040981292725 27 4.902193e+02 1.176895e+00 * time: 3.39121413230896 28 4.902189e+02 1.170642e+00 * time: 3.4006741046905518 29 4.902186e+02 1.167624e+00 * time: 3.509965181350708 30 4.902145e+02 1.110377e+00 * time: 3.517914056777954 31 4.902079e+02 1.010507e+00 * time: 3.525811195373535 32 4.901917e+02 9.619218e-01 * time: 3.5339860916137695 33 4.901683e+02 1.001006e+00 * time: 3.541361093521118 34 4.901473e+02 6.138233e-01 * time: 3.549687147140503 35 4.901412e+02 1.754342e-01 * time: 3.5580902099609375 36 4.901406e+02 2.617009e-02 * time: 3.5657901763916016 37 4.901405e+02 4.585882e-03 * time: 3.5727341175079346 38 4.901405e+02 7.668184e-04 * time: 3.5788381099700928
FittedPumasModel
Dynamical system type: Closed form
Number of subjects: 32
Observation records: Active Missing
dv: 256 0
Total: 256 0
Number of parameters: Constant Optimized
0 5
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -490.14052
---------------
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: 1.5974044799804688e-5 1 9.056373e+02 5.541920e+02 * time: 1.2158150672912598 2 7.787963e+02 3.988904e+02 * time: 1.2395801544189453 3 5.915054e+02 1.495777e+02 * time: 1.2585251331329346 4 5.533120e+02 8.826583e+01 * time: 1.2763581275939941 5 5.389239e+02 9.144086e+01 * time: 1.2938711643218994 6 5.323499e+02 1.000933e+02 * time: 1.311155080795288 7 5.270252e+02 8.423844e+01 * time: 1.3280529975891113 8 5.233813e+02 5.194402e+01 * time: 1.34433913230896 9 5.213366e+02 3.461331e+01 * time: 1.3603341579437256 10 5.200972e+02 3.888113e+01 * time: 1.3762969970703125 11 5.191933e+02 3.556605e+01 * time: 1.3922669887542725 12 5.181335e+02 3.624436e+01 * time: 1.408890962600708 13 5.161626e+02 4.322775e+01 * time: 1.425595998764038 14 5.133202e+02 3.722515e+01 * time: 1.4413621425628662 15 5.107758e+02 3.401586e+01 * time: 1.457395076751709 16 5.095157e+02 2.854997e+01 * time: 1.4750151634216309 17 5.090165e+02 2.644560e+01 * time: 1.4917449951171875 18 5.085184e+02 2.744429e+01 * time: 1.5080780982971191 19 5.074309e+02 2.793918e+01 * time: 1.524137020111084 20 5.053757e+02 2.616169e+01 * time: 1.540666103363037 21 5.018507e+02 2.257667e+01 * time: 1.5578510761260986 22 4.942495e+02 3.832878e+01 * time: 1.5767691135406494 23 4.940229e+02 5.518159e+01 * time: 1.7241740226745605 24 4.909110e+02 3.042064e+01 * time: 1.7428641319274902 25 4.900234e+02 6.929306e+00 * time: 1.7600150108337402 26 4.897974e+02 1.087865e+00 * time: 1.7771029472351074 27 4.897942e+02 6.456402e-01 * time: 1.793013095855713 28 4.897940e+02 6.467689e-01 * time: 1.807805061340332 29 4.897939e+02 6.463480e-01 * time: 1.8227570056915283 30 4.897935e+02 6.408914e-01 * time: 1.838042974472046 31 4.897924e+02 6.208208e-01 * time: 1.8533251285552979 32 4.897900e+02 1.035462e+00 * time: 1.8690121173858643 33 4.897850e+02 1.452099e+00 * time: 1.8850760459899902 34 4.897776e+02 1.482593e+00 * time: 1.9009480476379395 35 4.897718e+02 8.420646e-01 * time: 1.9169609546661377 36 4.897702e+02 2.023876e-01 * time: 1.9330871105194092 37 4.897700e+02 1.885486e-02 * time: 1.9482049942016602 38 4.897700e+02 2.343932e-03 * time: 1.9623229503631592 39 4.897700e+02 4.417566e-04 * time: 1.9750261306762695
FittedPumasModel
Dynamical system type: Matrix exponential
Number of subjects: 32
Observation records: Active Missing
dv: 256 0
Total: 256 0
Number of parameters: Constant Optimized
0 5
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -489.77002
---------------
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 ~ @. ProportionalNormal(μ, σ)
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 ~ @. ProportionalNormal(μ, σ)
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: 2.193450927734375e-5 1 2.005689e+02 1.252637e+02 * time: 1.0497369766235352 2 1.875728e+02 4.225518e+01 * time: 1.2884299755096436 3 1.863803e+02 3.429785e+01 * time: 2.0901389122009277 4 1.845581e+02 3.694230e+01 * time: 2.2628209590911865 5 1.828416e+02 1.170915e+01 * time: 2.4344639778137207 6 1.823277e+02 1.004017e+01 * time: 2.6077709197998047 7 1.810629e+02 5.326780e+00 * time: 2.7856898307800293 8 1.810479e+02 1.767987e+00 * time: 2.9626529216766357 9 1.810272e+02 3.852039e+00 * time: 3.1361849308013916 10 1.809414e+02 1.196626e+01 * time: 3.3177568912506104 11 1.806489e+02 1.861442e+01 * time: 3.499868869781494 12 1.801984e+02 1.361825e+01 * time: 4.166946887969971 13 1.800426e+02 2.175391e+01 * time: 4.3445470333099365 14 1.796557e+02 7.078229e+00 * time: 4.526432991027832 15 1.795834e+02 1.585242e+01 * time: 4.702404975891113 16 1.795221e+02 6.546886e+00 * time: 4.879318952560425 17 1.794624e+02 6.606140e+00 * time: 5.0510499477386475 18 1.793643e+02 6.628238e+00 * time: 5.224779844284058 19 1.793501e+02 1.024318e+01 * time: 5.399677038192749 20 1.793176e+02 1.199050e+01 * time: 5.659090042114258 21 1.792428e+02 1.615226e+01 * time: 5.8331379890441895 22 1.791565e+02 1.123676e+01 * time: 6.004519939422607 23 1.790992e+02 7.611851e+00 * time: 6.170881986618042 24 1.790758e+02 8.668272e+00 * time: 6.334403038024902 25 1.790636e+02 4.846912e+00 * time: 6.50288987159729 26 1.790414e+02 5.822103e+00 * time: 6.67062783241272 27 1.789730e+02 1.126845e+01 * time: 6.842267036437988 28 1.789084e+02 1.139115e+01 * time: 7.013418912887573 29 1.788319e+02 7.107729e+00 * time: 7.208477020263672 30 1.787893e+02 5.457761e+00 * time: 7.3818018436431885 31 1.787809e+02 7.359149e+00 * time: 7.547767877578735 32 1.787674e+02 9.549387e-01 * time: 7.712594985961914 33 1.787660e+02 5.381940e-01 * time: 7.876867055892944 34 1.787656e+02 3.650524e-01 * time: 8.035398006439209 35 1.787641e+02 8.314986e-01 * time: 8.197913885116577 36 1.787610e+02 2.111502e+00 * time: 8.380579948425293 37 1.787535e+02 3.928538e+00 * time: 8.5432710647583 38 1.787371e+02 6.009016e+00 * time: 8.70745587348938 39 1.787064e+02 7.236735e+00 * time: 8.87142300605774 40 1.786634e+02 5.848360e+00 * time: 9.117348909378052 41 1.786296e+02 2.108942e+00 * time: 9.348917007446289 42 1.786189e+02 3.791390e-01 * time: 9.530720949172974 43 1.786156e+02 1.190402e+00 * time: 9.694607973098755 44 1.786128e+02 1.326856e+00 * time: 9.85732388496399 45 1.786101e+02 8.410299e-01 * time: 10.015598058700562 46 1.786087e+02 1.557514e-01 * time: 10.16713285446167 47 1.786082e+02 1.827705e-01 * time: 10.336817026138306 48 1.786080e+02 2.803568e-01 * time: 10.48990797996521 49 1.786078e+02 2.036475e-01 * time: 10.641746044158936 50 1.786077e+02 6.205898e-02 * time: 10.799782991409302 51 1.786076e+02 3.523014e-02 * time: 10.967859983444214 52 1.786076e+02 6.000919e-02 * time: 11.123373031616211 53 1.786076e+02 4.776445e-02 * time: 11.278180837631226 54 1.786076e+02 1.444654e-02 * time: 11.42906904220581 55 1.786076e+02 6.551888e-03 * time: 11.615489959716797 56 1.786076e+02 1.321199e-02 * time: 11.76119589805603 57 1.786076e+02 9.970889e-03 * time: 11.907744884490967 58 1.786076e+02 3.085645e-03 * time: 12.060862064361572 59 1.786076e+02 1.559154e-03 * time: 12.236327886581421 60 1.786076e+02 2.912983e-03 * time: 12.387266874313354 61 1.786076e+02 2.912983e-03 * time: 12.638726949691772 62 1.786076e+02 2.912980e-03 * time: 12.892318964004517 63 1.786076e+02 2.912980e-03 * time: 13.138473987579346 64 1.786076e+02 2.912980e-03 * time: 13.422587871551514 65 1.786076e+02 3.862035e-03 * time: 16.391402006149292 66 1.786076e+02 2.077673e-03 * time: 16.54551887512207 67 1.786076e+02 1.097656e-03 * time: 16.697502851486206 68 1.786076e+02 1.098484e-03 * time: 16.870916843414307 69 1.786076e+02 1.098484e-03 * time: 17.122769832611084 70 1.786076e+02 1.098484e-03 * time: 17.35396695137024
FittedPumasModel
Dynamical system type: Matrix exponential
Number of subjects: 40
Observation records: Active Missing
dv: 240 0
Total: 240 0
Number of parameters: Constant Optimized
0 9
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: NoObjectiveChange
Log-likelihood value: -178.60759
--------------------
Estimate
--------------------
θkaFast 0.91098
θkaSlow 0.13112
θCL 1.0854
θVc 7.1008
θbioav 0.4802
ωCL 0.088113
ωVc 0.12133
ξbioav 1.8793e-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: 3.1948089599609375e-5 1 2.014577e+02 1.265001e+02 * time: 1.6165671348571777 2 1.880885e+02 4.190708e+01 * time: 2.4033570289611816 3 1.870317e+02 8.825666e+01 * time: 2.573599100112915 4 1.846027e+02 4.400156e+01 * time: 2.797227144241333 5 1.834445e+02 1.906624e+01 * time: 2.970062017440796 6 1.828599e+02 1.113882e+01 * time: 3.1464650630950928 7 1.815719e+02 7.449354e+00 * time: 3.327877998352051 8 1.815131e+02 2.164677e+00 * time: 3.50537109375 9 1.814896e+02 2.167319e+00 * time: 3.7928240299224854 10 1.814458e+02 4.615740e+00 * time: 3.9593570232391357 11 1.813173e+02 9.576971e+00 * time: 4.122282028198242 12 1.809756e+02 2.052077e+01 * time: 4.288645029067993 13 1.807625e+02 4.553362e+01 * time: 4.466207027435303 14 1.802224e+02 6.550877e+00 * time: 4.678085088729858 15 1.800862e+02 2.865558e+00 * time: 4.865941047668457 16 1.800780e+02 1.164627e+00 * time: 5.046549081802368 17 1.800737e+02 7.952373e-01 * time: 5.329648017883301 18 1.800352e+02 4.860542e+00 * time: 5.502113103866577 19 1.800089e+02 5.176831e+00 * time: 5.674833059310913 20 1.799679e+02 4.304132e+00 * time: 5.849648952484131 21 1.799153e+02 4.612400e+00 * time: 6.023215055465698 22 1.798423e+02 1.209424e+01 * time: 6.199089050292969 23 1.796821e+02 1.712227e+01 * time: 6.381019115447998 24 1.794275e+02 1.434611e+01 * time: 6.594538927078247 25 1.793771e+02 4.126316e+00 * time: 6.770749092102051 26 1.793487e+02 1.832912e+00 * time: 6.970639944076538 27 1.793330e+02 5.487581e+00 * time: 7.136390924453735 28 1.793233e+02 2.892631e+00 * time: 7.295491933822632 29 1.793119e+02 1.442037e+00 * time: 7.4974141120910645 30 1.792878e+02 4.112151e+00 * time: 7.685621976852417 31 1.792598e+02 4.604723e+00 * time: 7.894318103790283 32 1.792384e+02 4.252666e+00 * time: 8.05872893333435 33 1.792251e+02 4.413416e+00 * time: 8.230139970779419 34 1.792026e+02 3.215326e+00 * time: 8.403770923614502 35 1.791840e+02 3.434682e+00 * time: 8.60013198852539 36 1.791706e+02 1.613814e+00 * time: 8.780808925628662 37 1.791555e+02 3.468524e+00 * time: 8.953464984893799 38 1.791295e+02 5.601626e+00 * time: 9.126863956451416 39 1.790711e+02 9.972754e+00 * time: 9.316033124923706 40 1.789890e+02 1.140063e+01 * time: 9.482248067855835 41 1.788579e+02 1.271236e+01 * time: 9.645834922790527 42 1.787400e+02 5.348221e+00 * time: 9.814543962478638 43 1.786630e+02 5.312995e+00 * time: 10.008567094802856 44 1.786404e+02 2.703170e+00 * time: 10.171971082687378 45 1.786340e+02 2.740488e+00 * time: 10.334705114364624 46 1.786249e+02 1.915438e+00 * time: 10.49786901473999 47 1.786180e+02 1.467641e+00 * time: 10.662492036819458 48 1.786126e+02 4.657356e-01 * time: 10.838572978973389 49 1.786100e+02 3.512926e-01 * time: 10.99977707862854 50 1.786089e+02 2.606974e-01 * time: 11.162852048873901 51 1.786083e+02 1.858270e-01 * time: 11.328552007675171 52 1.786079e+02 9.828374e-02 * time: 11.512434005737305 53 1.786078e+02 2.453256e-02 * time: 11.682290077209473 54 1.786077e+02 4.221503e-02 * time: 11.848559141159058 55 1.786076e+02 2.776132e-02 * time: 12.021990060806274 56 1.786076e+02 9.874549e-03 * time: 12.209609031677246 57 1.786076e+02 5.349016e-03 * time: 12.385318994522095 58 1.786076e+02 7.577612e-03 * time: 12.551990032196045 59 1.786076e+02 4.424377e-03 * time: 12.718317031860352 60 1.786076e+02 3.969607e-03 * time: 12.916400909423828 61 1.786076e+02 1.861958e-04 * time: 13.082236051559448
FittedPumasModel
Dynamical system type: Matrix exponential
Number of subjects: 40
Observation records: Active Missing
dv: 240 0
Total: 240 0
Number of parameters: Constant Optimized
0 9
Likelihood approximation: FOCE
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -178.60759
-------------------
Estimate
-------------------
θkaFast 0.91098
θkaSlow 0.13112
θCL 1.0854
θVc 7.1008
θbioav 0.4802
ωCL 0.088113
ωVc 0.12133
nbioav 6.4766e7
σ 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.91098 | 0.910979 |
| 2 | θkaSlow | 0.131118 | 0.131118 |
| 3 | θCL | 1.0854 | 1.0854 |
| 4 | θVc | 7.10079 | 7.10081 |
| 5 | θbioav | 0.480201 | 0.480202 |
| 6 | ωCL | 0.0881134 | 0.0881134 |
| 7 | ωVc | 0.121334 | 0.121334 |
| 8 | σ | 0.105449 | 0.105449 |
| 9 | ξbioav | 1.87929e-5 | missing |
| 10 | nbioav | missing | 6.47662e7 |
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.