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.063201904296875 1 9.433464e+02 6.079483e+02 * time: 2.5907089710235596 2 8.189627e+02 4.423725e+02 * time: 2.6082170009613037 3 5.917683e+02 1.819248e+02 * time: 2.620884895324707 4 5.421783e+02 1.121313e+02 * time: 2.6312670707702637 5 5.255651e+02 7.407230e+01 * time: 2.6404850482940674 6 5.208427e+02 8.699271e+01 * time: 2.6491520404815674 7 5.174883e+02 8.974584e+01 * time: 2.657654047012329 8 5.138523e+02 7.328235e+01 * time: 2.6660568714141846 9 5.109883e+02 4.155805e+01 * time: 2.674194097518921 10 5.094359e+02 3.170517e+01 * time: 2.6820859909057617 11 5.086172e+02 3.327331e+01 * time: 2.689913034439087 12 5.080941e+02 2.942077e+01 * time: 2.6976559162139893 13 5.074009e+02 2.839941e+01 * time: 2.7053050994873047 14 5.059302e+02 3.330093e+01 * time: 2.7125680446624756 15 5.036399e+02 3.172884e+01 * time: 2.7202999591827393 16 5.017004e+02 3.160020e+01 * time: 2.7283830642700195 17 5.008553e+02 2.599524e+01 * time: 2.736406087875366 18 5.005913e+02 2.139314e+01 * time: 2.7440900802612305 19 5.003573e+02 2.134778e+01 * time: 2.7516849040985107 20 4.997249e+02 2.069868e+01 * time: 2.759363889694214 21 4.984453e+02 1.859010e+01 * time: 2.7671759128570557 22 4.959584e+02 2.156209e+01 * time: 2.7750539779663086 23 4.923347e+02 3.030833e+01 * time: 2.7837820053100586 24 4.906916e+02 1.652278e+01 * time: 2.794179916381836 25 4.902955e+02 6.360800e+00 * time: 2.803999900817871 26 4.902870e+02 7.028603e+00 * time: 2.8142449855804443 27 4.902193e+02 1.176895e+00 * time: 2.8241219520568848 28 4.902189e+02 1.170642e+00 * time: 2.8322269916534424 29 4.902186e+02 1.167624e+00 * time: 2.839421033859253 30 4.902145e+02 1.110377e+00 * time: 2.847939968109131 31 4.902079e+02 1.010507e+00 * time: 2.9751529693603516 32 4.901917e+02 9.619218e-01 * time: 2.9828169345855713 33 4.901683e+02 1.001006e+00 * time: 2.9900190830230713 34 4.901473e+02 6.138233e-01 * time: 2.9979419708251953 35 4.901412e+02 1.754342e-01 * time: 3.0057499408721924 36 4.901406e+02 2.617009e-02 * time: 3.0129308700561523 37 4.901405e+02 4.585882e-03 * time: 3.0193898677825928 38 4.901405e+02 7.668184e-04 * time: 3.0249040126800537
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.5020370483398438e-5 1 9.056373e+02 5.541920e+02 * time: 0.9572899341583252 2 7.787963e+02 3.988904e+02 * time: 0.9791219234466553 3 5.915054e+02 1.495777e+02 * time: 0.996955156326294 4 5.533120e+02 8.826583e+01 * time: 1.0139939785003662 5 5.389239e+02 9.144086e+01 * time: 1.0301790237426758 6 5.323499e+02 1.000933e+02 * time: 1.0460541248321533 7 5.270252e+02 8.423844e+01 * time: 1.0612709522247314 8 5.233813e+02 5.194402e+01 * time: 1.0758750438690186 9 5.213366e+02 3.461331e+01 * time: 1.0907680988311768 10 5.200972e+02 3.888113e+01 * time: 1.1057260036468506 11 5.191933e+02 3.556605e+01 * time: 1.1206159591674805 12 5.181335e+02 3.624436e+01 * time: 1.1350131034851074 13 5.161626e+02 4.322775e+01 * time: 1.1493380069732666 14 5.133202e+02 3.722515e+01 * time: 1.1636359691619873 15 5.107758e+02 3.401586e+01 * time: 1.177915096282959 16 5.095157e+02 2.854997e+01 * time: 1.1921839714050293 17 5.090165e+02 2.644560e+01 * time: 1.2068979740142822 18 5.085184e+02 2.744429e+01 * time: 1.2216670513153076 19 5.074309e+02 2.793918e+01 * time: 1.2361819744110107 20 5.053757e+02 2.616169e+01 * time: 1.2519581317901611 21 5.018507e+02 2.257667e+01 * time: 1.2675719261169434 22 4.942495e+02 3.832878e+01 * time: 1.2848000526428223 23 4.940229e+02 5.518159e+01 * time: 1.3062200546264648 24 4.909110e+02 3.042064e+01 * time: 1.327111005783081 25 4.900234e+02 6.929306e+00 * time: 1.3458940982818604 26 4.897974e+02 1.087865e+00 * time: 1.4497900009155273 27 4.897942e+02 6.456402e-01 * time: 1.4645581245422363 28 4.897940e+02 6.467689e-01 * time: 1.4783289432525635 29 4.897939e+02 6.463480e-01 * time: 1.4921369552612305 30 4.897935e+02 6.408914e-01 * time: 1.5061581134796143 31 4.897924e+02 6.208208e-01 * time: 1.5206749439239502 32 4.897900e+02 1.035462e+00 * time: 1.5353069305419922 33 4.897850e+02 1.452099e+00 * time: 1.550157070159912 34 4.897776e+02 1.482593e+00 * time: 1.564682960510254 35 4.897718e+02 8.420646e-01 * time: 1.57973313331604 36 4.897702e+02 2.023876e-01 * time: 1.5948259830474854 37 4.897700e+02 1.885486e-02 * time: 1.6088240146636963 38 4.897700e+02 2.343932e-03 * time: 1.6216390132904053 39 4.897700e+02 4.417566e-04 * time: 1.6331241130828857
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: 1.811981201171875e-5 1 2.005689e+02 1.252637e+02 * time: 1.210677146911621 2 1.875728e+02 4.225518e+01 * time: 1.3963000774383545 3 1.863803e+02 3.429785e+01 * time: 1.5734021663665771 4 1.845581e+02 3.694230e+01 * time: 1.730100154876709 5 1.828416e+02 1.170915e+01 * time: 1.8871140480041504 6 1.823277e+02 1.004017e+01 * time: 2.051359176635742 7 1.810629e+02 5.326781e+00 * time: 2.2129061222076416 8 1.810479e+02 1.767987e+00 * time: 2.3691749572753906 9 1.810272e+02 3.852038e+00 * time: 2.5397541522979736 10 1.809414e+02 1.196625e+01 * time: 2.696727991104126 11 1.806489e+02 1.861441e+01 * time: 2.8544280529022217 12 1.801984e+02 1.361800e+01 * time: 3.420073986053467 13 1.800427e+02 2.176186e+01 * time: 3.5733561515808105 14 1.796555e+02 7.078531e+00 * time: 3.7278101444244385 15 1.795833e+02 1.583428e+01 * time: 3.8817241191864014 16 1.795221e+02 6.549436e+00 * time: 4.038365125656128 17 1.794623e+02 6.609278e+00 * time: 4.192476034164429 18 1.793643e+02 6.616470e+00 * time: 4.349868059158325 19 1.793502e+02 1.025028e+01 * time: 4.510790109634399 20 1.793177e+02 1.200647e+01 * time: 4.6705100536346436 21 1.792426e+02 1.620362e+01 * time: 4.895181179046631 22 1.791563e+02 1.126230e+01 * time: 5.051059007644653 23 1.790991e+02 7.618307e+00 * time: 5.205412149429321 24 1.790757e+02 8.614130e+00 * time: 5.3558290004730225 25 1.790635e+02 4.847798e+00 * time: 5.512222051620483 26 1.790411e+02 5.839338e+00 * time: 5.665977954864502 27 1.789732e+02 1.117508e+01 * time: 5.8198561668396 28 1.789078e+02 1.140932e+01 * time: 5.986653089523315 29 1.788320e+02 7.102422e+00 * time: 6.172871112823486 30 1.787894e+02 5.575299e+00 * time: 6.343005180358887 31 1.787809e+02 7.418901e+00 * time: 6.510292053222656 32 1.787673e+02 8.440524e-01 * time: 6.672630071640015 33 1.787660e+02 5.201083e-01 * time: 6.978147029876709 34 1.787656e+02 3.728389e-01 * time: 7.276082992553711 35 1.787640e+02 8.646511e-01 * time: 7.514747142791748 36 1.787608e+02 2.106472e+00 * time: 7.754310131072998 37 1.787529e+02 3.930618e+00 * time: 7.946609973907471 38 1.787359e+02 5.982352e+00 * time: 8.12300705909729 39 1.787038e+02 7.139001e+00 * time: 8.300923109054565 40 1.786602e+02 5.600927e+00 * time: 8.479238986968994 41 1.786281e+02 1.847314e+00 * time: 8.655007123947144 42 1.786185e+02 4.597875e-01 * time: 8.846814155578613 43 1.786154e+02 1.194584e+00 * time: 9.01254916191101 44 1.786126e+02 1.309331e+00 * time: 9.165341138839722 45 1.786100e+02 8.113507e-01 * time: 9.315551042556763 46 1.786086e+02 1.436346e-01 * time: 9.465598106384277 47 1.786082e+02 1.862690e-01 * time: 9.625348091125488 48 1.786080e+02 2.732386e-01 * time: 9.769129037857056 49 1.786078e+02 1.985965e-01 * time: 9.91295599937439 50 1.786077e+02 5.740849e-02 * time: 10.0592360496521 51 1.786076e+02 3.227551e-02 * time: 10.217489004135132 52 1.786076e+02 5.665847e-02 * time: 10.362248182296753 53 1.786076e+02 4.525676e-02 * time: 10.506667137145996 54 1.786076e+02 1.344593e-02 * time: 10.64712405204773 55 1.786076e+02 5.720009e-03 * time: 10.791321992874146 56 1.786076e+02 1.346492e-02 * time: 10.948012113571167 57 1.786076e+02 9.891198e-03 * time: 11.089081048965454 58 1.786076e+02 3.530323e-03 * time: 11.232403993606567 59 1.786076e+02 1.555915e-03 * time: 11.373042106628418 60 1.786076e+02 2.917438e-03 * time: 11.529258966445923
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.91098
θkaSlow 0.13112
θCL 1.0854
θVc 7.1008
θbioav 0.4802
ωCL 0.088113
ωVc 0.12133
ξbioav 8.1468e-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: 1.2874603271484375e-5 1 2.014577e+02 1.265001e+02 * time: 0.8555738925933838 2 1.880885e+02 4.190708e+01 * time: 1.1580479145050049 3 1.870317e+02 8.825666e+01 * time: 1.8242778778076172 4 1.846027e+02 4.400156e+01 * time: 2.0181238651275635 5 1.834445e+02 1.906624e+01 * time: 2.1711039543151855 6 1.828599e+02 1.113882e+01 * time: 2.3190529346466064 7 1.815719e+02 7.449354e+00 * time: 2.5230519771575928 8 1.815131e+02 2.164677e+00 * time: 2.674468994140625 9 1.814896e+02 2.167318e+00 * time: 2.823741912841797 10 1.814458e+02 4.615739e+00 * time: 2.97320294380188 11 1.813173e+02 9.576970e+00 * time: 3.1231698989868164 12 1.809756e+02 2.052074e+01 * time: 3.2772839069366455 13 1.807625e+02 4.553382e+01 * time: 3.466399908065796 14 1.802224e+02 6.550955e+00 * time: 3.618635892868042 15 1.800862e+02 2.865350e+00 * time: 3.771721839904785 16 1.800780e+02 1.164556e+00 * time: 3.9240498542785645 17 1.800737e+02 7.952740e-01 * time: 4.072427988052368 18 1.800352e+02 4.860851e+00 * time: 4.240604877471924 19 1.800089e+02 5.176186e+00 * time: 4.394202947616577 20 1.799679e+02 4.303034e+00 * time: 4.550156831741333 21 1.799153e+02 4.614276e+00 * time: 4.703222036361694 22 1.798422e+02 1.209297e+01 * time: 4.857235908508301 23 1.796820e+02 1.712382e+01 * time: 5.02333402633667 24 1.794276e+02 1.436482e+01 * time: 5.173588037490845 25 1.793778e+02 4.168256e+00 * time: 5.322816848754883 26 1.793489e+02 1.884194e+00 * time: 5.4954140186309814 27 1.793332e+02 5.545209e+00 * time: 5.658416032791138 28 1.793234e+02 2.897988e+00 * time: 5.806437969207764 29 1.793122e+02 1.485326e+00 * time: 5.953146934509277 30 1.792880e+02 4.103519e+00 * time: 6.101562023162842 31 1.792600e+02 4.635403e+00 * time: 6.26385498046875 32 1.792384e+02 4.259489e+00 * time: 6.413678884506226 33 1.792252e+02 4.420120e+00 * time: 6.563277959823608 34 1.792027e+02 3.265166e+00 * time: 6.713786840438843 35 1.791841e+02 3.379082e+00 * time: 6.877103805541992 36 1.791705e+02 1.661786e+00 * time: 7.028661012649536 37 1.791554e+02 3.454495e+00 * time: 7.173839807510376 38 1.791291e+02 5.535319e+00 * time: 7.320298910140991 39 1.790716e+02 9.819219e+00 * time: 7.466820955276489 40 1.789893e+02 1.114214e+01 * time: 7.6288769245147705 41 1.788602e+02 1.302217e+01 * time: 7.7767698764801025 42 1.787428e+02 5.130620e+00 * time: 7.9244468212127686 43 1.786638e+02 5.464444e+00 * time: 8.072978019714355 44 1.786405e+02 2.240885e+00 * time: 8.232082843780518 45 1.786337e+02 2.660805e+00 * time: 8.380676031112671 46 1.786258e+02 2.067884e+00 * time: 8.534345865249634 47 1.786185e+02 1.627038e+00 * time: 8.685406923294067 48 1.786128e+02 5.834464e-01 * time: 8.867132902145386 49 1.786101e+02 4.091102e-01 * time: 9.019842863082886 50 1.786090e+02 3.054657e-01 * time: 9.170641899108887 51 1.786084e+02 2.265193e-01 * time: 9.321534872055054 52 1.786080e+02 1.005369e-01 * time: 9.472095966339111 53 1.786078e+02 3.557031e-02 * time: 9.63689398765564 54 1.786077e+02 5.001353e-02 * time: 9.783740997314453 55 1.786076e+02 3.412459e-02 * time: 9.931684970855713 56 1.786076e+02 1.192369e-02 * time: 10.081422805786133 57 1.786076e+02 5.846312e-03 * time: 10.243006944656372 58 1.786076e+02 8.968311e-03 * time: 10.388491868972778 59 1.786076e+02 6.350981e-03 * time: 10.532469034194946 60 1.786076e+02 2.725040e-03 * time: 10.678202867507935
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.91098
θkaSlow 0.13112
θCL 1.0854
θVc 7.1008
θbioav 0.4802
ωCL 0.088114
ωVc 0.12133
nbioav 3.8107e7
σ 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.910983 | 0.910984 |
| 2 | θkaSlow | 0.131118 | 0.131117 |
| 3 | θCL | 1.0854 | 1.0854 |
| 4 | θVc | 7.10078 | 7.10076 |
| 5 | θbioav | 0.480199 | 0.480196 |
| 6 | ωCL | 0.0881132 | 0.0881136 |
| 7 | ωVc | 0.121334 | 0.121335 |
| 8 | σ | 0.105449 | 0.105449 |
| 9 | ξbioav | 8.14685e-5 | missing |
| 10 | nbioav | missing | 3.81069e7 |
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.