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.022639989852905273
1 9.433464e+02 6.079483e+02
* time: 0.4672970771789551
2 8.189627e+02 4.423725e+02
* time: 0.5423541069030762
3 5.917683e+02 1.819248e+02
* time: 0.5495471954345703
4 5.421783e+02 1.121313e+02
* time: 0.5554101467132568
5 5.255651e+02 7.407230e+01
* time: 0.5607972145080566
6 5.208427e+02 8.699271e+01
* time: 0.5662209987640381
7 5.174883e+02 8.974584e+01
* time: 0.571274995803833
8 5.138523e+02 7.328235e+01
* time: 0.5761439800262451
9 5.109883e+02 4.155805e+01
* time: 0.5803589820861816
10 5.094359e+02 3.170517e+01
* time: 0.5847651958465576
11 5.086172e+02 3.327331e+01
* time: 0.5893440246582031
12 5.080941e+02 2.942077e+01
* time: 0.5940890312194824
13 5.074009e+02 2.839941e+01
* time: 0.6293261051177979
14 5.059302e+02 3.330093e+01
* time: 0.6335711479187012
15 5.036399e+02 3.172884e+01
* time: 0.6380481719970703
16 5.017004e+02 3.160020e+01
* time: 0.6425681114196777
17 5.008553e+02 2.599524e+01
* time: 0.6470451354980469
18 5.005913e+02 2.139314e+01
* time: 0.6515271663665771
19 5.003573e+02 2.134778e+01
* time: 0.6560111045837402
20 4.997249e+02 2.069868e+01
* time: 0.6606221199035645
21 4.984453e+02 1.859010e+01
* time: 0.6651790142059326
22 4.959584e+02 2.156209e+01
* time: 0.6696720123291016
23 4.923347e+02 3.030833e+01
* time: 0.6744470596313477
24 4.906916e+02 1.652278e+01
* time: 0.6794600486755371
25 4.902955e+02 6.360800e+00
* time: 0.6848950386047363
26 4.902870e+02 7.028603e+00
* time: 0.7130320072174072
27 4.902193e+02 1.176895e+00
* time: 0.7176849842071533
28 4.902189e+02 1.170642e+00
* time: 0.7214779853820801
29 4.902186e+02 1.167624e+00
* time: 0.7249040603637695
30 4.902145e+02 1.110377e+00
* time: 0.7290530204772949
31 4.902079e+02 1.010507e+00
* time: 0.7331171035766602
32 4.901917e+02 9.619218e-01
* time: 0.7372860908508301
33 4.901683e+02 1.001006e+00
* time: 0.7411642074584961
34 4.901473e+02 6.138233e-01
* time: 0.7456061840057373
35 4.901412e+02 1.754342e-01
* time: 0.749920129776001
36 4.901406e+02 2.617009e-02
* time: 0.7538762092590332
37 4.901405e+02 4.585882e-03
* time: 0.7576210498809814
38 4.901405e+02 7.668184e-04
* time: 0.7610540390014648FittedPumasModel
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.219519e+03 5.960578e+03
* time: 5.2928924560546875e-5
1 9.644493e+02 5.541920e+02
* time: 0.01297903060913086
2 8.376083e+02 3.988904e+02
* time: 0.022876977920532227
3 6.503174e+02 1.495777e+02
* time: 0.03140592575073242
4 6.121241e+02 8.826583e+01
* time: 0.07075190544128418
5 5.977360e+02 9.144086e+01
* time: 0.07838606834411621
6 5.911620e+02 1.000933e+02
* time: 0.08572793006896973
7 5.858372e+02 8.423844e+01
* time: 0.09274697303771973
8 5.821934e+02 5.194402e+01
* time: 0.09958600997924805
9 5.801487e+02 3.461331e+01
* time: 0.10676884651184082
10 5.789092e+02 3.888113e+01
* time: 0.11370301246643066
11 5.780054e+02 3.556605e+01
* time: 0.12068986892700195
12 5.769455e+02 3.624436e+01
* time: 0.12772703170776367
13 5.749747e+02 4.322775e+01
* time: 0.1347949504852295
14 5.721322e+02 3.722515e+01
* time: 0.1415410041809082
15 5.695879e+02 3.401586e+01
* time: 0.14921092987060547
16 5.683277e+02 2.854997e+01
* time: 0.18210697174072266
17 5.678285e+02 2.644560e+01
* time: 0.1887669563293457
18 5.673305e+02 2.744429e+01
* time: 0.1951289176940918
19 5.662430e+02 2.793918e+01
* time: 0.20142698287963867
20 5.641877e+02 2.616169e+01
* time: 0.20791196823120117
21 5.606628e+02 2.257667e+01
* time: 0.21489691734313965
22 5.530616e+02 3.832878e+01
* time: 0.22240304946899414
23 5.528349e+02 5.518159e+01
* time: 0.2317349910736084
24 5.497231e+02 3.042064e+01
* time: 0.2405259609222412
25 5.488355e+02 6.929306e+00
* time: 0.2485489845275879
26 5.486095e+02 1.087865e+00
* time: 0.2574489116668701
27 5.486062e+02 6.456402e-01
* time: 0.28629088401794434
28 5.486061e+02 6.467689e-01
* time: 0.2930169105529785
29 5.486060e+02 6.463480e-01
* time: 0.29962992668151855
30 5.486055e+02 6.408914e-01
* time: 0.30635499954223633
31 5.486045e+02 6.208208e-01
* time: 0.31331491470336914
32 5.486020e+02 1.035462e+00
* time: 0.3204069137573242
33 5.485971e+02 1.452099e+00
* time: 0.32751893997192383
34 5.485897e+02 1.482593e+00
* time: 0.33457303047180176
35 5.485839e+02 8.420646e-01
* time: 0.3418459892272949
36 5.485822e+02 2.023876e-01
* time: 0.34912705421447754
37 5.485821e+02 1.885486e-02
* time: 0.3560800552368164
38 5.485821e+02 2.343932e-03
* time: 0.37398290634155273
39 5.485821e+02 4.417566e-04
* time: 0.3797950744628906FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Matrix exponential
Log-likelihood value: -548.58208
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):
simpop =
Subject.(
simobs(
model_beta,
skeletonpop,
trueparam;
obstimes = simtimes,
rng = Random.seed!(Random.default_rng(), 128),
)
)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 3.179853e+02 3.723266e+02
* time: 6.818771362304688e-5
1 2.618605e+02 1.558889e+02
* time: 0.230363130569458
2 2.484481e+02 5.652093e+01
* time: 0.32152414321899414
3 2.466691e+02 3.819960e+01
* time: 0.4127671718597412
4 2.430995e+02 2.281161e+01
* time: 0.48877406120300293
5 2.421492e+02 1.230731e+01
* time: 0.5672860145568848
6 2.412891e+02 1.022578e+01
* time: 0.6767940521240234
7 2.398220e+02 4.018354e+00
* time: 0.7541241645812988
8 2.398069e+02 2.747661e+00
* time: 0.8262591361999512
9 2.398011e+02 1.787833e+00
* time: 0.8992671966552734
10 2.397895e+02 2.923095e+00
* time: 0.9753530025482178
11 2.397511e+02 4.766103e+00
* time: 1.0540721416473389
12 2.397274e+02 2.427078e+00
* time: 1.1599199771881104
13 2.396969e+02 2.317178e+00
* time: 1.233046054840088
14 2.396432e+02 6.390959e+00
* time: 1.306913137435913
15 2.395437e+02 1.431396e+01
* time: 1.3838679790496826
16 2.394242e+02 1.660862e+01
* time: 1.465512990951538
17 2.392883e+02 7.918201e+00
* time: 1.549138069152832
18 2.392648e+02 5.359476e+00
* time: 1.6643791198730469
19 2.392551e+02 8.504104e-01
* time: 1.7406351566314697
20 2.392536e+02 8.658366e-01
* time: 1.8199529647827148
21 2.392502e+02 1.359491e+00
* time: 1.902726173400879
22 2.392455e+02 1.908257e+00
* time: 1.986847162246704
23 2.392352e+02 3.434282e+00
* time: 2.073047161102295
24 2.392094e+02 4.019057e+00
* time: 2.186134099960327
25 2.391904e+02 2.360120e+00
* time: 2.282292127609253
26 2.391639e+02 2.232625e+00
* time: 2.3794641494750977
27 2.390985e+02 7.831580e+00
* time: 2.4656729698181152
28 2.390462e+02 1.220982e+01
* time: 2.5519449710845947
29 2.389849e+02 1.372374e+01
* time: 2.662506103515625
30 2.389295e+02 1.012388e+01
* time: 2.7423291206359863
31 2.389050e+02 7.364830e+00
* time: 2.8356921672821045
32 2.388675e+02 3.493531e+00
* time: 2.919481039047241
33 2.388348e+02 3.601972e+00
* time: 3.0026121139526367
34 2.388092e+02 9.896791e+00
* time: 3.099238157272339
35 2.387816e+02 5.816471e+00
* time: 3.1775050163269043
36 2.387667e+02 3.995305e+00
* time: 3.255885124206543
37 2.387624e+02 2.111619e+00
* time: 3.3375420570373535
38 2.387542e+02 3.340445e+00
* time: 3.4189090728759766
39 2.387468e+02 3.742750e+00
* time: 3.4994029998779297
40 2.387148e+02 8.270335e+00
* time: 3.5891401767730713
41 2.386920e+02 8.042091e+00
* time: 3.661177158355713
42 2.386311e+02 1.786257e+00
* time: 3.7374050617218018
43 2.385982e+02 4.295253e+00
* time: 3.8166260719299316
44 2.385769e+02 2.426177e+00
* time: 3.8975741863250732
45 2.385642e+02 4.507279e-01
* time: 3.979091167449951
46 2.385598e+02 1.783922e+00
* time: 4.068833112716675
47 2.385584e+02 4.345023e-01
* time: 4.142912149429321
48 2.385582e+02 3.504116e-01
* time: 4.216340065002441
49 2.385581e+02 9.328221e-02
* time: 4.293819189071655
50 2.385580e+02 1.265801e-02
* time: 4.370881080627441
51 2.385580e+02 2.371756e-03
* time: 4.4463629722595215
52 2.385580e+02 2.371756e-03
* time: 4.548573017120361
53 2.385580e+02 2.371756e-03
* time: 4.646895170211792
54 2.385580e+02 2.371756e-03
* time: 4.7547080516815186
55 2.385580e+02 2.371756e-03
* time: 4.873974084854126
56 2.385580e+02 2.371756e-03
* time: 4.974316120147705
57 2.385580e+02 2.371756e-03
* time: 5.073625087738037
58 2.385580e+02 2.371756e-03
* time: 5.117396116256714FittedPumasModel
Successful minimization: true
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Matrix exponential
Log-likelihood value: -238.55805
Number of subjects: 40
Number of parameters: Fixed Optimized
0 9
Observation records: Active Missing
dv: 240 0
Total: 240 0
---------------------
Estimate
---------------------
θkaFast 1.8907
θkaSlow 0.60771
θCL 1.0797
θVc 11.375
θbioav 0.14394
ωCL 0.08151
ωVc 0.10157
ξbioav 0.12513
σ 0.10441
---------------------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 3.554658e+02 3.718530e+02
* time: 5.412101745605469e-5
1 2.991160e+02 1.342892e+02
* time: 0.2258930206298828
2 2.851580e+02 4.511285e+01
* time: 0.3876941204071045
3 2.835992e+02 4.343268e+01
* time: 0.4952681064605713
4 2.821916e+02 4.373720e+01
* time: 0.5836749076843262
5 2.794524e+02 1.253857e+01
* time: 0.6774489879608154
6 2.787960e+02 1.102093e+01
* time: 0.8085079193115234
7 2.770092e+02 1.114240e+01
* time: 0.8974440097808838
8 2.769202e+02 5.089120e+00
* time: 0.9857780933380127
9 2.769088e+02 3.041932e+00
* time: 1.0754330158233643
10 2.768987e+02 1.603805e+00
* time: 1.164910078048706
11 2.767766e+02 1.411814e+01
* time: 1.260329008102417
12 2.766170e+02 2.719464e+01
* time: 1.380458116531372
13 2.764580e+02 3.014173e+01
* time: 1.469231128692627
14 2.762515e+02 1.669149e+01
* time: 1.5588390827178955
15 2.761434e+02 1.957297e+00
* time: 1.6520869731903076
16 2.761353e+02 3.756442e+00
* time: 1.7460319995880127
17 2.761237e+02 1.341491e+00
* time: 1.853945016860962
18 2.761152e+02 2.182772e+00
* time: 1.9394879341125488
19 2.761036e+02 3.209881e+00
* time: 2.0262770652770996
20 2.760912e+02 3.070991e+00
* time: 2.1159961223602295
21 2.760762e+02 1.990495e+00
* time: 2.2077889442443848
22 2.760587e+02 1.586627e+00
* time: 2.2998900413513184
23 2.760290e+02 2.471592e+00
* time: 2.4040980339050293
24 2.759737e+02 1.040809e+01
* time: 2.491948127746582
25 2.758948e+02 9.671126e+00
* time: 2.5789389610290527
26 2.758264e+02 7.207553e+00
* time: 2.6710479259490967
27 2.757728e+02 3.300455e+00
* time: 2.800842046737671
28 2.757344e+02 4.723183e+00
* time: 2.884537935256958
29 2.756731e+02 6.515451e+00
* time: 2.970367908477783
30 2.755936e+02 4.181865e+00
* time: 3.0588579177856445
31 2.755684e+02 5.163221e+00
* time: 3.169851064682007
32 2.755454e+02 8.132621e+00
* time: 3.278921127319336
33 2.755346e+02 2.188011e+00
* time: 3.363079071044922
34 2.755279e+02 3.127135e+00
* time: 3.4465129375457764
35 2.755221e+02 3.813163e+00
* time: 3.530811071395874
36 2.754865e+02 6.610197e+00
* time: 3.6212151050567627
37 2.754663e+02 5.265615e+00
* time: 3.711851119995117
38 2.754486e+02 9.998959e-01
* time: 3.817291021347046
39 2.754421e+02 1.112426e+00
* time: 3.9027609825134277
40 2.754395e+02 3.014869e-01
* time: 3.989914894104004
41 2.754386e+02 5.351895e-01
* time: 4.080051898956299
42 2.754384e+02 1.402404e-01
* time: 4.169112920761108
43 2.754383e+02 1.288533e-01
* time: 4.27583909034729
44 2.754381e+02 3.845371e-01
* time: 4.360318899154663
45 2.754376e+02 9.172184e-01
* time: 4.44579291343689
46 2.754363e+02 1.766701e+00
* time: 4.5349860191345215
47 2.754331e+02 3.007560e+00
* time: 4.630875110626221
48 2.754256e+02 4.573254e+00
* time: 4.726067066192627
49 2.754117e+02 5.874317e+00
* time: 4.834105968475342
50 2.753921e+02 6.633185e+00
* time: 4.921669960021973
51 2.753744e+02 2.205111e+00
* time: 5.011238098144531
52 2.753702e+02 2.363930e+00
* time: 5.104016065597534
53 2.753546e+02 8.713730e-01
* time: 5.198376893997192
54 2.753543e+02 3.795710e+00
* time: 5.307397127151489
55 2.753492e+02 1.784042e+00
* time: 5.392431020736694
56 2.753441e+02 8.514295e-01
* time: 5.474895000457764
57 2.753339e+02 1.341633e+00
* time: 5.562907934188843
58 2.753266e+02 2.387419e+00
* time: 5.653064966201782
59 2.753233e+02 1.246148e+00
* time: 5.740912914276123
60 2.753208e+02 2.973353e-01
* time: 5.839653968811035
61 2.753202e+02 3.696236e-01
* time: 5.917881965637207
62 2.753193e+02 4.553497e-01
* time: 5.999064922332764
63 2.753186e+02 2.081281e-01
* time: 6.083395957946777
64 2.753183e+02 8.672527e-02
* time: 6.167145013809204
65 2.753182e+02 1.015286e-01
* time: 6.263946056365967
66 2.753181e+02 1.003887e-01
* time: 6.339065074920654
67 2.753181e+02 5.060689e-02
* time: 6.4188148975372314
68 2.753180e+02 2.051541e-02
* time: 6.497474908828735
69 2.753180e+02 1.864290e-02
* time: 6.575284957885742
70 2.753180e+02 1.911113e-02
* time: 6.662434101104736
71 2.753180e+02 8.991085e-03
* time: 6.765403985977173
72 2.753180e+02 4.925215e-03
* time: 6.845324993133545
73 2.753180e+02 2.135748e-03
* time: 6.925084114074707
74 2.753180e+02 2.135926e-03
* time: 7.047389984130859FittedPumasModel
Successful minimization: false
Likelihood approximation: FOCE
Likelihood Optimizer: BFGS
Dynamical system type: Matrix exponential
Log-likelihood value: -275.31802
Number of subjects: 40
Number of parameters: Fixed Optimized
0 9
Observation records: Active Missing
dv: 240 0
Total: 240 0
----------------------
Estimate
----------------------
θkaFast 1.9882
θkaSlow 0.61333
θCL 1.0797
θVc 11.381
θbioav 0.13274
ωCL 0.081471
ωVc 0.10178
nbioav 1.4823e7
σ 0.10464
----------------------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 | 1.89067 | 1.98819 |
| 2 | θkaSlow | 0.607708 | 0.613333 |
| 3 | θCL | 1.07968 | 1.07967 |
| 4 | θVc | 11.3745 | 11.3811 |
| 5 | θbioav | 0.143942 | 0.132743 |
| 6 | ωCL | 0.0815101 | 0.081471 |
| 7 | ωVc | 0.101565 | 0.101784 |
| 8 | σ | 0.10441 | 0.104644 |
| 9 | ξbioav | 0.125128 | missing |
| 10 | nbioav | missing | 1.48228e7 |
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 | 1.54002e-269 | 0.0 |
| 4 | 0.003003 | 4.49003e-222 | 0.0 |
| 5 | 0.004004 | 3.80322e-191 | 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.