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:
*
*
# 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)