Discrepancies in PRED values between Lognormal and Gaussian PK models

Simulation
1

Hello everyone,

I am comparing lognormal vs. Gaussian distributions for random effects modelling and noticed significant differences in the PRED and IPRED values between the two models when i use subject_fits() to plot.

I then tried replicating this using the non-Gaussian random effects example of the warfarin data from the tutorial (Why are non-Gaussian random effects relevant?), and observed similar discrepancies in PRED values. Nevertheless the parameter estimates are similar.

Has anyone else experienced this? Any suggestions on why these differences occur or how to resolve them?

Thanks!

Here is the code:
*

gamma_sim
gamma_sim489×367 18.3 KB

lognormal_sim
lognormal_sim479×361 39.3 KB

*


# generate population for simulation to fit model on 
Random.seed!(123)

ev1 = DosageRegimen(100, time = 0, cmt = 1, evid = 1, addl = 5, ii = 12)
pop1 = map(i -> Subject(; id = i, 
                        events = ev1), 1:25)


# lognormal 
model_lognormal = @model begin
    @param begin
        θCL ∈ RealDomain(; lower = 0, init = 0.16025)
        θVc ∈ RealDomain(; lower = 0, init = 10.262)

        ωCL ∈ RealDomain(; lower = 0, init =  0.23505)
        ωVc ∈ RealDomain(; lower = 0, init = 0.10449)

        σ ∈ RealDomain(; lower = 0, init = 0.3582)
    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

param_wo_sigma =
    (init_params(model_lognormal)
    ...,σ= 1e-2)

sim_test1 = simobs(
    model_lognormal,
    pop1,
    param_wo_sigma,
    obstimes = 0.25:1:100,
)

# fit the model to created simulation data 
fit_model_lognormal = fit(
    model_lognormal,
    Subject.(sim_test1),
    (init_params(model_lognormal)
    ...,σ= 1e-2),
    FOCE(),
    optim_options = (; iterations = 10, show_trace = true, allow_f_increases = true),
)

# plot the individual plots 
figures = subject_fits(
    inspect(fit_model_lognormal);
    separate = true,
    paginate = true,
    observations = [:dv],
    facet = (; combinelabels = true),
    axis = (;)
)

figures[1]


# fit gamma model
model_gamma = @model begin
    @param begin
        θCL ∈ RealDomain(; lower = 1e-3, init = 0.16466)
        θVc ∈ RealDomain(; lower = 1e-3, init = 10.329)

        ωCL ∈ RealDomain(; lower = 1e-3, init = 0.23348)
        ωVc ∈ RealDomain(; lower = 1e-3, init = 0.10661)

        σ ∈ RealDomain(; lower = 1e-3, init = 0.35767)
    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
        μ := @. abs(Central) / Vc
        dv ~ @. Normal(μ, μ * (σ))
    end
end

param_wo_sigma =
(init_params(model_gamma)
...,σ= 1e-2)

sim_test1 = simobs(
    model_gamma,
    pop1,
    param_wo_sigma,
    obstimes = 0.25:1:100,
)

# fit the model to created simulation data 
fit_model_gamma = fit(
    model_gamma,
    Subject.(sim_test1),
    (init_params(model_gamma)
    ...,σ= 1e-1),
    FOCE(),
    optim_options = (; iterations = 100, show_trace = true, allow_f_increases = true),
)

# plot the individual plots 
figures = subject_fits(
    inspect(fit_model_gamma);
    separate = true,
    paginate = true,
    observations = [:dv],
    facet = (; combinelabels = true),
    axis = (;)
)

figures[1]

# compare estimates 
compare_estimates(; fit_model_lognormal, fit_model_gamma)
2

This is indeed a bug in the current released version of Pumas. We recently identified and and fixed it so it will work correctly in the upcoming release of Pumas.

Meanwhile you should be able to work around the issue by computing population predictions with

simobs(
    model_lognormal,
    pop1,
    coef(fit_model_lognormal),
    fill(
        init_randeffs(
            model_lognormal,
            coef(model_lognormal)
        ),
        length(pop1)
    );
    simulate_error = false
)
2 Likes
3

Hi Andreas,

Thank you for the quick replies and good to hear that its a bug that will get fixed.

I got the suggested code working with a minor adjustment should others try the same:

simobs(
    model_lognormal,
    pop1,
    coef(fit_model_lognormal),
    fill(
        init_randeffs(
            model_lognormal,
            coef(**fit_model_lognormal**)
        ),
        length(pop1)
    );
    simulate_error = false
)
3 Likes