Why are non-Gaussian random effects relevant?

Authors

Jose Storopoli

Andreas Noack

using Dates
using Pumas
using PumasUtilities
using DataFramesMeta
using PharmaDatasets
using CairoMakie
using AlgebraOfGraphics
using Random

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
end

If 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
end

As 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),)
end

The 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.7

We’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: dv

4 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
end
PumasModel
  Parameters: θCL, θVc, ωCL, ωVc, σ
  Random effects: _CL, _Vc
  Covariates: 
  Dynamical system variables: Central
  Dynamical system type: Closed form
  Derived: dv
  Observed: dv
model_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
end
PumasModel
  Parameters: θCL, θVc, ωCL, ωVc, σ
  Random effects: _CL, _Vc
  Covariates: 
  Dynamical system variables: Central
  Dynamical system type: Matrix exponential
  Derived: dv
  Observed: dv

We 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.02602100372314453
     1     9.433464e+02     6.079483e+02
 * time: 1.0018219947814941
     2     8.189627e+02     4.423725e+02
 * time: 1.0778021812438965
     3     5.917683e+02     1.819248e+02
 * time: 1.1044740676879883
     4     5.421783e+02     1.121313e+02
 * time: 1.125730037689209
     5     5.255651e+02     7.407230e+01
 * time: 1.144388198852539
     6     5.208427e+02     8.699271e+01
 * time: 1.1616590023040771
     7     5.174883e+02     8.974584e+01
 * time: 1.1784040927886963
     8     5.138523e+02     7.328235e+01
 * time: 1.195033073425293
     9     5.109883e+02     4.155805e+01
 * time: 1.2110869884490967
    10     5.094359e+02     3.170517e+01
 * time: 1.225978136062622
    11     5.086172e+02     3.327331e+01
 * time: 1.2407701015472412
    12     5.080941e+02     2.942077e+01
 * time: 1.2558159828186035
    13     5.074009e+02     2.839941e+01
 * time: 1.2705881595611572
    14     5.059302e+02     3.330093e+01
 * time: 1.2850661277770996
    15     5.036399e+02     3.172884e+01
 * time: 1.30059814453125
    16     5.017004e+02     3.160020e+01
 * time: 1.3166542053222656
    17     5.008553e+02     2.599524e+01
 * time: 1.360314130783081
    18     5.005913e+02     2.139314e+01
 * time: 1.375521183013916
    19     5.003573e+02     2.134778e+01
 * time: 1.3901910781860352
    20     4.997249e+02     2.069868e+01
 * time: 1.405493974685669
    21     4.984453e+02     1.859010e+01
 * time: 1.4211821556091309
    22     4.959584e+02     2.156209e+01
 * time: 1.4371252059936523
    23     4.923347e+02     3.030833e+01
 * time: 1.4541740417480469
    24     4.906916e+02     1.652278e+01
 * time: 1.4722859859466553
    25     4.902955e+02     6.360800e+00
 * time: 1.4892420768737793
    26     4.902870e+02     7.028603e+00
 * time: 1.5114049911499023
    27     4.902193e+02     1.176895e+00
 * time: 1.5289530754089355
    28     4.902189e+02     1.170642e+00
 * time: 1.54256010055542
    29     4.902186e+02     1.167624e+00
 * time: 1.5540320873260498
    30     4.902145e+02     1.110377e+00
 * time: 1.5685679912567139
    31     4.902079e+02     1.010507e+00
 * time: 1.582970142364502
    32     4.901917e+02     9.619218e-01
 * time: 1.5974969863891602
    33     4.901683e+02     1.001006e+00
 * time: 1.6227531433105469
    34     4.901473e+02     6.138233e-01
 * time: 1.6379129886627197
    35     4.901412e+02     1.754342e-01
 * time: 1.6530921459197998
    36     4.901406e+02     2.617009e-02
 * time: 1.6666510105133057
    37     4.901405e+02     4.585882e-03
 * time: 1.6785120964050293
    38     4.901405e+02     7.668184e-04
 * time: 1.6881821155548096
FittedPumasModel

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: 6.318092346191406e-5
     1     9.056373e+02     5.541920e+02
 * time: 0.03891611099243164
     2     7.787963e+02     3.988904e+02
 * time: 0.0725240707397461
     3     5.915054e+02     1.495777e+02
 * time: 0.09677004814147949
     4     5.533120e+02     8.826583e+01
 * time: 0.11900115013122559
     5     5.389239e+02     9.144086e+01
 * time: 0.14025497436523438
     6     5.323499e+02     1.000933e+02
 * time: 0.16268706321716309
     7     5.270252e+02     8.423844e+01
 * time: 0.18620514869689941
     8     5.233813e+02     5.194402e+01
 * time: 0.20844817161560059
     9     5.213366e+02     3.461331e+01
 * time: 0.2733621597290039
    10     5.200972e+02     3.888113e+01
 * time: 0.2935221195220947
    11     5.191933e+02     3.556605e+01
 * time: 0.31372809410095215
    12     5.181335e+02     3.624436e+01
 * time: 0.33366918563842773
    13     5.161626e+02     4.322775e+01
 * time: 0.3528480529785156
    14     5.133202e+02     3.722515e+01
 * time: 0.37081408500671387
    15     5.107758e+02     3.401586e+01
 * time: 0.38916015625
    16     5.095157e+02     2.854997e+01
 * time: 0.4082601070404053
    17     5.090165e+02     2.644560e+01
 * time: 0.42623114585876465
    18     5.085184e+02     2.744429e+01
 * time: 0.44341015815734863
    19     5.074309e+02     2.793918e+01
 * time: 0.4605841636657715
    20     5.053757e+02     2.616169e+01
 * time: 0.47840213775634766
    21     5.018507e+02     2.257667e+01
 * time: 0.4974551200866699
    22     4.942495e+02     3.832878e+01
 * time: 0.519212007522583
    23     4.940229e+02     5.518159e+01
 * time: 0.5466179847717285
    24     4.909110e+02     3.042064e+01
 * time: 0.5746350288391113
    25     4.900234e+02     6.929306e+00
 * time: 0.6005170345306396
    26     4.897974e+02     1.087865e+00
 * time: 0.6263010501861572
    27     4.897942e+02     6.456402e-01
 * time: 0.6489071846008301
    28     4.897940e+02     6.467689e-01
 * time: 0.6981141567230225
    29     4.897939e+02     6.463480e-01
 * time: 0.7166881561279297
    30     4.897935e+02     6.408914e-01
 * time: 0.7354199886322021
    31     4.897924e+02     6.208208e-01
 * time: 0.7547152042388916
    32     4.897900e+02     1.035462e+00
 * time: 0.7742931842803955
    33     4.897850e+02     1.452099e+00
 * time: 0.7939980030059814
    34     4.897776e+02     1.482593e+00
 * time: 0.8134751319885254
    35     4.897718e+02     8.420646e-01
 * time: 0.8339540958404541
    36     4.897702e+02     2.023876e-01
 * time: 0.8545091152191162
    37     4.897700e+02     1.885486e-02
 * time: 0.873035192489624
    38     4.897700e+02     2.343932e-03
 * time: 0.8896360397338867
    39     4.897700e+02     4.417566e-04
 * time: 0.9039561748504639
FittedPumasModel

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
end
PumasModel
  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: dv

Now 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
end
PumasModel
  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: dv

We 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.0
trueparam = (;
    θ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: 7.581710815429688e-5
     1     2.005689e+02     1.252637e+02
 * time: 0.3474099636077881
     2     1.875728e+02     4.225518e+01
 * time: 0.5161159038543701
     3     1.863803e+02     3.429785e+01
 * time: 0.6789588928222656
     4     1.845581e+02     3.694229e+01
 * time: 0.868034839630127
     5     1.828416e+02     1.170915e+01
 * time: 1.0379559993743896
     6     1.823277e+02     1.004017e+01
 * time: 1.1548748016357422
     7     1.810629e+02     5.326780e+00
 * time: 1.2712910175323486
     8     1.810479e+02     1.767987e+00
 * time: 1.385589838027954
     9     1.810272e+02     3.852037e+00
 * time: 1.509922981262207
    10     1.809414e+02     1.196625e+01
 * time: 1.6359548568725586
    11     1.806489e+02     1.861440e+01
 * time: 1.7637288570404053
    12     1.801984e+02     1.361774e+01
 * time: 2.254415988922119
    13     1.800427e+02     2.177039e+01
 * time: 2.3771588802337646
    14     1.796554e+02     7.079039e+00
 * time: 2.504149913787842
    15     1.795832e+02     1.581499e+01
 * time: 2.633692979812622
    16     1.795220e+02     6.552120e+00
 * time: 2.768831968307495
    17     1.794621e+02     6.612645e+00
 * time: 2.89694881439209
    18     1.793643e+02     6.603903e+00
 * time: 3.0205190181732178
    19     1.793502e+02     1.025787e+01
 * time: 3.165799856185913
    20     1.793178e+02     1.202199e+01
 * time: 3.2765138149261475
    21     1.792424e+02     1.625661e+01
 * time: 3.3882968425750732
    22     1.791560e+02     1.128856e+01
 * time: 3.5046188831329346
    23     1.790991e+02     7.625637e+00
 * time: 3.6228508949279785
    24     1.790756e+02     8.557258e+00
 * time: 3.740238904953003
    25     1.790633e+02     4.848482e+00
 * time: 3.859757900238037
    26     1.790408e+02     5.859306e+00
 * time: 3.978538990020752
    27     1.789734e+02     1.106035e+01
 * time: 4.121018886566162
    28     1.789070e+02     1.143361e+01
 * time: 4.231564998626709
    29     1.788321e+02     7.098448e+00
 * time: 4.341272830963135
    30     1.787896e+02     5.709857e+00
 * time: 4.454293966293335
    31     1.787809e+02     7.456555e+00
 * time: 4.573001861572266
    32     1.787671e+02     7.320931e-01
 * time: 4.692580938339233
    33     1.787660e+02     5.023990e-01
 * time: 4.810844898223877
    34     1.787656e+02     3.813994e-01
 * time: 4.925073862075806
    35     1.787639e+02     8.608909e-01
 * time: 5.069129943847656
    36     1.787607e+02     2.047321e+00
 * time: 5.17998480796814
    37     1.787525e+02     3.859529e+00
 * time: 5.289959907531738
    38     1.787349e+02     5.864920e+00
 * time: 5.413597822189331
    39     1.787017e+02     6.966045e+00
 * time: 5.548888921737671
    40     1.786574e+02     5.348939e+00
 * time: 5.685351848602295
    41     1.786270e+02     1.642025e+00
 * time: 5.811041831970215
    42     1.786181e+02     5.211823e-01
 * time: 5.934624910354614
    43     1.786153e+02     1.186061e+00
 * time: 6.062311887741089
    44     1.786125e+02     1.292005e+00
 * time: 6.205183982849121
    45     1.786099e+02     7.814598e-01
 * time: 6.315224885940552
    46     1.786086e+02     1.369456e-01
 * time: 6.42535400390625
    47     1.786082e+02     1.912170e-01
 * time: 6.535862922668457
    48     1.786080e+02     2.670802e-01
 * time: 6.654233932495117
    49     1.786078e+02     1.979262e-01
 * time: 6.768675804138184
    50     1.786077e+02     5.177918e-02
 * time: 6.885594844818115
    51     1.786076e+02     2.998328e-02
 * time: 7.003453969955444
    52     1.786076e+02     5.799706e-02
 * time: 7.1177818775177
    53     1.786076e+02     4.280434e-02
 * time: 7.259394884109497
    54     1.786076e+02     1.467829e-02
 * time: 7.3675549030303955
    55     1.786076e+02     6.928178e-03
 * time: 7.468969821929932
    56     1.786076e+02     1.236015e-02
 * time: 7.573523998260498
    57     1.786076e+02     9.950826e-03
 * time: 7.686716794967651
    58     1.786076e+02     2.974370e-03
 * time: 7.793972969055176
    59     1.786076e+02     1.374576e-03
 * time: 7.90009880065918
    60     1.786076e+02     2.860078e-03
 * time: 7.998617887496948
    61     1.786076e+02     2.100127e-03
 * time: 8.097692966461182
    62     1.786076e+02     2.100127e-03
 * time: 8.266541957855225
    63     1.786076e+02     6.412834e-03
 * time: 8.383128881454468
    64     1.786076e+02     6.412834e-03
 * time: 8.541384935379028
    65     1.786076e+02     6.412834e-03
 * time: 8.77962589263916
    66     1.786076e+02     6.412834e-03
 * time: 9.194366931915283
    67     1.786076e+02     6.412834e-03
 * time: 9.422312021255493
FittedPumasModel

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: 6.103515625e-5
     1     2.014577e+02     1.265001e+02
 * time: 0.3875720500946045
     2     1.880885e+02     4.190708e+01
 * time: 0.5650589466094971
     3     1.870317e+02     8.825666e+01
 * time: 0.7014479637145996
     4     1.846027e+02     4.400156e+01
 * time: 0.9199090003967285
     5     1.834445e+02     1.906624e+01
 * time: 1.0387461185455322
     6     1.828599e+02     1.113882e+01
 * time: 1.15706205368042
     7     1.815719e+02     7.449355e+00
 * time: 1.280817985534668
     8     1.815131e+02     2.164678e+00
 * time: 1.4090170860290527
     9     1.814896e+02     2.167319e+00
 * time: 1.5333330631256104
    10     1.814458e+02     4.615738e+00
 * time: 1.6596269607543945
    11     1.813173e+02     9.576967e+00
 * time: 1.7890710830688477
    12     1.809756e+02     2.052077e+01
 * time: 1.9889180660247803
    13     1.807625e+02     4.553366e+01
 * time: 2.1098780632019043
    14     1.802224e+02     6.550892e+00
 * time: 2.2313129901885986
    15     1.800862e+02     2.865509e+00
 * time: 2.350567102432251
    16     1.800780e+02     1.164611e+00
 * time: 2.4725489616394043
    17     1.800737e+02     7.952462e-01
 * time: 2.600260019302368
    18     1.800352e+02     4.860618e+00
 * time: 2.7309439182281494
    19     1.800089e+02     5.176689e+00
 * time: 2.863070011138916
    20     1.799679e+02     4.303892e+00
 * time: 2.990499973297119
    21     1.799153e+02     4.612832e+00
 * time: 3.1363110542297363
    22     1.798423e+02     1.209387e+01
 * time: 3.249553918838501
    23     1.796821e+02     1.712256e+01
 * time: 3.360459089279175
    24     1.794275e+02     1.435100e+01
 * time: 3.473684072494507
    25     1.793773e+02     4.137313e+00
 * time: 3.5931200981140137
    26     1.793488e+02     1.846545e+00
 * time: 3.745600938796997
    27     1.793331e+02     5.502533e+00
 * time: 3.865586996078491
    28     1.793234e+02     2.894037e+00
 * time: 3.9850361347198486
    29     1.793119e+02     1.453372e+00
 * time: 4.131155967712402
    30     1.792879e+02     4.109884e+00
 * time: 4.242964029312134
    31     1.792599e+02     4.613173e+00
 * time: 4.359745025634766
    32     1.792384e+02     4.254549e+00
 * time: 4.469506025314331
    33     1.792251e+02     4.415024e+00
 * time: 4.589903116226196
    34     1.792026e+02     3.229145e+00
 * time: 4.709100961685181
    35     1.791841e+02     3.422094e+00
 * time: 4.8633880615234375
    36     1.791705e+02     1.621228e+00
 * time: 5.015721082687378
    37     1.791555e+02     3.462667e+00
 * time: 5.143255949020386
    38     1.791293e+02     5.586685e+00
 * time: 5.28898811340332
    39     1.790713e+02     9.927638e+00
 * time: 5.483910083770752
    40     1.789891e+02     1.133271e+01
 * time: 5.598509073257446
    41     1.788585e+02     1.279172e+01
 * time: 5.716385126113892
    42     1.787407e+02     5.291681e+00
 * time: 5.8444530963897705
    43     1.786633e+02     5.367971e+00
 * time: 5.975500106811523
    44     1.786405e+02     2.571548e+00
 * time: 6.099683046340942
    45     1.786339e+02     2.720314e+00
 * time: 6.22140908241272
    46     1.786252e+02     1.930583e+00
 * time: 6.37310004234314
    47     1.786182e+02     1.523164e+00
 * time: 6.483062028884888
    48     1.786126e+02     5.027920e-01
 * time: 6.592617034912109
    49     1.786100e+02     3.733023e-01
 * time: 6.7067201137542725
    50     1.786089e+02     2.749553e-01
 * time: 6.823368072509766
    51     1.786083e+02     1.981772e-01
 * time: 6.9436609745025635
    52     1.786079e+02     1.005262e-01
 * time: 7.060589075088501
    53     1.786078e+02     2.549641e-02
 * time: 7.177897930145264
    54     1.786077e+02     4.452931e-02
 * time: 7.293743133544922
    55     1.786076e+02     2.967125e-02
 * time: 7.407896041870117
    56     1.786076e+02     1.068274e-02
 * time: 7.5468199253082275
    57     1.786076e+02     5.355447e-03
 * time: 7.65038800239563
    58     1.786076e+02     8.180920e-03
 * time: 7.752283096313477
    59     1.786076e+02     4.935439e-03
 * time: 7.855996131896973
    60     1.786076e+02     4.387986e-03
 * time: 7.992188930511475
    61     1.786076e+02     4.388023e-03
 * time: 8.174855947494507
    62     1.786076e+02     4.387934e-03
 * time: 8.339471101760864
    63     1.786076e+02     4.387934e-03
 * time: 8.506487131118774
FittedPumasModel

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.