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.03019094467163086
     1     9.433464e+02     6.079483e+02
 * time: 1.353376865386963
     2     8.189627e+02     4.423725e+02
 * time: 1.3939549922943115
     3     5.917683e+02     1.819248e+02
 * time: 1.4251210689544678
     4     5.421783e+02     1.121313e+02
 * time: 1.4504680633544922
     5     5.255651e+02     7.407230e+01
 * time: 1.4725499153137207
     6     5.208427e+02     8.699271e+01
 * time: 1.492473840713501
     7     5.174883e+02     8.974584e+01
 * time: 1.509335994720459
     8     5.138523e+02     7.328235e+01
 * time: 1.526608943939209
     9     5.109883e+02     4.155805e+01
 * time: 1.542827844619751
    10     5.094359e+02     3.170517e+01
 * time: 1.5589170455932617
    11     5.086172e+02     3.327331e+01
 * time: 1.5766479969024658
    12     5.080941e+02     2.942077e+01
 * time: 1.5954580307006836
    13     5.074009e+02     2.839941e+01
 * time: 1.614182949066162
    14     5.059302e+02     3.330093e+01
 * time: 1.7384769916534424
    15     5.036399e+02     3.172884e+01
 * time: 1.7559590339660645
    16     5.017004e+02     3.160020e+01
 * time: 1.7732889652252197
    17     5.008553e+02     2.599524e+01
 * time: 1.7898850440979004
    18     5.005913e+02     2.139314e+01
 * time: 1.8053419589996338
    19     5.003573e+02     2.134778e+01
 * time: 1.820483922958374
    20     4.997249e+02     2.069868e+01
 * time: 1.8363189697265625
    21     4.984453e+02     1.859010e+01
 * time: 1.8533010482788086
    22     4.959584e+02     2.156209e+01
 * time: 1.870826005935669
    23     4.923347e+02     3.030833e+01
 * time: 1.888974905014038
    24     4.906916e+02     1.652278e+01
 * time: 1.9124529361724854
    25     4.902955e+02     6.360800e+00
 * time: 1.9307138919830322
    26     4.902870e+02     7.028603e+00
 * time: 1.9501080513000488
    27     4.902193e+02     1.176895e+00
 * time: 1.9682848453521729
    28     4.902189e+02     1.170642e+00
 * time: 1.9828920364379883
    29     4.902186e+02     1.167624e+00
 * time: 1.9954710006713867
    30     4.902145e+02     1.110377e+00
 * time: 2.010745048522949
    31     4.902079e+02     1.010507e+00
 * time: 2.0260698795318604
    32     4.901917e+02     9.619218e-01
 * time: 2.041718006134033
    33     4.901683e+02     1.001006e+00
 * time: 2.0571420192718506
    34     4.901473e+02     6.138233e-01
 * time: 2.0761020183563232
    35     4.901412e+02     1.754342e-01
 * time: 2.095212936401367
    36     4.901406e+02     2.617009e-02
 * time: 2.1125588417053223
    37     4.901405e+02     4.585882e-03
 * time: 2.1276509761810303
    38     4.901405e+02     7.668184e-04
 * time: 2.1401848793029785
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: 7.605552673339844e-5
     1     9.056373e+02     5.541920e+02
 * time: 0.04127001762390137
     2     7.787963e+02     3.988904e+02
 * time: 0.0713341236114502
     3     5.915054e+02     1.495777e+02
 * time: 0.09560513496398926
     4     5.533120e+02     8.826583e+01
 * time: 0.11838912963867188
     5     5.389239e+02     9.144086e+01
 * time: 0.14207100868225098
     6     5.323499e+02     1.000933e+02
 * time: 0.16589617729187012
     7     5.270252e+02     8.423844e+01
 * time: 0.18761801719665527
     8     5.233813e+02     5.194402e+01
 * time: 0.20717597007751465
     9     5.213366e+02     3.461331e+01
 * time: 0.22741103172302246
    10     5.200972e+02     3.888113e+01
 * time: 0.24939513206481934
    11     5.191933e+02     3.556605e+01
 * time: 0.2716820240020752
    12     5.181335e+02     3.624436e+01
 * time: 0.29453206062316895
    13     5.161626e+02     4.322775e+01
 * time: 0.3178591728210449
    14     5.133202e+02     3.722515e+01
 * time: 0.3403491973876953
    15     5.107758e+02     3.401586e+01
 * time: 0.36373114585876465
    16     5.095157e+02     2.854997e+01
 * time: 0.3889780044555664
    17     5.090165e+02     2.644560e+01
 * time: 0.414959192276001
    18     5.085184e+02     2.744429e+01
 * time: 0.4393889904022217
    19     5.074309e+02     2.793918e+01
 * time: 0.4628779888153076
    20     5.053757e+02     2.616169e+01
 * time: 0.48591017723083496
    21     5.018507e+02     2.257667e+01
 * time: 0.5973429679870605
    22     4.942495e+02     3.832878e+01
 * time: 0.6196150779724121
    23     4.940229e+02     5.518159e+01
 * time: 0.6472940444946289
    24     4.909110e+02     3.042064e+01
 * time: 0.6733942031860352
    25     4.900234e+02     6.929306e+00
 * time: 0.6974291801452637
    26     4.897974e+02     1.087865e+00
 * time: 0.7206649780273438
    27     4.897942e+02     6.456402e-01
 * time: 0.7410562038421631
    28     4.897940e+02     6.467689e-01
 * time: 0.7590069770812988
    29     4.897939e+02     6.463480e-01
 * time: 0.7779541015625
    30     4.897935e+02     6.408914e-01
 * time: 0.7976460456848145
    31     4.897924e+02     6.208208e-01
 * time: 0.816892147064209
    32     4.897900e+02     1.035462e+00
 * time: 0.8364500999450684
    33     4.897850e+02     1.452099e+00
 * time: 0.8562309741973877
    34     4.897776e+02     1.482593e+00
 * time: 0.8750159740447998
    35     4.897718e+02     8.420646e-01
 * time: 0.8949220180511475
    36     4.897702e+02     2.023876e-01
 * time: 0.9152951240539551
    37     4.897700e+02     1.885486e-02
 * time: 0.9348080158233643
    38     4.897700e+02     2.343932e-03
 * time: 0.9532670974731445
    39     4.897700e+02     4.417566e-04
 * time: 0.9695911407470703
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.414817810058594e-5
     1     2.005689e+02     1.252637e+02
 * time: 0.4500119686126709
     2     1.875728e+02     4.225518e+01
 * time: 0.5945520401000977
     3     1.863803e+02     3.429785e+01
 * time: 0.7343041896820068
     4     1.845581e+02     3.694229e+01
 * time: 0.8472030162811279
     5     1.828416e+02     1.170915e+01
 * time: 1.0117571353912354
     6     1.823277e+02     1.004017e+01
 * time: 1.1256740093231201
     7     1.810629e+02     5.326780e+00
 * time: 1.2429859638214111
     8     1.810479e+02     1.767987e+00
 * time: 1.3514091968536377
     9     1.810272e+02     3.852037e+00
 * time: 1.458204984664917
    10     1.809414e+02     1.196625e+01
 * time: 1.5720210075378418
    11     1.806489e+02     1.861440e+01
 * time: 1.708198070526123
    12     1.801984e+02     1.361774e+01
 * time: 2.1012110710144043
    13     1.800427e+02     2.177039e+01
 * time: 2.2236781120300293
    14     1.796554e+02     7.079039e+00
 * time: 2.3736460208892822
    15     1.795832e+02     1.581499e+01
 * time: 2.4939141273498535
    16     1.795220e+02     6.552120e+00
 * time: 2.6164820194244385
    17     1.794621e+02     6.612645e+00
 * time: 2.7344610691070557
    18     1.793643e+02     6.603903e+00
 * time: 2.852046012878418
    19     1.793502e+02     1.025787e+01
 * time: 2.973180055618286
    20     1.793178e+02     1.202199e+01
 * time: 3.1122090816497803
    21     1.792424e+02     1.625661e+01
 * time: 3.227485179901123
    22     1.791560e+02     1.128856e+01
 * time: 3.3419041633605957
    23     1.790991e+02     7.625637e+00
 * time: 3.456753969192505
    24     1.790756e+02     8.557258e+00
 * time: 3.574429988861084
    25     1.790633e+02     4.848482e+00
 * time: 3.7008211612701416
    26     1.790408e+02     5.859306e+00
 * time: 3.842241048812866
    27     1.789734e+02     1.106035e+01
 * time: 3.970172166824341
    28     1.789070e+02     1.143361e+01
 * time: 4.086485147476196
    29     1.788321e+02     7.098448e+00
 * time: 4.204787015914917
    30     1.787896e+02     5.709857e+00
 * time: 4.33295202255249
    31     1.787809e+02     7.456555e+00
 * time: 4.457854986190796
    32     1.787671e+02     7.320931e-01
 * time: 4.583729028701782
    33     1.787660e+02     5.023990e-01
 * time: 4.6835291385650635
    34     1.787656e+02     3.813994e-01
 * time: 4.789961099624634
    35     1.787639e+02     8.608909e-01
 * time: 4.897998094558716
    36     1.787607e+02     2.047321e+00
 * time: 5.011054992675781
    37     1.787525e+02     3.859529e+00
 * time: 5.121412038803101
    38     1.787349e+02     5.864920e+00
 * time: 5.248653173446655
    39     1.787017e+02     6.966045e+00
 * time: 5.35902214050293
    40     1.786574e+02     5.348939e+00
 * time: 5.476291179656982
    41     1.786270e+02     1.642025e+00
 * time: 5.591240167617798
    42     1.786181e+02     5.211823e-01
 * time: 5.711748123168945
    43     1.786153e+02     1.186061e+00
 * time: 5.831022024154663
    44     1.786125e+02     1.292005e+00
 * time: 5.949275970458984
    45     1.786099e+02     7.814598e-01
 * time: 6.046231031417847
    46     1.786086e+02     1.369456e-01
 * time: 6.159016132354736
    47     1.786082e+02     1.912170e-01
 * time: 6.269965171813965
    48     1.786080e+02     2.670802e-01
 * time: 6.385892152786255
    49     1.786078e+02     1.979262e-01
 * time: 6.5037829875946045
    50     1.786077e+02     5.177918e-02
 * time: 6.645445108413696
    51     1.786076e+02     2.998328e-02
 * time: 6.7523791790008545
    52     1.786076e+02     5.799706e-02
 * time: 6.868740081787109
    53     1.786076e+02     4.280434e-02
 * time: 6.984434127807617
    54     1.786076e+02     1.467829e-02
 * time: 7.107870101928711
    55     1.786076e+02     6.928178e-03
 * time: 7.227673053741455
    56     1.786076e+02     1.236015e-02
 * time: 7.3375630378723145
    57     1.786076e+02     9.950826e-03
 * time: 7.4569830894470215
    58     1.786076e+02     2.974370e-03
 * time: 7.552105188369751
    59     1.786076e+02     1.374576e-03
 * time: 7.650382041931152
    60     1.786076e+02     2.860078e-03
 * time: 7.738620042800903
    61     1.786076e+02     2.100127e-03
 * time: 7.834052085876465
    62     1.786076e+02     2.100127e-03
 * time: 7.984111070632935
    63     1.786076e+02     6.412834e-03
 * time: 8.117985010147095
    64     1.786076e+02     6.412834e-03
 * time: 8.299153089523315
    65     1.786076e+02     6.412834e-03
 * time: 8.533292055130005
    66     1.786076e+02     6.412834e-03
 * time: 8.871594190597534
    67     1.786076e+02     6.412834e-03
 * time: 9.074562072753906
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: 0.00010704994201660156
     1     2.014577e+02     1.265001e+02
 * time: 0.7635269165039062
     2     1.880885e+02     4.190708e+01
 * time: 0.9961369037628174
     3     1.870317e+02     8.825666e+01
 * time: 1.1609859466552734
     4     1.846027e+02     4.400156e+01
 * time: 1.3571081161499023
     5     1.834445e+02     1.906624e+01
 * time: 1.531541109085083
     6     1.828599e+02     1.113882e+01
 * time: 1.6907451152801514
     7     1.815719e+02     7.449355e+00
 * time: 1.9133079051971436
     8     1.815131e+02     2.164678e+00
 * time: 2.1093130111694336
     9     1.814896e+02     2.167319e+00
 * time: 2.335047960281372
    10     1.814458e+02     4.615738e+00
 * time: 2.4730210304260254
    11     1.813173e+02     9.576967e+00
 * time: 2.6187009811401367
    12     1.809756e+02     2.052077e+01
 * time: 2.795315980911255
    13     1.807625e+02     4.553366e+01
 * time: 2.9571480751037598
    14     1.802224e+02     6.550892e+00
 * time: 3.1346209049224854
    15     1.800862e+02     2.865509e+00
 * time: 3.2632651329040527
    16     1.800780e+02     1.164611e+00
 * time: 3.3947110176086426
    17     1.800737e+02     7.952462e-01
 * time: 3.5190060138702393
    18     1.800352e+02     4.860618e+00
 * time: 3.662774085998535
    19     1.800089e+02     5.176689e+00
 * time: 3.8033759593963623
    20     1.799679e+02     4.303892e+00
 * time: 3.9650189876556396
    21     1.799153e+02     4.612832e+00
 * time: 4.096179008483887
    22     1.798423e+02     1.209387e+01
 * time: 4.224668979644775
    23     1.796821e+02     1.712256e+01
 * time: 4.354512929916382
    24     1.794275e+02     1.435100e+01
 * time: 4.488548994064331
    25     1.793773e+02     4.137313e+00
 * time: 4.630553960800171
    26     1.793488e+02     1.846545e+00
 * time: 4.8373329639434814
    27     1.793331e+02     5.502533e+00
 * time: 4.971108913421631
    28     1.793234e+02     2.894037e+00
 * time: 5.108158111572266
    29     1.793119e+02     1.453372e+00
 * time: 5.237026929855347
    30     1.792879e+02     4.109884e+00
 * time: 5.377150058746338
    31     1.792599e+02     4.613173e+00
 * time: 5.529743909835815
    32     1.792384e+02     4.254549e+00
 * time: 5.6778481006622314
    33     1.792251e+02     4.415024e+00
 * time: 5.839111089706421
    34     1.792026e+02     3.229145e+00
 * time: 5.9711339473724365
    35     1.791841e+02     3.422094e+00
 * time: 6.096899032592773
    36     1.791705e+02     1.621228e+00
 * time: 6.22212290763855
    37     1.791555e+02     3.462667e+00
 * time: 6.356561899185181
    38     1.791293e+02     5.586685e+00
 * time: 6.4936981201171875
    39     1.790713e+02     9.927638e+00
 * time: 6.648448944091797
    40     1.789891e+02     1.133271e+01
 * time: 6.787213087081909
    41     1.788585e+02     1.279172e+01
 * time: 6.924952030181885
    42     1.787407e+02     5.291681e+00
 * time: 7.057162046432495
    43     1.786633e+02     5.367971e+00
 * time: 7.19856595993042
    44     1.786405e+02     2.571548e+00
 * time: 7.342253923416138
    45     1.786339e+02     2.720314e+00
 * time: 7.494275093078613
    46     1.786252e+02     1.930583e+00
 * time: 7.671706914901733
    47     1.786182e+02     1.523164e+00
 * time: 7.812493085861206
    48     1.786126e+02     5.027920e-01
 * time: 7.9629480838775635
    49     1.786100e+02     3.733023e-01
 * time: 8.127568006515503
    50     1.786089e+02     2.749553e-01
 * time: 8.292948007583618
    51     1.786083e+02     1.981772e-01
 * time: 8.450971126556396
    52     1.786079e+02     1.005262e-01
 * time: 8.604620933532715
    53     1.786078e+02     2.549641e-02
 * time: 8.77908992767334
    54     1.786077e+02     4.452931e-02
 * time: 8.92313289642334
    55     1.786076e+02     2.967125e-02
 * time: 9.064459085464478
    56     1.786076e+02     1.068274e-02
 * time: 9.207327127456665
    57     1.786076e+02     5.355447e-03
 * time: 9.332042932510376
    58     1.786076e+02     8.180920e-03
 * time: 9.45642900466919
    59     1.786076e+02     4.935439e-03
 * time: 9.589504957199097
    60     1.786076e+02     4.387986e-03
 * time: 9.769643068313599
    61     1.786076e+02     4.388023e-03
 * time: 9.946478128433228
    62     1.786076e+02     4.387934e-03
 * time: 10.146687030792236
    63     1.786076e+02     4.387934e-03
 * time: 10.300440073013306
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.