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 variables: Central
  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 variables: Central
  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.0201261043548584
     1     9.433464e+02     6.079483e+02
 * time: 0.3814280033111572
     2     8.189627e+02     4.423725e+02
 * time: 0.435168981552124
     3     5.917683e+02     1.819248e+02
 * time: 0.4399840831756592
     4     5.421783e+02     1.121313e+02
 * time: 0.444133996963501
     5     5.255651e+02     7.407230e+01
 * time: 0.44841814041137695
     6     5.208427e+02     8.699271e+01
 * time: 0.4524421691894531
     7     5.174883e+02     8.974584e+01
 * time: 0.4565620422363281
     8     5.138523e+02     7.328235e+01
 * time: 0.4609229564666748
     9     5.109883e+02     4.155805e+01
 * time: 0.46535205841064453
    10     5.094359e+02     3.170517e+01
 * time: 0.46944403648376465
    11     5.086172e+02     3.327331e+01
 * time: 0.5055789947509766
    12     5.080941e+02     2.942077e+01
 * time: 0.5088729858398438
    13     5.074009e+02     2.839941e+01
 * time: 0.5118250846862793
    14     5.059302e+02     3.330093e+01
 * time: 0.514415979385376
    15     5.036399e+02     3.172884e+01
 * time: 0.5172109603881836
    16     5.017004e+02     3.160020e+01
 * time: 0.5199899673461914
    17     5.008553e+02     2.599524e+01
 * time: 0.5228710174560547
    18     5.005913e+02     2.139314e+01
 * time: 0.5255730152130127
    19     5.003573e+02     2.134778e+01
 * time: 0.5286381244659424
    20     4.997249e+02     2.069868e+01
 * time: 0.531775951385498
    21     4.984453e+02     1.859010e+01
 * time: 0.535491943359375
    22     4.959584e+02     2.156209e+01
 * time: 0.539546012878418
    23     4.923347e+02     3.030833e+01
 * time: 0.5438539981842041
    24     4.906916e+02     1.652278e+01
 * time: 0.574660062789917
    25     4.902955e+02     6.360800e+00
 * time: 0.578071117401123
    26     4.902870e+02     7.028603e+00
 * time: 0.5814800262451172
    27     4.902193e+02     1.176895e+00
 * time: 0.5845470428466797
    28     4.902189e+02     1.170642e+00
 * time: 0.5869951248168945
    29     4.902186e+02     1.167624e+00
 * time: 0.5891871452331543
    30     4.902145e+02     1.110377e+00
 * time: 0.5918381214141846
    31     4.902079e+02     1.010507e+00
 * time: 0.5947401523590088
    32     4.901917e+02     9.619218e-01
 * time: 0.5976910591125488
    33     4.901683e+02     1.001006e+00
 * time: 0.6008269786834717
    34     4.901473e+02     6.138233e-01
 * time: 0.6045181751251221
    35     4.901412e+02     1.754342e-01
 * time: 0.6082191467285156
    36     4.901406e+02     2.617009e-02
 * time: 0.622840166091919
    37     4.901405e+02     4.585882e-03
 * time: 0.6254091262817383
    38     4.901405e+02     7.668184e-04
 * time: 0.6274089813232422
FittedPumasModel

Successful minimization:                      true

Likelihood approximation:                     FOCE
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: 0.00011205673217773438
     1     9.644493e+02     5.541920e+02
 * time: 0.009450912475585938
     2     8.376083e+02     3.988904e+02
 * time: 0.01663994789123535
     3     6.503174e+02     1.495777e+02
 * time: 0.022736072540283203
     4     6.121241e+02     8.826583e+01
 * time: 0.028399944305419922
     5     5.977360e+02     9.144086e+01
 * time: 0.034002065658569336
     6     5.911620e+02     1.000933e+02
 * time: 0.08651590347290039
     7     5.858372e+02     8.423844e+01
 * time: 0.09163188934326172
     8     5.821934e+02     5.194402e+01
 * time: 0.09604096412658691
     9     5.801487e+02     3.461331e+01
 * time: 0.10044693946838379
    10     5.789092e+02     3.888113e+01
 * time: 0.10505104064941406
    11     5.780054e+02     3.556605e+01
 * time: 0.1096200942993164
    12     5.769455e+02     3.624436e+01
 * time: 0.11417508125305176
    13     5.749747e+02     4.322775e+01
 * time: 0.11935901641845703
    14     5.721322e+02     3.722515e+01
 * time: 0.12450003623962402
    15     5.695879e+02     3.401586e+01
 * time: 0.1299269199371338
    16     5.683277e+02     2.854997e+01
 * time: 0.13556790351867676
    17     5.678285e+02     2.644560e+01
 * time: 0.1765739917755127
    18     5.673305e+02     2.744429e+01
 * time: 0.1812279224395752
    19     5.662430e+02     2.793918e+01
 * time: 0.18545103073120117
    20     5.641877e+02     2.616169e+01
 * time: 0.1896369457244873
    21     5.606628e+02     2.257667e+01
 * time: 0.19402694702148438
    22     5.530616e+02     3.832878e+01
 * time: 0.19926905632019043
    23     5.528349e+02     5.518159e+01
 * time: 0.20571398735046387
    24     5.497231e+02     3.042064e+01
 * time: 0.2118690013885498
    25     5.488355e+02     6.929306e+00
 * time: 0.2174530029296875
    26     5.486095e+02     1.087865e+00
 * time: 0.22309494018554688
    27     5.486062e+02     6.456402e-01
 * time: 0.2496960163116455
    28     5.486061e+02     6.467689e-01
 * time: 0.25463390350341797
    29     5.486060e+02     6.463480e-01
 * time: 0.25903892517089844
    30     5.486055e+02     6.408914e-01
 * time: 0.2634730339050293
    31     5.486045e+02     6.208208e-01
 * time: 0.26807689666748047
    32     5.486020e+02     1.035462e+00
 * time: 0.2728419303894043
    33     5.485971e+02     1.452099e+00
 * time: 0.2775731086730957
    34     5.485897e+02     1.482593e+00
 * time: 0.2820761203765869
    35     5.485839e+02     8.420646e-01
 * time: 0.2870190143585205
    36     5.485822e+02     2.023876e-01
 * time: 0.2919299602508545
    37     5.485821e+02     1.885486e-02
 * time: 0.2962961196899414
    38     5.485821e+02     2.343932e-03
 * time: 0.3114321231842041
    39     5.485821e+02     4.417566e-04
 * time: 0.31543898582458496
FittedPumasModel

Successful minimization:                      true

Likelihood approximation:                     FOCE
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
end
PumasModel
  Parameters: θkaFast, θkaSlow, θCL, θVc, θbioav, ωCL, ωVc, ξbioav, σ
  Random effects: _CL, _Vc, ηbioavLogit
  Covariates: 
  Dynamical variables: DepotFast, DepotSlow, Central
  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 variables: DepotFast, DepotSlow, Central
  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):

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: 7.605552673339844e-5
     1     2.618605e+02     1.558889e+02
 * time: 0.10933899879455566
     2     2.484481e+02     5.652093e+01
 * time: 0.1750199794769287
     3     2.466691e+02     3.819960e+01
 * time: 0.24143695831298828
     4     2.430995e+02     2.281161e+01
 * time: 0.3574190139770508
     5     2.421492e+02     1.230731e+01
 * time: 0.4085831642150879
     6     2.412891e+02     1.022578e+01
 * time: 0.46239304542541504
     7     2.398220e+02     4.018353e+00
 * time: 0.5180420875549316
     8     2.398069e+02     2.747659e+00
 * time: 0.5729420185089111
     9     2.398011e+02     1.787831e+00
 * time: 0.6267940998077393
    10     2.397895e+02     2.923097e+00
 * time: 0.7245621681213379
    11     2.397511e+02     4.766107e+00
 * time: 0.7757611274719238
    12     2.397274e+02     2.427088e+00
 * time: 0.8295199871063232
    13     2.396969e+02     2.317178e+00
 * time: 0.8846850395202637
    14     2.396432e+02     6.390949e+00
 * time: 0.9398651123046875
    15     2.395437e+02     1.431395e+01
 * time: 0.9958970546722412
    16     2.394242e+02     1.660863e+01
 * time: 1.069460153579712
    17     2.392883e+02     7.918244e+00
 * time: 1.1207640171051025
    18     2.392648e+02     5.359505e+00
 * time: 1.18402099609375
    19     2.392551e+02     8.504113e-01
 * time: 1.2382371425628662
    20     2.392536e+02     8.658360e-01
 * time: 1.2925951480865479
    21     2.392502e+02     1.359484e+00
 * time: 1.3683440685272217
    22     2.392455e+02     1.908268e+00
 * time: 1.4192450046539307
    23     2.392352e+02     3.434293e+00
 * time: 1.4721519947052002
    24     2.392094e+02     4.019050e+00
 * time: 1.5273971557617188
    25     2.391904e+02     2.360105e+00
 * time: 1.5922420024871826
    26     2.391639e+02     2.232579e+00
 * time: 1.678631067276001
    27     2.390985e+02     7.831657e+00
 * time: 1.7303359508514404
    28     2.390462e+02     1.220981e+01
 * time: 1.7819900512695312
    29     2.389849e+02     1.372376e+01
 * time: 1.8462719917297363
    30     2.389295e+02     1.012359e+01
 * time: 1.9024131298065186
    31     2.389050e+02     7.364624e+00
 * time: 1.9668641090393066
    32     2.388675e+02     3.493348e+00
 * time: 2.041321039199829
    33     2.388348e+02     3.602127e+00
 * time: 2.0922040939331055
    34     2.388092e+02     9.896348e+00
 * time: 2.1466290950775146
    35     2.387816e+02     5.815979e+00
 * time: 2.201741933822632
    36     2.387667e+02     3.995259e+00
 * time: 2.2562379837036133
    37     2.387624e+02     2.110864e+00
 * time: 2.3114640712738037
    38     2.387542e+02     3.340382e+00
 * time: 2.384658098220825
    39     2.387468e+02     3.743004e+00
 * time: 2.4347140789031982
    40     2.387148e+02     8.270607e+00
 * time: 2.4878780841827393
    41     2.386920e+02     8.043639e+00
 * time: 2.5416719913482666
    42     2.386311e+02     1.788043e+00
 * time: 2.597146987915039
    43     2.385981e+02     4.286869e+00
 * time: 2.651772975921631
    44     2.385768e+02     2.412656e+00
 * time: 2.724640130996704
    45     2.385642e+02     4.499325e-01
 * time: 2.7754950523376465
    46     2.385598e+02     1.789469e+00
 * time: 2.8282470703125
    47     2.385584e+02     4.262494e-01
 * time: 2.882390022277832
    48     2.385582e+02     3.514292e-01
 * time: 2.9356930255889893
    49     2.385581e+02     9.399545e-02
 * time: 2.989333152770996
    50     2.385580e+02     1.262250e-02
 * time: 3.062472105026245
    51     2.385580e+02     2.333105e-03
 * time: 3.1101269721984863
    52     2.385580e+02     2.333105e-03
 * time: 3.174201011657715
    53     2.385580e+02     2.333105e-03
 * time: 3.1969170570373535
FittedPumasModel

Successful minimization:                      true

Likelihood approximation:                     FOCE
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: 6.29425048828125e-5
     1     2.991160e+02     1.342892e+02
 * time: 0.1877729892730713
     2     2.851580e+02     4.511285e+01
 * time: 0.2556428909301758
     3     2.835992e+02     4.343268e+01
 * time: 0.32502102851867676
     4     2.821916e+02     4.373720e+01
 * time: 0.38392210006713867
     5     2.794524e+02     1.253857e+01
 * time: 0.44292497634887695
     6     2.787961e+02     1.102093e+01
 * time: 0.5432360172271729
     7     2.770092e+02     1.114239e+01
 * time: 0.5977859497070312
     8     2.769202e+02     5.089120e+00
 * time: 0.6552069187164307
     9     2.769088e+02     3.041932e+00
 * time: 0.7139780521392822
    10     2.768987e+02     1.603805e+00
 * time: 0.7706279754638672
    11     2.767766e+02     1.411813e+01
 * time: 0.864063024520874
    12     2.766170e+02     2.719462e+01
 * time: 0.9191699028015137
    13     2.764580e+02     3.014175e+01
 * time: 0.9756309986114502
    14     2.762515e+02     1.669154e+01
 * time: 1.034451961517334
    15     2.761434e+02     1.957266e+00
 * time: 1.0935769081115723
    16     2.761353e+02     3.756416e+00
 * time: 1.1517620086669922
    17     2.761237e+02     1.341505e+00
 * time: 1.2298040390014648
    18     2.761152e+02     2.182776e+00
 * time: 1.2836289405822754
    19     2.761036e+02     3.209868e+00
 * time: 1.3400259017944336
    20     2.760912e+02     3.070994e+00
 * time: 1.3990750312805176
    21     2.760762e+02     1.990495e+00
 * time: 1.4576199054718018
    22     2.760587e+02     1.586630e+00
 * time: 1.516191005706787
    23     2.760290e+02     2.471581e+00
 * time: 1.5946969985961914
    24     2.759737e+02     1.040810e+01
 * time: 1.6506729125976562
    25     2.758948e+02     9.671069e+00
 * time: 1.7094769477844238
    26     2.758264e+02     7.208229e+00
 * time: 1.769197940826416
    27     2.757728e+02     3.300150e+00
 * time: 1.8386609554290771
    28     2.757344e+02     4.723361e+00
 * time: 1.916698932647705
    29     2.756731e+02     6.515164e+00
 * time: 1.9717249870300293
    30     2.755937e+02     4.176876e+00
 * time: 2.0306739807128906
    31     2.755684e+02     5.161233e+00
 * time: 2.101623058319092
    32     2.755453e+02     8.103754e+00
 * time: 2.1604950428009033
    33     2.755346e+02     2.183052e+00
 * time: 2.2404870986938477
    34     2.755279e+02     3.129360e+00
 * time: 2.295146942138672
    35     2.755221e+02     3.813180e+00
 * time: 2.350477933883667
    36     2.754865e+02     6.625383e+00
 * time: 2.408895969390869
    37     2.754662e+02     5.264434e+00
 * time: 2.4669718742370605
    38     2.754486e+02     9.924709e-01
 * time: 2.5259618759155273
    39     2.754420e+02     1.102256e+00
 * time: 2.6027870178222656
    40     2.754395e+02     3.092900e-01
 * time: 2.6565239429473877
    41     2.754386e+02     5.369253e-01
 * time: 2.714401960372925
    42     2.754384e+02     1.423783e-01
 * time: 2.7705979347229004
    43     2.754383e+02     1.288451e-01
 * time: 2.82674503326416
    44     2.754381e+02     3.887239e-01
 * time: 2.90631103515625
    45     2.754376e+02     9.266434e-01
 * time: 2.959779977798462
    46     2.754363e+02     1.786193e+00
 * time: 3.013880968093872
    47     2.754331e+02     3.042628e+00
 * time: 3.071894884109497
    48     2.754256e+02     4.634171e+00
 * time: 3.1305549144744873
    49     2.754117e+02     5.954940e+00
 * time: 3.189218044281006
    50     2.753920e+02     6.732566e+00
 * time: 3.2681539058685303
    51     2.753738e+02     2.208977e+00
 * time: 3.3221750259399414
    52     2.753698e+02     2.481125e+00
 * time: 3.3791260719299316
    53     2.753538e+02     9.838345e-01
 * time: 3.4377779960632324
    54     2.753535e+02     3.815520e+00
 * time: 3.4955050945281982
    55     2.753483e+02     1.735894e+00
 * time: 3.554702043533325
    56     2.753436e+02     8.278884e-01
 * time: 3.6424169540405273
    57     2.753339e+02     1.256117e+00
 * time: 3.6961560249328613
    58     2.753262e+02     2.343600e+00
 * time: 3.7537600994110107
    59     2.753231e+02     1.270829e+00
 * time: 3.810537099838257
    60     2.753207e+02     3.103092e-01
 * time: 3.867335081100464
    61     2.753202e+02     3.488786e-01
 * time: 3.9231810569763184
    62     2.753193e+02     4.518882e-01
 * time: 3.9989800453186035
    63     2.753186e+02     2.027309e-01
 * time: 4.051621913909912
    64     2.753183e+02     9.363635e-02
 * time: 4.108617067337036
    65     2.753182e+02     9.589112e-02
 * time: 4.163969039916992
    66     2.753181e+02     9.484170e-02
 * time: 4.2190680503845215
    67     2.753181e+02     4.831472e-02
 * time: 4.297069072723389
    68     2.753180e+02     2.328249e-02
 * time: 4.349025011062622
    69     2.753180e+02     1.740503e-02
 * time: 4.399848937988281
    70     2.753180e+02     1.793621e-02
 * time: 4.4543139934539795
    71     2.753180e+02     8.549994e-03
 * time: 4.508888006210327
    72     2.753180e+02     5.684889e-03
 * time: 4.563337087631226
    73     2.753180e+02     4.082116e-03
 * time: 4.622344970703125
    74     2.753180e+02     7.506601e-03
 * time: 4.696007966995239
    75     2.753180e+02     6.945870e-03
 * time: 4.746423959732056
    76     2.753180e+02     5.205184e-03
 * time: 4.804795980453491
    77     2.753180e+02     2.313568e-03
 * time: 4.858575105667114
    78     2.753180e+02     2.313568e-03
 * time: 4.911381006240845
FittedPumasModel

Successful minimization:                      true

Likelihood approximation:                     FOCE
Log-likelihood value:                   -275.31801
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.13275
ωCL          0.081471
ωVc          0.10178
nbioav       1.3974e8
σ            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.9882
2 θkaSlow 0.607708 0.613328
3 θCL 1.07968 1.07967
4 θVc 11.3745 11.381
5 θbioav 0.143942 0.132746
6 ωCL 0.0815101 0.0814711
7 ωVc 0.101565 0.101784
8 σ 0.10441 0.104644
9 ξbioav 0.125128 missing
10 nbioav missing 1.39739e8

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.54e-269 0.0
4 0.003003 4.48998e-222 0.0
5 0.004004 3.80318e-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.