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.023321866989135742
1 9.433464e+02 6.079483e+02
* time: 0.9495799541473389
2 8.189627e+02 4.423725e+02
* time: 0.9782040119171143
3 5.917683e+02 1.819248e+02
* time: 1.0000720024108887
4 5.421783e+02 1.121313e+02
* time: 1.0176048278808594
5 5.255651e+02 7.407230e+01
* time: 1.032956838607788
6 5.208427e+02 8.699271e+01
* time: 1.0477538108825684
7 5.174883e+02 8.974584e+01
* time: 1.0622279644012451
8 5.138523e+02 7.328235e+01
* time: 1.0760118961334229
9 5.109883e+02 4.155805e+01
* time: 1.089210033416748
10 5.094359e+02 3.170517e+01
* time: 1.1013619899749756
11 5.086172e+02 3.327331e+01
* time: 1.1141078472137451
12 5.080941e+02 2.942077e+01
* time: 1.1262528896331787
13 5.074009e+02 2.839941e+01
* time: 1.1385838985443115
14 5.059302e+02 3.330093e+01
* time: 1.1502878665924072
15 5.036399e+02 3.172884e+01
* time: 1.163031816482544
16 5.017004e+02 3.160020e+01
* time: 1.1760468482971191
17 5.008553e+02 2.599524e+01
* time: 1.1895718574523926
18 5.005913e+02 2.139314e+01
* time: 1.2018039226531982
19 5.003573e+02 2.134778e+01
* time: 1.2140109539031982
20 4.997249e+02 2.069868e+01
* time: 1.2262978553771973
21 4.984453e+02 1.859010e+01
* time: 1.2388498783111572
22 4.959584e+02 2.156209e+01
* time: 1.2520709037780762
23 4.923347e+02 3.030833e+01
* time: 1.2660958766937256
24 4.906916e+02 1.652278e+01
* time: 1.281385898590088
25 4.902955e+02 6.360800e+00
* time: 1.2959959506988525
26 4.902870e+02 7.028603e+00
* time: 1.3102738857269287
27 4.902193e+02 1.176895e+00
* time: 1.3251848220825195
28 4.902189e+02 1.170642e+00
* time: 1.3358509540557861
29 4.902186e+02 1.167624e+00
* time: 1.3451600074768066
30 4.902145e+02 1.110377e+00
* time: 1.3559179306030273
31 4.902079e+02 1.010507e+00
* time: 1.3673088550567627
32 4.901917e+02 9.619218e-01
* time: 1.437675952911377
33 4.901683e+02 1.001006e+00
* time: 1.4488029479980469
34 4.901473e+02 6.138233e-01
* time: 1.4614408016204834
35 4.901412e+02 1.754342e-01
* time: 1.4742319583892822
36 4.901406e+02 2.617009e-02
* time: 1.4856019020080566
37 4.901405e+02 4.585882e-03
* time: 1.4955790042877197
38 4.901405e+02 7.668184e-04
* time: 1.5039420127868652FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Closed form
Log-likelihood value: -490.14052
Number of subjects: 32
Number of parameters: Fixed Optimized
0 5
Observation records: Active Missing
dv: 256 0
Total: 256 0
-----------------
Estimate
-----------------
θCL 0.16025
θVc 10.262
ωCL 0.23505
ωVc 0.10449
σ 0.3582
-----------------fit_gamma = fit(model_gamma, pop, iparams_pk, FOCE())[ Info: Checking the initial parameter values.
[ Info: The initial negative log likelihood and its gradient are finite. Check passed.
Iter Function value Gradient norm
0 4.160707e+03 5.960578e+03
* time: 5.602836608886719e-5
1 9.056373e+02 5.541920e+02
* time: 0.031707048416137695
2 7.787963e+02 3.988904e+02
* time: 0.05679917335510254
3 5.915054e+02 1.495777e+02
* time: 0.07663607597351074
4 5.533120e+02 8.826583e+01
* time: 0.09508109092712402
5 5.389239e+02 9.144086e+01
* time: 0.11274504661560059
6 5.323499e+02 1.000933e+02
* time: 0.12959504127502441
7 5.270252e+02 8.423844e+01
* time: 0.14529013633728027
8 5.233813e+02 5.194402e+01
* time: 0.2552320957183838
9 5.213366e+02 3.461331e+01
* time: 0.2718980312347412
10 5.200972e+02 3.888113e+01
* time: 0.28837013244628906
11 5.191933e+02 3.556605e+01
* time: 0.3045051097869873
12 5.181335e+02 3.624436e+01
* time: 0.3205680847167969
13 5.161626e+02 4.322775e+01
* time: 0.33653903007507324
14 5.133202e+02 3.722515e+01
* time: 0.35202813148498535
15 5.107758e+02 3.401586e+01
* time: 0.3684370517730713
16 5.095157e+02 2.854997e+01
* time: 0.3846549987792969
17 5.090165e+02 2.644560e+01
* time: 0.3999161720275879
18 5.085184e+02 2.744429e+01
* time: 0.4146890640258789
19 5.074309e+02 2.793918e+01
* time: 0.4295060634613037
20 5.053757e+02 2.616169e+01
* time: 0.44443297386169434
21 5.018507e+02 2.257667e+01
* time: 0.46131420135498047
22 4.942495e+02 3.832878e+01
* time: 0.4790160655975342
23 4.940229e+02 5.518159e+01
* time: 0.502190113067627
24 4.909110e+02 3.042064e+01
* time: 0.5241701602935791
25 4.900234e+02 6.929306e+00
* time: 0.543673038482666
26 4.897974e+02 1.087865e+00
* time: 0.5637311935424805
27 4.897942e+02 6.456402e-01
* time: 0.5803930759429932
28 4.897940e+02 6.467689e-01
* time: 0.595656156539917
29 4.897939e+02 6.463480e-01
* time: 0.6108222007751465
30 4.897935e+02 6.408914e-01
* time: 0.6257810592651367
31 4.897924e+02 6.208208e-01
* time: 0.6421020030975342
32 4.897900e+02 1.035462e+00
* time: 0.658735990524292
33 4.897850e+02 1.452099e+00
* time: 0.6769599914550781
34 4.897776e+02 1.482593e+00
* time: 0.69504714012146
35 4.897718e+02 8.420646e-01
* time: 0.7143480777740479
36 4.897702e+02 2.023876e-01
* time: 0.7343130111694336
37 4.897700e+02 1.885486e-02
* time: 0.7551779747009277
38 4.897700e+02 2.343932e-03
* time: 0.7714481353759766
39 4.897700e+02 4.417566e-04
* time: 0.7853829860687256FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Matrix exponential
Log-likelihood value: -489.77002
Number of subjects: 32
Number of parameters: Fixed Optimized
0 5
Observation records: Active Missing
dv: 256 0
Total: 256 0
-----------------
Estimate
-----------------
θCL 0.16466
θVc 10.329
ωCL 0.23348
ωVc 0.10661
σ 0.35767
-----------------Finally, let’s compare the estimates:
compare_estimates(; lognormal = fit_lognormal, gamma = fit_gamma)5×3 DataFrame
| Row | parameter | lognormal | gamma |
|---|---|---|---|
| String | Float64? | Float64? | |
| 1 | θCL | 0.160253 | 0.164658 |
| 2 | θVc | 10.2617 | 10.3288 |
| 3 | ωCL | 0.235046 | 0.233484 |
| 4 | ωVc | 0.10449 | 0.106611 |
| 5 | σ | 0.358205 | 0.357667 |
As mention above, the mean of a log-normal is \(\exp \left\{ \mu + \frac{\sigma^2}{2} \right\}\).
So let’s compare that with the gamma typical values:
DataFrame(;
parameter = ["θCL", "θVc"],
θLogNormal = [coef(fit_lognormal).θCL, coef(fit_lognormal).θVc],
ELogNormal = [
exp(log(coef(fit_lognormal).θCL) + coef(fit_lognormal).ωCL^2 / 2),
exp(log(coef(fit_lognormal).θVc) + coef(fit_lognormal).ωVc^2 / 2),
],
θGamma = [coef(fit_gamma).θCL, coef(fit_gamma).θVc],
)2×4 DataFrame
| Row | parameter | θLogNormal | ELogNormal | θGamma |
|---|---|---|---|---|
| String | Float64 | Float64 | Float64 | |
| 1 | θCL | 0.160253 | 0.164741 | 0.164658 |
| 2 | θVc | 10.2617 | 10.3178 | 10.3288 |
As you can see the Gaussian model has a slight bias in the estimation of both θCL and θVc.
Let’s also plot the two probability density functions (PDF) for θCL:
plotdataPK = @chain DataFrame(; x = range(0, 0.5; length = 1_000)) begin
@rtransform begin
:LogNormal =
pdf(LogNormal(log(coef(fit_lognormal).θCL), coef(fit_lognormal).ωCL), :x)
:Gamma = pdf(LogNormal(log(coef(fit_gamma).θCL), coef(fit_gamma).ωCL), :x)
end
end
first(plotdataPK, 5)5×3 DataFrame
| Row | x | LogNormal | Gamma |
|---|---|---|---|
| Float64 | Float64 | Float64 | |
| 1 | 0.0 | 0.0 | 0.0 |
| 2 | 0.000500501 | 5.27258e-128 | 5.2502e-131 |
| 3 | 0.001001 | 9.25648e-99 | 3.24235e-101 |
| 4 | 0.0015015 | 2.10132e-83 | 1.45529e-85 |
| 5 | 0.002002 | 2.71651e-73 | 2.9785e-75 |
data(stack(plotdataPK, [:LogNormal, :Gamma])) *
mapping(:x, :value; color = :variable) *
visual(Lines) |> draw
4.2 Bioavaliability Parallel Absorption Simulation
This is a parallel absorption model with bioavaliabity in both the “fast” as the “slow” depots.
First, the traditional approach with a logistic transformation of a Gaussian random variable. This makes the individual relative bioavailibility logit-normally distributed.
model_logitnormal = @model begin
@param begin
θkaFast ∈ RealDomain(; lower = 0)
θkaSlow ∈ RealDomain(; lower = 0)
θCL ∈ RealDomain(; lower = 0)
θVc ∈ RealDomain(; lower = 0)
θbioav ∈ RealDomain(; lower = 0.0, upper = 1.0)
ωCL ∈ RealDomain(; lower = 0)
ωVc ∈ RealDomain(; lower = 0)
# I call this one ξ to distinguish it from ω since the interpretation is NOT a relative error (coefficient of variation)
ξbioav ∈ RealDomain(; lower = 0, init = 0.1)
σ ∈ RealDomain(; lower = 0)
end
@random begin
_CL ~ Gamma(1 / ωCL^2, θCL * ωCL^2)
_Vc ~ Gamma(1 / ωVc^2, θVc * ωVc^2)
# define the latent Gaussian random effect. Notice the logit transform
ηbioavLogit ~ Normal(logit(θbioav), ξbioav)
end
@pre begin
kaFast = θkaFast
kaSlow = θkaSlow
CL = _CL
Vc = _Vc
end
@dosecontrol begin
# _bioav is LogitNormal distributed
_bioav = logistic(ηbioavLogit)
bioav = (; DepotFast = _bioav, DepotSlow = 1 - _bioav)
end
@dynamics begin
DepotFast' = -kaFast * DepotFast
DepotSlow' = -kaSlow * DepotSlow
Central' = kaFast * DepotFast + kaSlow * DepotSlow - CL / Vc * Central
end
@derived begin
μ := @. Central / Vc
dv ~ @. Normal(μ, abs(μ) * σ)
end
endPumasModel
Parameters: θkaFast, θkaSlow, θCL, θVc, θbioav, ωCL, ωVc, ξbioav, σ
Random effects: _CL, _Vc, ηbioavLogit
Covariates:
Dynamical system variables: DepotFast, DepotSlow, Central
Dynamical system type: Matrix exponential
Derived: dv
Observed: dvNow the same model but with the non-Gaussian random-effects using a beta distribution instead of the logit parameterization of the Gaussian distribution:
model_beta = @model begin
@param begin
θkaFast ∈ RealDomain(; lower = 0)
θkaSlow ∈ RealDomain(; lower = 0)
θCL ∈ RealDomain(; lower = 0)
θVc ∈ RealDomain(; lower = 0)
θbioav ∈ RealDomain(; lower = 0.0, upper = 1.0)
ωCL ∈ RealDomain(; lower = 0)
ωVc ∈ RealDomain(; lower = 0)
# We call this one n since the interpretation is like the length of a Binomial distribution
nbioav ∈ RealDomain(; lower = 0, init = 10)
σ ∈ RealDomain(; lower = 0)
end
@random begin
ηCL ~ Gamma(1 / ωCL^2, θCL * ωCL^2)
ηVc ~ Gamma(1 / ωVc^2, θVc * ωVc^2)
# The makes E(_bioav) = θbioav
# See https://en.wikipedia.org/wiki/Beta_distribution
ηbioav ~ Beta(θbioav * nbioav, (1 - θbioav) * nbioav)
end
@pre begin
kaFast = θkaFast
kaSlow = θkaSlow
CL = ηCL
Vc = ηVc
end
@dosecontrol begin
bioav = (; DepotFast = ηbioav, DepotSlow = 1 - ηbioav)
end
@dynamics begin
DepotFast' = -kaFast * DepotFast
DepotSlow' = -kaSlow * DepotSlow
Central' = kaFast * DepotFast + kaSlow * DepotSlow - CL / Vc * Central
end
@derived begin
μ := @. Central / Vc
dv ~ @. Normal(μ, abs(μ) * σ)
end
endPumasModel
Parameters: θkaFast, θkaSlow, θCL, θVc, θbioav, ωCL, ωVc, nbioav, σ
Random effects: ηCL, ηVc, ηbioav
Covariates:
Dynamical system variables: DepotFast, DepotSlow, Central
Dynamical system type: Matrix exponential
Derived: dv
Observed: dvWe have two types of random effects here.
First, as you are already familiar from the previous example, the clearance (CL), volume of concentration (Vc), and absorption rate (ka) have typical values (i.e. fixed effects) and between-subject variability (i.e. random effects) modelled as a gamma distribution.
Second, bioavailability Bioav is modelled as a beta distribution. Generally the beta distribution is parametrized as:
\[\text{Beta}(\alpha, \beta)\]
where both parameters \(\alpha\) and \(\beta\) are shape parameters.
One nice thing about the beta distribution is that it only takes values between and including 0 and 1, i.e. \([0, 1]\). This makes it the perfect candidate to model bioavailability parameters which are generally bounded in that interval. So, we don’t need to do a logistic transformation.
Another nice thing about the beta distribution is that we can use the alternative \((\mu, n)\)-parametrization with with \(\mu\) serving as a mean-value parameter:
\[\text{Beta}(\mu, n)\]
where in the original beta parametrization:
- \(\alpha = \mu n\)
- \(\beta = (1 - \mu) n\)
Hence, our mean is:
\[\operatorname{E}[F] = \mu = \theta_F\]
which, again, does not depend on any other parameters. The variance is
\[\operatorname{Var}(F) = \frac{\mu(1 - \mu)}{n}\]
so similar to the mean of Bernoulli trials.
Now let’s generate some data for the simulation:
dr = DosageRegimen(
DosageRegimen(100; cmt = :DepotFast),
DosageRegimen(100; cmt = :DepotSlow),
)2×10 DataFrame
| Row | time | cmt | amt | evid | ii | addl | rate | duration | ss | route |
|---|---|---|---|---|---|---|---|---|---|---|
| Float64 | Symbol | Float64 | Int8 | Float64 | Int64 | Float64 | Float64 | Int8 | NCA.Route | |
| 1 | 0.0 | DepotFast | 100.0 | 1 | 0.0 | 0 | 0.0 | 0.0 | 0 | NullRoute |
| 2 | 0.0 | DepotSlow | 100.0 | 1 | 0.0 | 0 | 0.0 | 0.0 | 0 | NullRoute |
simtimes = [0.5, 1.0, 2.0, 4.0, 8.0, 24.0]6-element Vector{Float64}:
0.5
1.0
2.0
4.0
8.0
24.0trueparam = (;
θkaFast = 0.9,
θkaSlow = 0.2,
θCL = 1.1,
θVc = 10.0,
θbioav = 0.7,
ωCL = 0.1,
ωVc = 0.1,
nbioav = 40,
σ = 0.1,
)(θkaFast = 0.9,
θkaSlow = 0.2,
θCL = 1.1,
θVc = 10.0,
θbioav = 0.7,
ωCL = 0.1,
ωVc = 0.1,
nbioav = 40,
σ = 0.1,)For simplicity, we just add 20% to the true values for initial values:
initparamBeta = map(t -> 1.2 * t, trueparam)(θkaFast = 1.08,
θkaSlow = 0.24,
θCL = 1.32,
θVc = 12.0,
θbioav = 0.84,
ωCL = 0.12,
ωVc = 0.12,
nbioav = 48.0,
σ = 0.12,)The initial values for the LogitNormal need to have ξbioav defined instead of nbioav:
initparamLogitNormal =
(Base.structdiff(initparamBeta, NamedTuple{(:nbioav,)})..., ξbioav = 0.1)(θkaFast = 1.08,
θkaSlow = 0.24,
θCL = 1.32,
θVc = 12.0,
θbioav = 0.84,
ωCL = 0.12,
ωVc = 0.12,
σ = 0.12,
ξbioav = 0.1,)Setup empty Subjects with the dose information:
skeletonpop = map(i -> Subject(; id = i, events = dr), 1:40)Population
Subjects: 40
Observations: Next, we simulate the data (while setting the seed for reprocibility):
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: 8.606910705566406e-5
1 2.005689e+02 1.252637e+02
* time: 0.3123788833618164
2 1.875728e+02 4.225518e+01
* time: 0.4525790214538574
3 1.863803e+02 3.429785e+01
* time: 0.5735318660736084
4 1.845581e+02 3.694229e+01
* time: 0.6681790351867676
5 1.828416e+02 1.170915e+01
* time: 0.7723720073699951
6 1.823277e+02 1.004017e+01
* time: 0.8736698627471924
7 1.810629e+02 5.326780e+00
* time: 0.9710860252380371
8 1.810479e+02 1.767987e+00
* time: 1.0630109310150146
9 1.810272e+02 3.852037e+00
* time: 1.1556029319763184
10 1.809414e+02 1.196625e+01
* time: 1.252328872680664
11 1.806489e+02 1.861440e+01
* time: 1.4375710487365723
12 1.801984e+02 1.361774e+01
* time: 1.7836909294128418
13 1.800427e+02 2.177039e+01
* time: 1.8821258544921875
14 1.796554e+02 7.079039e+00
* time: 1.981132984161377
15 1.795832e+02 1.581499e+01
* time: 2.083469867706299
16 1.795220e+02 6.552120e+00
* time: 2.1899218559265137
17 1.794621e+02 6.612645e+00
* time: 2.2994489669799805
18 1.793643e+02 6.603903e+00
* time: 2.410783052444458
19 1.793502e+02 1.025787e+01
* time: 2.5205888748168945
20 1.793178e+02 1.202199e+01
* time: 2.6309380531311035
21 1.792424e+02 1.625661e+01
* time: 2.790705919265747
22 1.791560e+02 1.128856e+01
* time: 2.89349102973938
23 1.790991e+02 7.625637e+00
* time: 2.9939050674438477
24 1.790756e+02 8.557258e+00
* time: 3.092802047729492
25 1.790633e+02 4.848482e+00
* time: 3.2034640312194824
26 1.790408e+02 5.859306e+00
* time: 3.3146040439605713
27 1.789734e+02 1.106035e+01
* time: 3.4254848957061768
28 1.789070e+02 1.143361e+01
* time: 3.5353519916534424
29 1.788321e+02 7.098448e+00
* time: 3.652029037475586
30 1.787896e+02 5.709857e+00
* time: 3.767518997192383
31 1.787809e+02 7.456555e+00
* time: 3.8812530040740967
32 1.787671e+02 7.320931e-01
* time: 3.9962120056152344
33 1.787660e+02 5.023990e-01
* time: 4.109953880310059
34 1.787656e+02 3.813994e-01
* time: 4.278822898864746
35 1.787639e+02 8.608909e-01
* time: 4.372612953186035
36 1.787607e+02 2.047321e+00
* time: 4.464963912963867
37 1.787525e+02 3.859529e+00
* time: 4.563003063201904
38 1.787349e+02 5.864920e+00
* time: 4.669327974319458
39 1.787017e+02 6.966045e+00
* time: 4.780895948410034
40 1.786574e+02 5.348939e+00
* time: 4.896574020385742
41 1.786270e+02 1.642025e+00
* time: 5.008840084075928
42 1.786181e+02 5.211823e-01
* time: 5.117498874664307
43 1.786153e+02 1.186061e+00
* time: 5.223114967346191
44 1.786125e+02 1.292005e+00
* time: 5.327185869216919
45 1.786099e+02 7.814598e-01
* time: 5.432420015335083
46 1.786086e+02 1.369456e-01
* time: 5.537523031234741
47 1.786082e+02 1.912170e-01
* time: 5.6850409507751465
48 1.786080e+02 2.670802e-01
* time: 5.77911901473999
49 1.786078e+02 1.979262e-01
* time: 5.87620997428894
50 1.786077e+02 5.177918e-02
* time: 5.972462892532349
51 1.786076e+02 2.998328e-02
* time: 6.073517084121704
52 1.786076e+02 5.799706e-02
* time: 6.172092914581299
53 1.786076e+02 4.280434e-02
* time: 6.269625902175903
54 1.786076e+02 1.467829e-02
* time: 6.3693718910217285
55 1.786076e+02 6.928178e-03
* time: 6.46371603012085
56 1.786076e+02 1.236015e-02
* time: 6.5639870166778564
57 1.786076e+02 9.950826e-03
* time: 6.663352966308594
58 1.786076e+02 2.974370e-03
* time: 6.763549089431763
59 1.786076e+02 1.374576e-03
* time: 6.863864898681641
60 1.786076e+02 2.860078e-03
* time: 6.95699405670166
61 1.786076e+02 2.100127e-03
* time: 7.049622058868408
62 1.786076e+02 2.100127e-03
* time: 7.2214789390563965
63 1.786076e+02 6.412834e-03
* time: 7.32464599609375
64 1.786076e+02 6.412834e-03
* time: 7.47111701965332
65 1.786076e+02 6.412834e-03
* time: 7.67839789390564
66 1.786076e+02 6.412834e-03
* time: 7.979866981506348
67 1.786076e+02 6.412834e-03
* time: 8.190343856811523FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Matrix exponential
Log-likelihood value: -178.60759
Number of subjects: 40
Number of parameters: Fixed Optimized
0 9
Observation records: Active Missing
dv: 240 0
Total: 240 0
-----------------------
Estimate
-----------------------
θkaFast 0.91097
θkaSlow 0.13112
θCL 1.0854
θVc 7.1008
θbioav 0.4802
ωCL 0.088113
ωVc 0.12133
ξbioav 1.8429e-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: 5.888938903808594e-5
1 2.014577e+02 1.265001e+02
* time: 0.3328239917755127
2 1.880885e+02 4.190708e+01
* time: 0.4948408603668213
3 1.870317e+02 8.825666e+01
* time: 0.6129148006439209
4 1.846027e+02 4.400156e+01
* time: 0.7584187984466553
5 1.834445e+02 1.906624e+01
* time: 0.8778209686279297
6 1.828599e+02 1.113882e+01
* time: 0.9931468963623047
7 1.815719e+02 7.449355e+00
* time: 1.1115097999572754
8 1.815131e+02 2.164678e+00
* time: 1.2275359630584717
9 1.814896e+02 2.167319e+00
* time: 1.4205148220062256
10 1.814458e+02 4.615738e+00
* time: 1.5260868072509766
11 1.813173e+02 9.576967e+00
* time: 1.6331589221954346
12 1.809756e+02 2.052077e+01
* time: 1.7454698085784912
13 1.807625e+02 4.553366e+01
* time: 1.866729974746704
14 1.802224e+02 6.550892e+00
* time: 1.9833929538726807
15 1.800862e+02 2.865509e+00
* time: 2.102334976196289
16 1.800780e+02 1.164611e+00
* time: 2.2184908390045166
17 1.800737e+02 7.952462e-01
* time: 2.3358278274536133
18 1.800352e+02 4.860618e+00
* time: 2.4553959369659424
19 1.800089e+02 5.176689e+00
* time: 2.5746448040008545
20 1.799679e+02 4.303892e+00
* time: 2.695615768432617
21 1.799153e+02 4.612832e+00
* time: 2.8187828063964844
22 1.798423e+02 1.209387e+01
* time: 3.005941867828369
23 1.796821e+02 1.712256e+01
* time: 3.115203857421875
24 1.794275e+02 1.435100e+01
* time: 3.2260568141937256
25 1.793773e+02 4.137313e+00
* time: 3.33463191986084
26 1.793488e+02 1.846545e+00
* time: 3.4825668334960938
27 1.793331e+02 5.502533e+00
* time: 3.5994417667388916
28 1.793234e+02 2.894037e+00
* time: 3.7149438858032227
29 1.793119e+02 1.453372e+00
* time: 3.8290138244628906
30 1.792879e+02 4.109884e+00
* time: 3.943673849105835
31 1.792599e+02 4.613173e+00
* time: 4.058923959732056
32 1.792384e+02 4.254549e+00
* time: 4.175391912460327
33 1.792251e+02 4.415024e+00
* time: 4.290223836898804
34 1.792026e+02 3.229145e+00
* time: 4.406689882278442
35 1.791841e+02 3.422094e+00
* time: 4.573133945465088
36 1.791705e+02 1.621228e+00
* time: 4.680441856384277
37 1.791555e+02 3.462667e+00
* time: 4.7831408977508545
38 1.791293e+02 5.586685e+00
* time: 4.885271787643433
39 1.790713e+02 9.927638e+00
* time: 4.997097969055176
40 1.789891e+02 1.133271e+01
* time: 5.110766887664795
41 1.788585e+02 1.279172e+01
* time: 5.223315954208374
42 1.787407e+02 5.291681e+00
* time: 5.339637994766235
43 1.786633e+02 5.367971e+00
* time: 5.456313848495483
44 1.786405e+02 2.571548e+00
* time: 5.5684168338775635
45 1.786339e+02 2.720314e+00
* time: 5.682646989822388
46 1.786252e+02 1.930583e+00
* time: 5.7967259883880615
47 1.786182e+02 1.523164e+00
* time: 5.906875848770142
48 1.786126e+02 5.027920e-01
* time: 6.071589946746826
49 1.786100e+02 3.733023e-01
* time: 6.1722328662872314
50 1.786089e+02 2.749553e-01
* time: 6.273086786270142
51 1.786083e+02 1.981772e-01
* time: 6.372139930725098
52 1.786079e+02 1.005262e-01
* time: 6.472882986068726
53 1.786078e+02 2.549641e-02
* time: 6.5798749923706055
54 1.786077e+02 4.452931e-02
* time: 6.684991836547852
55 1.786076e+02 2.967125e-02
* time: 6.791165828704834
56 1.786076e+02 1.068274e-02
* time: 6.897058963775635
57 1.786076e+02 5.355447e-03
* time: 7.001224994659424
58 1.786076e+02 8.180920e-03
* time: 7.103579998016357
59 1.786076e+02 4.935439e-03
* time: 7.204744815826416
60 1.786076e+02 4.387986e-03
* time: 7.33115291595459
61 1.786076e+02 4.388023e-03
* time: 7.49640679359436
62 1.786076e+02 4.387934e-03
* time: 7.6858978271484375
63 1.786076e+02 4.387934e-03
* time: 7.793101787567139FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Matrix exponential
Log-likelihood value: -178.60759
Number of subjects: 40
Number of parameters: Fixed Optimized
0 9
Observation records: Active Missing
dv: 240 0
Total: 240 0
----------------------
Estimate
----------------------
θkaFast 0.91099
θkaSlow 0.13112
θCL 1.0854
θVc 7.1007
θbioav 0.48019
ωCL 0.088114
ωVc 0.12134
nbioav 2.7333e7
σ 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.910971 | 0.910988 |
| 2 | θkaSlow | 0.131115 | 0.131117 |
| 3 | θCL | 1.0854 | 1.0854 |
| 4 | θVc | 7.10076 | 7.10073 |
| 5 | θbioav | 0.480202 | 0.480191 |
| 6 | ωCL | 0.0881132 | 0.0881137 |
| 7 | ωVc | 0.121334 | 0.121335 |
| 8 | σ | 0.105448 | 0.105449 |
| 9 | ξbioav | 1.84292e-5 | missing |
| 10 | nbioav | missing | 2.73333e7 |
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.