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.04597115516662598 1 9.433464e+02 6.079483e+02 * time: 1.8780691623687744 2 8.189627e+02 4.423725e+02 * time: 1.899062156677246 3 5.917683e+02 1.819248e+02 * time: 1.912703037261963 4 5.421783e+02 1.121313e+02 * time: 1.9249820709228516 5 5.255651e+02 7.407230e+01 * time: 1.9365110397338867 6 5.208427e+02 8.699271e+01 * time: 1.9479711055755615 7 5.174883e+02 8.974584e+01 * time: 1.9592311382293701 8 5.138523e+02 7.328235e+01 * time: 1.9698700904846191 9 5.109883e+02 4.155805e+01 * time: 1.9796860218048096 10 5.094359e+02 3.170517e+01 * time: 1.989164113998413 11 5.086172e+02 3.327331e+01 * time: 2.000030040740967 12 5.080941e+02 2.942077e+01 * time: 2.0106441974639893 13 5.074009e+02 2.839941e+01 * time: 2.0211970806121826 14 5.059302e+02 3.330093e+01 * time: 2.0314440727233887 15 5.036399e+02 3.172884e+01 * time: 2.042119026184082 16 5.017004e+02 3.160020e+01 * time: 2.05240797996521 17 5.008553e+02 2.599524e+01 * time: 2.0624630451202393 18 5.005913e+02 2.139314e+01 * time: 2.0731611251831055 19 5.003573e+02 2.134778e+01 * time: 2.0841870307922363 20 4.997249e+02 2.069868e+01 * time: 2.094933032989502 21 4.984453e+02 1.859010e+01 * time: 2.10601806640625 22 4.959584e+02 2.156209e+01 * time: 2.117276191711426 23 4.923347e+02 3.030833e+01 * time: 2.1281700134277344 24 4.906916e+02 1.652278e+01 * time: 2.1394519805908203 25 4.902955e+02 6.360800e+00 * time: 2.1513590812683105 26 4.902870e+02 7.028603e+00 * time: 2.1643271446228027 27 4.902193e+02 1.176895e+00 * time: 2.1762921810150146 28 4.902189e+02 1.170642e+00 * time: 2.186432123184204 29 4.902186e+02 1.167624e+00 * time: 2.194944143295288 30 4.902145e+02 1.110377e+00 * time: 2.20440411567688 31 4.902079e+02 1.010507e+00 * time: 2.2139179706573486 32 4.901917e+02 9.619218e-01 * time: 2.2239840030670166 33 4.901683e+02 1.001006e+00 * time: 2.234485149383545 34 4.901473e+02 6.138233e-01 * time: 2.2453300952911377 35 4.901412e+02 1.754342e-01 * time: 2.2561330795288086 36 4.901406e+02 2.617009e-02 * time: 2.266223192214966 37 4.901405e+02 4.585882e-03 * time: 2.2749409675598145 38 4.901405e+02 7.668184e-04 * time: 2.281895160675049
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.6927719116210938e-5 1 9.056373e+02 5.541920e+02 * time: 0.5169129371643066 2 7.787963e+02 3.988904e+02 * time: 0.5380470752716064 3 5.915054e+02 1.495777e+02 * time: 0.5563468933105469 4 5.533120e+02 8.826583e+01 * time: 0.5736510753631592 5 5.389239e+02 9.144086e+01 * time: 0.5903630256652832 6 5.323499e+02 1.000933e+02 * time: 0.743056058883667 7 5.270252e+02 8.423844e+01 * time: 0.7575769424438477 8 5.233813e+02 5.194402e+01 * time: 0.77134108543396 9 5.213366e+02 3.461331e+01 * time: 0.7853739261627197 10 5.200972e+02 3.888113e+01 * time: 0.7994740009307861 11 5.191933e+02 3.556605e+01 * time: 0.8135089874267578 12 5.181335e+02 3.624436e+01 * time: 0.8275079727172852 13 5.161626e+02 4.322775e+01 * time: 0.841804027557373 14 5.133202e+02 3.722515e+01 * time: 0.8554840087890625 15 5.107758e+02 3.401586e+01 * time: 0.8693890571594238 16 5.095157e+02 2.854997e+01 * time: 0.8835248947143555 17 5.090165e+02 2.644560e+01 * time: 0.8970499038696289 18 5.085184e+02 2.744429e+01 * time: 0.9101541042327881 19 5.074309e+02 2.793918e+01 * time: 0.9232690334320068 20 5.053757e+02 2.616169e+01 * time: 0.9365460872650146 21 5.018507e+02 2.257667e+01 * time: 0.9504098892211914 22 4.942495e+02 3.832878e+01 * time: 0.9655919075012207 23 4.940229e+02 5.518159e+01 * time: 0.9840719699859619 24 4.909110e+02 3.042064e+01 * time: 1.0016670227050781 25 4.900234e+02 6.929306e+00 * time: 1.0177199840545654 26 4.897974e+02 1.087865e+00 * time: 1.0336899757385254 27 4.897942e+02 6.456402e-01 * time: 1.0480380058288574 28 4.897940e+02 6.467689e-01 * time: 1.061345100402832 29 4.897939e+02 6.463480e-01 * time: 1.0746419429779053 30 4.897935e+02 6.408914e-01 * time: 1.0882339477539062 31 4.897924e+02 6.208208e-01 * time: 1.1022520065307617 32 4.897900e+02 1.035462e+00 * time: 1.1182188987731934 33 4.897850e+02 1.452099e+00 * time: 1.1346418857574463 34 4.897776e+02 1.482593e+00 * time: 1.1506659984588623 35 4.897718e+02 8.420646e-01 * time: 1.1674160957336426 36 4.897702e+02 2.023876e-01 * time: 1.1841580867767334 37 4.897700e+02 1.885486e-02 * time: 1.1997380256652832 38 4.897700e+02 2.343932e-03 * time: 1.2141940593719482 39 4.897700e+02 4.417566e-04 * time: 1.2279369831085205
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.2901763916015625e-5 1 2.005689e+02 1.252637e+02 * time: 1.0234711170196533 2 1.875728e+02 4.225518e+01 * time: 1.24658203125 3 1.863803e+02 3.429785e+01 * time: 1.4539780616760254 4 1.845581e+02 3.694230e+01 * time: 1.6397490501403809 5 1.828416e+02 1.170915e+01 * time: 1.8272130489349365 6 1.823277e+02 1.004017e+01 * time: 2.0775349140167236 7 1.810629e+02 5.326780e+00 * time: 2.2590599060058594 8 1.810479e+02 1.767987e+00 * time: 2.434329032897949 9 1.810272e+02 3.852037e+00 * time: 2.609889030456543 10 1.809414e+02 1.196625e+01 * time: 2.789309024810791 11 1.806489e+02 1.861441e+01 * time: 2.9697859287261963 12 1.801984e+02 1.361805e+01 * time: 3.9093050956726074 13 1.800427e+02 2.176038e+01 * time: 4.099323034286499 14 1.796556e+02 7.078592e+00 * time: 4.2751219272613525 15 1.795833e+02 1.583768e+01 * time: 4.435286998748779 16 1.795221e+02 6.548930e+00 * time: 4.595155954360962 17 1.794623e+02 6.608680e+00 * time: 6.190507888793945 18 1.793643e+02 6.618623e+00 * time: 6.343122959136963 19 1.793502e+02 1.024896e+01 * time: 6.494726896286011 20 1.793177e+02 1.200354e+01 * time: 6.65095591545105 21 1.792426e+02 1.619402e+01 * time: 6.8053600788116455 22 1.791563e+02 1.125775e+01 * time: 6.958499908447266 23 1.790992e+02 7.616758e+00 * time: 7.1101720333099365 24 1.790757e+02 8.625107e+00 * time: 7.26874303817749 25 1.790635e+02 4.847543e+00 * time: 7.434875965118408 26 1.790412e+02 5.835623e+00 * time: 7.60021710395813 27 1.789731e+02 1.119532e+01 * time: 7.766736030578613 28 1.789079e+02 1.140511e+01 * time: 7.934931039810181 29 1.788320e+02 7.103608e+00 * time: 8.10484004020691 30 1.787894e+02 5.550526e+00 * time: 8.276628971099854 31 1.787809e+02 7.408701e+00 * time: 8.439433097839355 32 1.787673e+02 8.662710e-01 * time: 8.60557508468628 33 1.787660e+02 5.235721e-01 * time: 8.818289995193481 34 1.787656e+02 3.711878e-01 * time: 8.990493059158325 35 1.787640e+02 8.611814e-01 * time: 9.166373014450073 36 1.787608e+02 2.111949e+00 * time: 9.341784954071045 37 1.787530e+02 3.936251e+00 * time: 9.519098043441772 38 1.787361e+02 5.994893e+00 * time: 9.702872037887573 39 1.787043e+02 7.164098e+00 * time: 9.88258409500122 40 1.786607e+02 5.650144e+00 * time: 10.065239906311035 41 1.786284e+02 1.894970e+00 * time: 10.235699892044067 42 1.786185e+02 4.455262e-01 * time: 10.401416063308716 43 1.786154e+02 1.194909e+00 * time: 10.587754011154175 44 1.786126e+02 1.312730e+00 * time: 10.75137996673584 45 1.786100e+02 8.172273e-01 * time: 10.905668020248413 46 1.786086e+02 1.453543e-01 * time: 11.064626932144165 47 1.786082e+02 1.853559e-01 * time: 11.22005295753479 48 1.786080e+02 2.747125e-01 * time: 11.369112014770508 49 1.786078e+02 1.989652e-01 * time: 11.51590895652771 50 1.786077e+02 5.864785e-02 * time: 11.671787977218628 51 1.786076e+02 3.327371e-02 * time: 11.841831922531128 52 1.786076e+02 5.696473e-02 * time: 11.991991996765137 53 1.786076e+02 4.644636e-02 * time: 12.13944411277771 54 1.786076e+02 1.344483e-02 * time: 12.287133932113647 55 1.786076e+02 5.947516e-03 * time: 12.432610988616943 56 1.786076e+02 1.353345e-02 * time: 12.57513689994812 57 1.786076e+02 1.008868e-02 * time: 12.73903203010559 58 1.786076e+02 3.361428e-03 * time: 12.881474018096924 59 1.786076e+02 1.523993e-03 * time: 13.021722078323364 60 1.786076e+02 2.928621e-03 * time: 13.168906927108765 61 1.786076e+02 2.928621e-03 * time: 13.371624946594238 62 1.786076e+02 4.474249e-03 * time: 13.608237981796265 63 1.786076e+02 5.070930e-03 * time: 13.759991884231567 64 1.786076e+02 5.122940e-03 * time: 13.908951044082642 65 1.786076e+02 5.122940e-03 * time: 14.09579610824585 66 1.786076e+02 5.121982e-03 * time: 14.289611101150513 67 1.786076e+02 5.121982e-03 * time: 14.496097087860107 68 1.786076e+02 5.121982e-03 * time: 14.70067811012268 69 1.786076e+02 5.121982e-03 * time: 14.92599606513977 70 1.786076e+02 5.121982e-03 * time: 15.14700698852539
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.9087066650390625e-5 1 2.014577e+02 1.265001e+02 * time: 1.2183561325073242 2 1.880885e+02 4.190708e+01 * time: 1.4384980201721191 3 1.870317e+02 8.825666e+01 * time: 1.6190481185913086 4 1.846027e+02 4.400156e+01 * time: 1.8213460445404053 5 1.834445e+02 1.906624e+01 * time: 1.9997022151947021 6 1.828599e+02 1.113882e+01 * time: 2.17546010017395 7 1.815719e+02 7.449353e+00 * time: 2.423370122909546 8 1.815131e+02 2.164677e+00 * time: 2.611722230911255 9 1.814896e+02 2.167318e+00 * time: 2.7915332317352295 10 1.814458e+02 4.615740e+00 * time: 2.9743480682373047 11 1.813173e+02 9.576966e+00 * time: 3.153503179550171 12 1.809756e+02 2.052080e+01 * time: 3.338219165802002 13 1.807625e+02 4.553350e+01 * time: 3.523097038269043 14 1.802224e+02 6.550827e+00 * time: 3.7580761909484863 15 1.800862e+02 2.865669e+00 * time: 3.938246011734009 16 1.800780e+02 1.164666e+00 * time: 4.114613056182861 17 1.800737e+02 7.952193e-01 * time: 4.290215015411377 18 1.800352e+02 4.860366e+00 * time: 4.47108006477356 19 1.800089e+02 5.177171e+00 * time: 4.667839050292969 20 1.799679e+02 4.304716e+00 * time: 4.8471081256866455 21 1.799153e+02 4.611393e+00 * time: 5.0282440185546875 22 1.798423e+02 1.209482e+01 * time: 5.208768129348755 23 1.796822e+02 1.712120e+01 * time: 5.386247158050537 24 1.794274e+02 1.433687e+01 * time: 5.581980228424072 25 1.793768e+02 4.105433e+00 * time: 5.758422136306763 26 1.793486e+02 1.806678e+00 * time: 5.96169114112854 27 1.793329e+02 5.459074e+00 * time: 6.151771068572998 28 1.793233e+02 2.889752e+00 * time: 6.32658314704895 29 1.793117e+02 1.420025e+00 * time: 6.500113010406494 30 1.792878e+02 4.116831e+00 * time: 6.675004005432129 31 1.792597e+02 4.587919e+00 * time: 6.864106178283691 32 1.792384e+02 4.248794e+00 * time: 7.040878057479858 33 1.792250e+02 4.410458e+00 * time: 7.2147040367126465 34 1.792025e+02 3.186755e+00 * time: 7.394885063171387 35 1.791840e+02 3.455052e+00 * time: 7.590989112854004 36 1.791706e+02 1.601035e+00 * time: 7.768816232681274 37 1.791555e+02 3.485560e+00 * time: 7.944549083709717 38 1.791298e+02 5.624719e+00 * time: 8.122567176818848 39 1.790705e+02 1.007188e+01 * time: 8.318824052810669 40 1.789886e+02 1.153111e+01 * time: 8.499807119369507 41 1.788564e+02 1.256040e+01 * time: 8.668207168579102 42 1.787385e+02 5.450112e+00 * time: 8.841007232666016 43 1.786624e+02 5.174549e+00 * time: 9.021423101425171 44 1.786404e+02 2.969108e+00 * time: 9.192657232284546 45 1.786341e+02 2.776513e+00 * time: 9.364057064056396 46 1.786242e+02 1.847094e+00 * time: 9.541200160980225 47 1.786177e+02 1.331812e+00 * time: 9.721344232559204 48 1.786124e+02 5.074904e-01 * time: 9.891043186187744 49 1.786099e+02 2.928028e-01 * time: 10.06237506866455 50 1.786089e+02 2.292811e-01 * time: 10.237298011779785 51 1.786083e+02 1.602144e-01 * time: 10.423487186431885 52 1.786079e+02 8.960392e-02 * time: 10.592082023620605 53 1.786078e+02 2.250539e-02 * time: 10.75566816329956 54 1.786077e+02 3.750070e-02 * time: 10.925458192825317 55 1.786076e+02 2.547671e-02 * time: 11.080775022506714 56 1.786076e+02 7.939675e-03 * time: 11.25467300415039 57 1.786076e+02 5.044277e-03 * time: 11.411046028137207 58 1.786076e+02 6.897135e-03 * time: 11.560242176055908 59 1.786076e+02 4.595869e-03 * time: 11.708198070526123 60 1.786076e+02 4.342682e-03 * time: 11.88126802444458 61 1.786076e+02 4.342682e-03 * time: 11.96411919593811
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.