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.040734052658081055 1 9.433464e+02 6.079483e+02 * time: 1.6484060287475586 2 8.189627e+02 4.423725e+02 * time: 1.665153980255127 3 5.917683e+02 1.819248e+02 * time: 1.67787504196167 4 5.421783e+02 1.121313e+02 * time: 1.687757968902588 5 5.255651e+02 7.407230e+01 * time: 1.6969361305236816 6 5.208427e+02 8.699271e+01 * time: 1.7057580947875977 7 5.174883e+02 8.974584e+01 * time: 1.7142770290374756 8 5.138523e+02 7.328235e+01 * time: 1.7228219509124756 9 5.109883e+02 4.155805e+01 * time: 1.7321860790252686 10 5.094359e+02 3.170517e+01 * time: 1.7408039569854736 11 5.086172e+02 3.327331e+01 * time: 1.7493140697479248 12 5.080941e+02 2.942077e+01 * time: 1.9181909561157227 13 5.074009e+02 2.839941e+01 * time: 1.9258949756622314 14 5.059302e+02 3.330093e+01 * time: 1.9331130981445312 15 5.036399e+02 3.172884e+01 * time: 1.9408659934997559 16 5.017004e+02 3.160020e+01 * time: 1.9489710330963135 17 5.008553e+02 2.599524e+01 * time: 1.9570209980010986 18 5.005913e+02 2.139314e+01 * time: 1.9648590087890625 19 5.003573e+02 2.134778e+01 * time: 1.972567081451416 20 4.997249e+02 2.069868e+01 * time: 1.9806771278381348 21 4.984453e+02 1.859010e+01 * time: 1.9888219833374023 22 4.959584e+02 2.156209e+01 * time: 1.9966509342193604 23 4.923347e+02 3.030833e+01 * time: 2.0049431324005127 24 4.906916e+02 1.652278e+01 * time: 2.014130115509033 25 4.902955e+02 6.360800e+00 * time: 2.0224530696868896 26 4.902870e+02 7.028603e+00 * time: 2.0311830043792725 27 4.902193e+02 1.176895e+00 * time: 2.039703130722046 28 4.902189e+02 1.170642e+00 * time: 2.0464539527893066 29 4.902186e+02 1.167624e+00 * time: 2.052335023880005 30 4.902145e+02 1.110377e+00 * time: 2.0594921112060547 31 4.902079e+02 1.010507e+00 * time: 2.066772937774658 32 4.901917e+02 9.619218e-01 * time: 2.0743000507354736 33 4.901683e+02 1.001006e+00 * time: 2.081063985824585 34 4.901473e+02 6.138233e-01 * time: 2.0887341499328613 35 4.901412e+02 1.754342e-01 * time: 2.1005051136016846 36 4.901406e+02 2.617009e-02 * time: 2.107785940170288 37 4.901405e+02 4.585882e-03 * time: 2.1141209602355957 38 4.901405e+02 7.668184e-04 * time: 2.1193110942840576
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: 2.288818359375e-5 1 9.056373e+02 5.541920e+02 * time: 0.6233658790588379 2 7.787963e+02 3.988904e+02 * time: 0.6491479873657227 3 5.915054e+02 1.495777e+02 * time: 0.6698758602142334 4 5.533120e+02 8.826583e+01 * time: 0.6897730827331543 5 5.389239e+02 9.144086e+01 * time: 0.7090640068054199 6 5.323499e+02 1.000933e+02 * time: 0.7278919219970703 7 5.270252e+02 8.423844e+01 * time: 0.7460670471191406 8 5.233813e+02 5.194402e+01 * time: 0.7641470432281494 9 5.213366e+02 3.461331e+01 * time: 0.7817509174346924 10 5.200972e+02 3.888113e+01 * time: 0.8009359836578369 11 5.191933e+02 3.556605e+01 * time: 0.8201179504394531 12 5.181335e+02 3.624436e+01 * time: 0.83974289894104 13 5.161626e+02 4.322775e+01 * time: 0.8590378761291504 14 5.133202e+02 3.722515e+01 * time: 0.8763320446014404 15 5.107758e+02 3.401586e+01 * time: 0.8955700397491455 16 5.095157e+02 2.854997e+01 * time: 0.9150700569152832 17 5.090165e+02 2.644560e+01 * time: 0.933621883392334 18 5.085184e+02 2.744429e+01 * time: 0.9506728649139404 19 5.074309e+02 2.793918e+01 * time: 0.9672858715057373 20 5.053757e+02 2.616169e+01 * time: 0.9857730865478516 21 5.018507e+02 2.257667e+01 * time: 1.0059919357299805 22 4.942495e+02 3.832878e+01 * time: 1.0286738872528076 23 4.940229e+02 5.518159e+01 * time: 1.055772066116333 24 4.909110e+02 3.042064e+01 * time: 1.210886001586914 25 4.900234e+02 6.929306e+00 * time: 1.2322299480438232 26 4.897974e+02 1.087865e+00 * time: 1.2533280849456787 27 4.897942e+02 6.456402e-01 * time: 1.2735838890075684 28 4.897940e+02 6.467689e-01 * time: 1.2920639514923096 29 4.897939e+02 6.463480e-01 * time: 1.3106920719146729 30 4.897935e+02 6.408914e-01 * time: 1.329679012298584 31 4.897924e+02 6.208208e-01 * time: 1.3488428592681885 32 4.897900e+02 1.035462e+00 * time: 1.3678650856018066 33 4.897850e+02 1.452099e+00 * time: 1.386612892150879 34 4.897776e+02 1.482593e+00 * time: 1.4053668975830078 35 4.897718e+02 8.420646e-01 * time: 1.4246890544891357 36 4.897702e+02 2.023876e-01 * time: 1.4436509609222412 37 4.897700e+02 1.885486e-02 * time: 1.461332082748413 38 4.897700e+02 2.343932e-03 * time: 1.4773800373077393 39 4.897700e+02 4.417566e-04 * time: 1.4914588928222656
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 ~ @. Normal(μ, abs(μ) * σ)
end
endPumasModel
Parameters: θkaFast, θkaSlow, θCL, θVc, θbioav, ωCL, ωVc, ξbioav, σ
Random effects: _CL, _Vc, ηbioavLogit
Covariates:
Dynamical system variables: DepotFast, DepotSlow, Central
Dynamical system type: Matrix exponential
Derived: dv
Observed: dvNow the same model but with the non-Gaussian random-effects using a beta distribution instead of the logit parameterization of the Gaussian distribution:
model_beta = @model begin
@param begin
θkaFast ∈ RealDomain(; lower = 0)
θkaSlow ∈ RealDomain(; lower = 0)
θCL ∈ RealDomain(; lower = 0)
θVc ∈ RealDomain(; lower = 0)
θbioav ∈ RealDomain(; lower = 0.0, upper = 1.0)
ωCL ∈ RealDomain(; lower = 0)
ωVc ∈ RealDomain(; lower = 0)
# We call this one n since the interpretation is like the length of a Binomial distribution
nbioav ∈ RealDomain(; lower = 0, init = 10)
σ ∈ RealDomain(; lower = 0)
end
@random begin
ηCL ~ Gamma(1 / ωCL^2, θCL * ωCL^2)
ηVc ~ Gamma(1 / ωVc^2, θVc * ωVc^2)
# The makes E(_bioav) = θbioav
# See https://en.wikipedia.org/wiki/Beta_distribution
ηbioav ~ Beta(θbioav * nbioav, (1 - θbioav) * nbioav)
end
@pre begin
kaFast = θkaFast
kaSlow = θkaSlow
CL = ηCL
Vc = ηVc
end
@dosecontrol begin
bioav = (; DepotFast = ηbioav, DepotSlow = 1 - ηbioav)
end
@dynamics begin
DepotFast' = -kaFast * DepotFast
DepotSlow' = -kaSlow * DepotSlow
Central' = kaFast * DepotFast + kaSlow * DepotSlow - CL / Vc * Central
end
@derived begin
μ := @. Central / Vc
dv ~ @. Normal(μ, abs(μ) * σ)
end
endPumasModel
Parameters: θkaFast, θkaSlow, θCL, θVc, θbioav, ωCL, ωVc, nbioav, σ
Random effects: ηCL, ηVc, ηbioav
Covariates:
Dynamical system variables: DepotFast, DepotSlow, Central
Dynamical system type: Matrix exponential
Derived: dv
Observed: dvWe have two types of random effects here.
First, as you are already familiar from the previous example, the clearance (CL), volume of concentration (Vc), and absorption rate (ka) have typical values (i.e. fixed effects) and between-subject variability (i.e. random effects) modelled as a gamma distribution.
Second, bioavailability Bioav is modelled as a beta distribution. Generally the beta distribution is parametrized as:
\[\text{Beta}(\alpha, \beta)\]
where both parameters \(\alpha\) and \(\beta\) are shape parameters.
One nice thing about the beta distribution is that it only takes values between and including 0 and 1, i.e. \([0, 1]\). This makes it the perfect candidate to model bioavailability parameters which are generally bounded in that interval. So, we don’t need to do a logistic transformation.
Another nice thing about the beta distribution is that we can use the alternative \((\mu, n)\)-parametrization with with \(\mu\) serving as a mean-value parameter:
\[\text{Beta}(\mu, n)\]
where in the original beta parametrization:
- \(\alpha = \mu n\)
- \(\beta = (1 - \mu) n\)
Hence, our mean is:
\[\operatorname{E}[F] = \mu = \theta_F\]
which, again, does not depend on any other parameters. The variance is
\[\operatorname{Var}(F) = \frac{\mu(1 - \mu)}{n}\]
so similar to the mean of Bernoulli trials.
Now let’s generate some data for the simulation:
dr = DosageRegimen(
DosageRegimen(100; cmt = :DepotFast),
DosageRegimen(100; cmt = :DepotSlow),
)2×10 DataFrame
| Row | time | cmt | amt | evid | ii | addl | rate | duration | ss | route |
|---|---|---|---|---|---|---|---|---|---|---|
| Float64 | Symbol | Float64 | Int8 | Float64 | Int64 | Float64 | Float64 | Int8 | NCA.Route | |
| 1 | 0.0 | DepotFast | 100.0 | 1 | 0.0 | 0 | 0.0 | 0.0 | 0 | NullRoute |
| 2 | 0.0 | DepotSlow | 100.0 | 1 | 0.0 | 0 | 0.0 | 0.0 | 0 | NullRoute |
simtimes = [0.5, 1.0, 2.0, 4.0, 8.0, 24.0]6-element Vector{Float64}:
0.5
1.0
2.0
4.0
8.0
24.0trueparam = (;
θkaFast = 0.9,
θkaSlow = 0.2,
θCL = 1.1,
θVc = 10.0,
θbioav = 0.7,
ωCL = 0.1,
ωVc = 0.1,
nbioav = 40,
σ = 0.1,
)(θkaFast = 0.9,
θkaSlow = 0.2,
θCL = 1.1,
θVc = 10.0,
θbioav = 0.7,
ωCL = 0.1,
ωVc = 0.1,
nbioav = 40,
σ = 0.1,)For simplicity, we just add 20% to the true values for initial values:
initparamBeta = map(t -> 1.2 * t, trueparam)(θkaFast = 1.08,
θkaSlow = 0.24,
θCL = 1.32,
θVc = 12.0,
θbioav = 0.84,
ωCL = 0.12,
ωVc = 0.12,
nbioav = 48.0,
σ = 0.12,)The initial values for the LogitNormal need to have ξbioav defined instead of nbioav:
initparamLogitNormal =
(Base.structdiff(initparamBeta, NamedTuple{(:nbioav,)})..., ξbioav = 0.1)(θkaFast = 1.08,
θkaSlow = 0.24,
θCL = 1.32,
θVc = 12.0,
θbioav = 0.84,
ωCL = 0.12,
ωVc = 0.12,
σ = 0.12,
ξbioav = 0.1,)Setup empty Subjects with the dose information:
skeletonpop = map(i -> Subject(; id = i, events = dr), 1:40)Population
Subjects: 40
Observations: Next, we simulate the data (while setting the seed for reprocibility):
Random.seed!(128)
simpop = Subject.(simobs(model_beta, skeletonpop, trueparam; obstimes = simtimes))Finally let’s fit both models:
fit_logitnormal = fit(model_logitnormal, simpop, initparamLogitNormal, FOCE())[ Info: Checking the initial parameter values. [ Info: The initial negative log likelihood and its gradient are finite. Check passed. Iter Function value Gradient norm 0 2.512964e+02 3.650794e+02 * time: 3.0040740966796875e-5 1 2.005689e+02 1.252637e+02 * time: 0.8852999210357666 2 1.875728e+02 4.225518e+01 * time: 1.0686330795288086 3 1.863803e+02 3.429785e+01 * time: 1.239806890487671 4 1.845581e+02 3.694230e+01 * time: 1.385936975479126 5 1.828416e+02 1.170915e+01 * time: 1.5339958667755127 6 1.823277e+02 1.004017e+01 * time: 1.6820728778839111 7 1.810629e+02 5.326780e+00 * time: 1.8251008987426758 8 1.810479e+02 1.767987e+00 * time: 1.9725298881530762 9 1.810272e+02 3.852037e+00 * time: 2.1210150718688965 10 1.809414e+02 1.196625e+01 * time: 2.264986038208008 11 1.806489e+02 1.861441e+01 * time: 2.416796922683716 12 1.801984e+02 1.361805e+01 * time: 3.0460450649261475 13 1.800427e+02 2.176038e+01 * time: 3.1904430389404297 14 1.796556e+02 7.078592e+00 * time: 3.334941864013672 15 1.795833e+02 1.583768e+01 * time: 3.4783689975738525 16 1.795221e+02 6.548930e+00 * time: 3.626533031463623 17 1.794623e+02 6.608680e+00 * time: 3.7702219486236572 18 1.793643e+02 6.618623e+00 * time: 3.914803981781006 19 1.793502e+02 1.024896e+01 * time: 4.060198068618774 20 1.793177e+02 1.200354e+01 * time: 4.208891868591309 21 1.792426e+02 1.619402e+01 * time: 4.353918075561523 22 1.791563e+02 1.125775e+01 * time: 4.502980947494507 23 1.790992e+02 7.616758e+00 * time: 4.732562065124512 24 1.790757e+02 8.625107e+00 * time: 4.873489856719971 25 1.790635e+02 4.847543e+00 * time: 5.017596960067749 26 1.790412e+02 5.835623e+00 * time: 5.162026882171631 27 1.789731e+02 1.119532e+01 * time: 5.304728984832764 28 1.789079e+02 1.140511e+01 * time: 5.454603910446167 29 1.788320e+02 7.103608e+00 * time: 5.607102870941162 30 1.787894e+02 5.550526e+00 * time: 5.753566026687622 31 1.787809e+02 7.408701e+00 * time: 5.89891505241394 32 1.787673e+02 8.662710e-01 * time: 6.071402072906494 33 1.787660e+02 5.235721e-01 * time: 6.211725950241089 34 1.787656e+02 3.711878e-01 * time: 6.346076965332031 35 1.787640e+02 8.611814e-01 * time: 6.482311964035034 36 1.787608e+02 2.111949e+00 * time: 6.61766505241394 37 1.787530e+02 3.936251e+00 * time: 6.753859043121338 38 1.787361e+02 5.994893e+00 * time: 6.89070200920105 39 1.787043e+02 7.164098e+00 * time: 7.0302958488464355 40 1.786607e+02 5.650144e+00 * time: 7.189091920852661 41 1.786284e+02 1.894970e+00 * time: 7.326766014099121 42 1.786185e+02 4.455262e-01 * time: 7.463728904724121 43 1.786154e+02 1.194909e+00 * time: 7.600420951843262 44 1.786126e+02 1.312730e+00 * time: 7.738098859786987 45 1.786100e+02 8.172273e-01 * time: 7.877038955688477 46 1.786086e+02 1.453543e-01 * time: 8.081285953521729 47 1.786082e+02 1.853559e-01 * time: 8.384360074996948 48 1.786080e+02 2.747125e-01 * time: 8.667639970779419 49 1.786078e+02 1.989652e-01 * time: 8.85991907119751 50 1.786077e+02 5.864785e-02 * time: 9.06082797050476 51 1.786076e+02 3.327371e-02 * time: 9.253113031387329 52 1.786076e+02 5.696473e-02 * time: 9.422330856323242 53 1.786076e+02 4.644636e-02 * time: 9.591341972351074 54 1.786076e+02 1.344483e-02 * time: 9.760030031204224 55 1.786076e+02 5.947516e-03 * time: 9.9441499710083 56 1.786076e+02 1.353345e-02 * time: 10.11112093925476 57 1.786076e+02 1.008868e-02 * time: 10.265622854232788 58 1.786076e+02 3.361428e-03 * time: 10.411036968231201 59 1.786076e+02 1.523993e-03 * time: 10.55524206161499 60 1.786076e+02 2.928621e-03 * time: 10.711039066314697 61 1.786076e+02 2.928621e-03 * time: 10.890878915786743 62 1.786076e+02 4.474249e-03 * time: 11.085495948791504 63 1.786076e+02 5.070930e-03 * time: 11.224091053009033 64 1.786076e+02 5.122940e-03 * time: 11.37916088104248 65 1.786076e+02 5.122940e-03 * time: 11.544913053512573 66 1.786076e+02 5.121982e-03 * time: 11.694450855255127 67 1.786076e+02 5.121982e-03 * time: 11.87617301940918 68 1.786076e+02 5.121982e-03 * time: 12.048439979553223 69 1.786076e+02 5.121982e-03 * time: 12.228302001953125 70 1.786076e+02 5.121982e-03 * time: 12.434552907943726
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 3.8363e-6
σ 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: 2.288818359375e-5 1 2.014577e+02 1.265001e+02 * time: 0.805635929107666 2 1.880885e+02 4.190708e+01 * time: 0.9805049896240234 3 1.870317e+02 8.825666e+01 * time: 1.1354548931121826 4 1.846027e+02 4.400156e+01 * time: 1.3182048797607422 5 1.834445e+02 1.906624e+01 * time: 1.6261930465698242 6 1.828599e+02 1.113882e+01 * time: 1.775735855102539 7 1.815719e+02 7.449353e+00 * time: 1.9279699325561523 8 1.815131e+02 2.164677e+00 * time: 2.07989501953125 9 1.814896e+02 2.167318e+00 * time: 2.2291269302368164 10 1.814458e+02 4.615740e+00 * time: 2.380488872528076 11 1.813173e+02 9.576966e+00 * time: 2.530979871749878 12 1.809756e+02 2.052080e+01 * time: 2.6898319721221924 13 1.807625e+02 4.553350e+01 * time: 2.8917760848999023 14 1.802224e+02 6.550827e+00 * time: 3.041895866394043 15 1.800862e+02 2.865669e+00 * time: 3.195106029510498 16 1.800780e+02 1.164666e+00 * time: 3.3446738719940186 17 1.800737e+02 7.952193e-01 * time: 3.4930219650268555 18 1.800352e+02 4.860366e+00 * time: 3.6461288928985596 19 1.800089e+02 5.177171e+00 * time: 3.8138010501861572 20 1.799679e+02 4.304716e+00 * time: 3.968172073364258 21 1.799153e+02 4.611393e+00 * time: 4.12302303314209 22 1.798423e+02 1.209482e+01 * time: 4.278738975524902 23 1.796822e+02 1.712120e+01 * time: 4.430605888366699 24 1.794274e+02 1.433687e+01 * time: 4.597321033477783 25 1.793768e+02 4.105433e+00 * time: 4.747212886810303 26 1.793486e+02 1.806678e+00 * time: 4.919404983520508 27 1.793329e+02 5.459074e+00 * time: 5.0702478885650635 28 1.793233e+02 2.889752e+00 * time: 5.233597040176392 29 1.793117e+02 1.420025e+00 * time: 5.381443023681641 30 1.792878e+02 4.116831e+00 * time: 5.538021087646484 31 1.792597e+02 4.587919e+00 * time: 5.691138029098511 32 1.792384e+02 4.248794e+00 * time: 5.843174934387207 33 1.792250e+02 4.410458e+00 * time: 6.009737014770508 34 1.792025e+02 3.186755e+00 * time: 6.162298917770386 35 1.791840e+02 3.455052e+00 * time: 6.3132500648498535 36 1.791706e+02 1.601035e+00 * time: 6.46602988243103 37 1.791555e+02 3.485560e+00 * time: 6.630998849868774 38 1.791298e+02 5.624719e+00 * time: 6.778714895248413 39 1.790705e+02 1.007188e+01 * time: 6.930199861526489 40 1.789886e+02 1.153111e+01 * time: 7.086639881134033 41 1.788564e+02 1.256040e+01 * time: 7.253225088119507 42 1.787385e+02 5.450112e+00 * time: 7.408982992172241 43 1.786624e+02 5.174549e+00 * time: 7.563453912734985 44 1.786404e+02 2.969108e+00 * time: 7.712846040725708 45 1.786341e+02 2.776513e+00 * time: 7.863362073898315 46 1.786242e+02 1.847094e+00 * time: 8.029717922210693 47 1.786177e+02 1.331812e+00 * time: 8.180491924285889 48 1.786124e+02 5.074904e-01 * time: 8.329438924789429 49 1.786099e+02 2.928028e-01 * time: 8.478920936584473 50 1.786089e+02 2.292811e-01 * time: 8.642493963241577 51 1.786083e+02 1.602144e-01 * time: 8.792470932006836 52 1.786079e+02 8.960392e-02 * time: 8.942414045333862 53 1.786078e+02 2.250539e-02 * time: 9.086040019989014 54 1.786077e+02 3.750070e-02 * time: 9.238970041275024 55 1.786076e+02 2.547671e-02 * time: 9.380466938018799 56 1.786076e+02 7.939675e-03 * time: 9.525609970092773 57 1.786076e+02 5.044277e-03 * time: 9.668581008911133 58 1.786076e+02 6.897135e-03 * time: 9.796664953231812 59 1.786076e+02 4.595869e-03 * time: 9.942738056182861 60 1.786076e+02 4.342682e-03 * time: 10.090517044067383 61 1.786076e+02 4.342682e-03 * time: 10.167239904403687
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: NoXChange
Log-likelihood value: -178.60759
-------------------
Estimate
-------------------
θkaFast 0.91099
θkaSlow 0.13112
θCL 1.0854
θVc 7.1007
θbioav 0.48019
ωCL 0.088114
ωVc 0.12134
nbioav 4.3311e7
σ 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.910979 | 0.910988 |
| 2 | θkaSlow | 0.131118 | 0.131117 |
| 3 | θCL | 1.0854 | 1.0854 |
| 4 | θVc | 7.10079 | 7.10072 |
| 5 | θbioav | 0.4802 | 0.480191 |
| 6 | ωCL | 0.0881133 | 0.0881138 |
| 7 | ωVc | 0.121334 | 0.121335 |
| 8 | σ | 0.105448 | 0.105449 |
| 9 | ξbioav | 3.83627e-6 | missing |
| 10 | nbioav | missing | 4.3311e7 |
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.