Variational Expectation Maximization (VEM)

Authors

Mohamed Tarek

Vijay Ivaturi

WarningExperimental Feature

VEM is an experimental estimation method introduced in Pumas 2.8. Its API may change in future releases without deprecation.

In this tutorial, we introduce Variational Expectation Maximization (VEM), a new estimation method available in Pumas.Experimental. We will cover the conceptual background, fit a warfarin PK model with VEM, compare the results to FOCE, and provide a comprehensive reference of VEM options.

using Pumas
using Pumas.Experimental: VEM
using PumasUtilities
using PharmaDatasets
using DataFramesMeta
using CairoMakie
using AlgebraOfGraphics
using SummaryTables
using Random: default_rng, seed!
seed!(123)
Random.TaskLocalRNG()

1 Conceptual Background

1.1 The Marginal Likelihood Challenge

In nonlinear mixed-effects (NLME) modeling, we want to estimate population parameters while accounting for between-subject variability through random effects. The marginal likelihood requires integrating over all possible random effect values for every subject, and these integrals are generally intractable.

Different estimation methods handle this challenge in different ways:

  • FOCE/LaplaceI approximate the integral using a second-order Taylor expansion around the mode of the random effects.
  • SAEM uses Markov Chain Monte Carlo (MCMC) sampling within an Expectation-Maximization framework.
  • VEM maximizes a lower bound on the log-likelihood using variational inference.

1.2 The Evidence Lower Bound (ELBO)

VEM works by maximizing the Evidence Lower Bound (ELBO), which is a lower bound on the log marginal likelihood. The ELBO equals the log-likelihood minus a KL divergence term that measures how far the approximate posterior is from the true posterior. A higher ELBO typically means a better fit.

Intuitively, VEM simultaneously finds:

  1. Population parameters (fixed effects, variance components, residual error), and
  2. An approximation to each subject’s random effect posterior distribution.

This joint optimization is performed in a single level, unlike FOCE/LaplaceI which alternates between inner (random effects) and outer (population parameters) optimization loops.

One key advantage of VEM is that it can handle any response data distribution, including discrete, censored, truncated, time-to-event and even custom distributions. FOCE by comparison is limited to a restricted set of distributions, and while LaplaceI can handle any response distribution, it can be very slow for complex models.

1.3 Presampling vs. Resampling

VEM has two optimization variants controlled by the presample keyword:

Presampling (presample=true, default) Resampling (presample=false)
Optimizer L-BFGS (deterministic) Adam (stochastic)
Phases 1 (default) 10 (default)
Samples Pre-drawn at start of phase, reused across iterations Re-drawn at every iteration

In the resampling variant, the number of samples increases (by default, unless overridden) and the learning rate of Adam decreases from one phase to the next.

1.4 Diagonal vs. Full Posterior Covariance

VEM approximates each subject’s random effect posterior with a multivariate normal distribution. The diagonal keyword controls the structure of this approximation:

  • diagonal=true (default): assumes independent random effects in the posterior. Faster, fewer parameters to optimize.
  • diagonal=false: captures correlations between random effects in the posterior. Slower but may be more accurate when random effects are correlated.

2 Warfarin PKPD Model with VEM

We will fit an oral PKPD model to the warfarin dataset, first using VEM and then comparing with FOCE.

2.1 Model Definition

warfarin_model = @model begin
    @param begin
        # PK
        pop_CL ∈ RealDomain(lower = 0.0, init = 0.134)
        pop_V ∈ RealDomain(lower = 0.0, init = 8.11)
        pop_tabs ∈ RealDomain(lower = 0.0, init = 0.523)
        pop_lag ∈ RealDomain(lower = 0.0, init = 0.1)
        pk_Ω ∈ PSDDomain([0.01, 0.01, 0.01])
        σ_prop ∈ RealDomain(lower = 0.0, init = 0.00752)
        σ_add ∈ RealDomain(lower = 0.0, init = 0.0661)
        # PD
        pop_e0 ∈ RealDomain(lower = 0.0, init = 100.0)
        pop_emax ∈ RealDomain(init = -1.0)
        pop_c50 ∈ RealDomain(lower = 0.0, init = 1.0)
        pop_tover ∈ RealDomain(lower = 0.0, init = 14.0)
        pd_Ω ∈ PSDDomain([0.01, 0.01, 0.01])
        σ_fx ∈ RealDomain(lower = 0.0, init = 0.01)
    end

    @random begin
        pk_η ~ MvNormal(pk_Ω)
        pd_η ~ MvNormal(pd_Ω)
    end

    @covariates begin
        FSZV
        FSZCL
    end

    @pre begin
        # PK
        CL = FSZCL * pop_CL * exp(pk_η[1])
        Vc = FSZV * pop_V * exp(pk_η[2])
        tabs = pop_tabs * exp(pk_η[3])
        Ka = log(2) / tabs
        # PD
        e0 = pop_e0 * exp(pd_η[1])
        emax = pop_emax
        c50 = pop_c50 * exp(pd_η[2])
        tover = pop_tover * exp(pd_η[3])
        kout = log(2) / tover
        rin = e0 * kout
    end

    @dosecontrol begin
        lags = (Depot = pop_lag,)
    end

    @init begin
        Turnover = e0
    end

    @vars begin
        cp := Central / Vc
        ratein := Ka * Depot
        pd := 1 + emax * cp / (c50 + cp)
    end

    @dynamics begin
        Depot' = -ratein
        Central' = ratein - CL * cp
        Turnover' = rin * pd - kout * Turnover
    end

    @derived begin
        conc ~ @. CombinedNormal(cp, σ_add, σ_prop)
        pca ~ @. Normal(Turnover, σ_fx)
    end
end
PumasModel
  Parameters: pop_CL, pop_V, pop_tabs, pop_lag, pk_Ω, σ_prop, σ_add, pop_e0, pop_emax, pop_c50, pop_tover, pd_Ω, σ_fx
  Random effects: pk_η, pd_η
  Covariates: FSZV, FSZCL
  Dynamical system variables: Depot, Central, Turnover
  Dynamical system type: Nonlinear ODE
  Derived: conc, pca
  Observed: conc, pca
TipSame Model for VEM and FOCE

The @model definition is identical regardless of whether you fit with VEM or FOCE/LaplaceI. VEM works with standard @model definitions – no special model macro is required.

NoteRandom effect on emax

Note that no random effect is used on emax. We tried adding a random effect on emax and estimating the model parameters. However, that resulted in the row and column of pd_Ω corresponding to that random effect to be all close to 0 which caused numerical issues for the asymptotic esimtator of the confidence interval. So no random effect is used on emax here for simplicity.

2.2 Data Loading

warfarin_data = @chain dataset("pumas/warfarin_pumas") begin
    @rtransform begin
        :FSZV = :wtbl / 70
        :FSZCL = (:wtbl / 70)^0.75
    end
end
first(warfarin_data, 5)
5×12 DataFrame
Row id time evid amt cmt conc pca wtbl age sex FSZV FSZCL
Int64 Float64 Int64 Float64? Int64? Float64? Float64? Float64 Int64 String1 Float64 Float64
1 1 0.0 1 100.0 1 missing missing 66.7 50 M 0.952857 0.96443
2 1 0.5 0 missing missing 0.0 missing 66.7 50 M 0.952857 0.96443
3 1 1.0 0 missing missing 1.9 missing 66.7 50 M 0.952857 0.96443
4 1 2.0 0 missing missing 3.3 missing 66.7 50 M 0.952857 0.96443
5 1 3.0 0 missing missing 6.6 missing 66.7 50 M 0.952857 0.96443 
pop = read_pumas(
    warfarin_data;
    id = :id,
    time = :time,
    amt = :amt,
    cmt = :cmt,
    evid = :evid,
    covariates = [:FSZV, :FSZCL],
    observations = [:conc, :pca],
)
Population
  Subjects: 32
  Covariates: FSZV, FSZCL
  Observations: conc, pca

2.3 Fitting with Default VEM

Fitting with VEM follows the same fit interface as all other estimation methods. Simply pass VEM() as the estimation algorithm:

@time fpm_vem = fit(
    warfarin_model,
    pop,
    init_params(warfarin_model),
    VEM(; optim_options = (; show_every = 20)),
)
[ 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     1.395960e+08     3.459221e+08
 * time: 2.6941299438476562e-5
    20     1.782428e+03     1.102616e+03
 * time: 10.418637037277222
    40     1.158299e+03     1.183852e+02
 * time: 17.910160064697266
    60     1.082972e+03     5.027322e+02
 * time: 25.838254928588867
    80     1.033477e+03     1.661968e+02
 * time: 34.972655057907104
   100     9.683433e+02     9.345296e+01
 * time: 45.8673620223999
   120     9.206464e+02     2.761997e+01
 * time: 57.22170615196228
   140     8.862514e+02     9.321403e+01
 * time: 70.1200110912323
   160     8.339402e+02     2.689859e+01
 * time: 83.05230498313904
   180     7.830751e+02     5.439420e+00
 * time: 97.00664114952087
   200     7.463957e+02     4.798360e+01
 * time: 111.47983503341675
   220     7.338703e+02     2.505895e+01
 * time: 126.80636501312256
   240     7.269872e+02     6.267607e+00
 * time: 142.92881512641907
   260     7.245847e+02     3.532142e+00
 * time: 160.94051504135132
   280     7.237545e+02     6.477516e+00
 * time: 177.29596710205078
195.636728 seconds (285.32 M allocations: 57.369 GiB, 6.37% gc time, 0.96% compilation time)

FittedPumasModel

Dynamical system type:               Nonlinear ODE
Solver(s): (OrdinaryDiffEqVerner.Vern7,OrdinaryDiffEqRosenbrock.Rodas5P)

Number of subjects:                             32

Observation records:         Active        Missing
    conc:                       251             47
    pca:                        232             66
    Total:                      483            113

Number of parameters:      Constant      Optimized
                                  0             23

Estimation method:              VEM (EXPERIMENTAL)

------------------------
            Estimate
------------------------
pop_CL       0.1344
pop_V        8.0181
pop_tabs     0.40771
pop_lag      0.88161
pk_Ω₁,₁      0.065883
pk_Ω₂,₁      0.0089982
pk_Ω₃,₁      0.034655
pk_Ω₂,₂      0.020545
pk_Ω₃,₂     -0.026903
pk_Ω₃,₃      0.92808
σ_prop       0.088454
σ_add        0.38888
pop_e0      96.476
pop_emax    -1.0714
pop_c50      1.5354
pop_tover   14.161
pd_Ω₁,₁      0.00031903
pd_Ω₂,₁      0.006396
pd_Ω₃,₁     -0.0013854
pd_Ω₂,₂      0.1291
pd_Ω₃,₂     -0.030921
pd_Ω₃,₃      0.017369
σ_fx         4.332
------------------------

During fitting, VEM displays a convergence trace showing the negative ELBO (the objective being minimized). The objective value is the negative of the lower bound on the log marginal likelihood, so lower values indicate a better fit.

2.4 Extracting Results

Population parameter estimates can be extracted with coef:

simple_table(coeftable(fpm_vem))
parameter constant estimate
pop_CL false 0.134
pop_V false 8.02
pop_tabs false 0.408
pop_lag false 0.882
pk_Ω₁,₁ false 0.0659
pk_Ω₂,₁ false 0.009
pk_Ω₃,₁ false 0.0347
pk_Ω₂,₂ false 0.0205
pk_Ω₃,₂ false -0.0269
pk_Ω₃,₃ false 0.928
σ_prop false 0.0885
σ_add false 0.389
pop_e0 false 96.5
pop_emax false -1.07
pop_c50 false 1.54
pop_tover false 14.2
pd_Ω₁,₁ false 0.000319
pd_Ω₂,₁ false 0.0064
pd_Ω₃,₁ false -0.00139
pd_Ω₂,₂ false 0.129
pd_Ω₃,₂ false -0.0309
pd_Ω₃,₃ false 0.0174
σ_fx false 4.33

The posterior approximation modes of the random effects for each subject are available via empirical_bayes:

ebes = empirical_bayes(fpm_vem)
first(ebes, 3)
3-element Vector{@NamedTuple{pk_η::Vector{Float64}, pd_η::Vector{Float64}}}:
 (pk_η = [0.7702024660066762, 0.06125565173800886, 1.522628941795522], pd_η = [-0.006178465748795943, -0.13407626817045565, 0.06422442650436823])
 (pk_η = [0.02537465589436046, 0.05748436584165914, 0.4135965686975201], pd_η = [0.020340484329945556, 0.4114666774513079, -0.10179111154486246])
 (pk_η = [-0.29756294829753516, -0.0372087147891645, -0.3448171786273172], pd_η = [-0.027228113648185465, -0.5470013133386835, 0.11987520069672504])

Note that these generally may not correspond to modes of the true posterior.

2.5 Log-Likelihood Computation

VEM does not compute the log-likelihood during fitting. To evaluate the log-likelihood after fitting, you must specify an approximation method:

ll_foce = loglikelihood(fpm_vem, FOCE())
-1072.1979779880155
WarningLog-Likelihood Approximation Starting Points

When using FOCE() (or LaplaceI()) for log-likelihood evaluation, the empirical Bayes estimates (EBEs) are computed for each subject starting from the mode of the prior (i.e., zero random effects in the Gaussian case) rather than the VEM posterior approximation’s mode.

For well-behaved models, this is unlikely to cause a significant difference in the log-likelihood approximation. However, this behavior will likely change in a future release to use the VEM posterior approximation’s mode as the starting point for the Laplace approximation.

2.6 Diagnostics

Model diagnostics are available through inspect. For models with Gaussian error models (like our proportional normal error), you must explicitly specify how to compute weighted residuals using wres_approx:

res = inspect(fpm_vem; wres_approx = FOCE())
[ Info: Calculating predictions.
[ Info: Calculating weighted residuals.
[ Info: Calculating empirical bayes.
[ Info: Evaluating dose control parameters.
[ Info: Evaluating individual parameters.
[ Info: Done.

FittedPumasModelInspection

Likelihood approximation used for weighted residuals: FOCE
NoteWeighted Residuals Linearization Point

When computing the weights of the weighted residuals, the linearization of the mean and variance are performed at the mode of the approximate posterior for each subject, learned by VEM, rather than the mode of the true posterior.

The predictions used for the residuals are also computed using the mode of the VEM posterior approximation, which may differ from the EBEs or the modes of the true posterior.

Goodness-of-fit plots work as other estimation methods in Pumas:

goodness_of_fit(res, observations = [:conc])

goodness_of_fit(res, observations = [:pca])

Warningwres_approx Required for Gaussian Errors

When using VEM with Gaussian-like error models (Normal, ProportionalNormal, etc.), you must pass wres_approx = FOCE() or wres_approx = LaplaceI() to inspect. Omitting it will raise an error. For other error models (Beta, Gamma, Poisson, Bernoulli, TimeToEvent), this is not needed.

2.7 Predictions and VPCs

Predictions and visual predictive checks work as other estimation methods in Pumas:

preds = predict(fpm_vem)
first(preds, 3)
Prediction
  Subjects: 3
  Predictions: conc, pca
  Covariates: FSZV, FSZCL
NoteVEM Posterior Approximation’s Mode in Predictions

The individual predictions are computed using the mode of the VEM posterior approximation for each subject, which may differ from the EBEs or the modes of the true posterior.

conc_vpc_res = vpc(fpm_vem, observations = [:conc]);
vpc_plot(conc_vpc_res)

pca_vpc_res = vpc(fpm_vem, observations = [:pca]);
vpc_plot(pca_vpc_res)

2.8 Inference Workaround

VEM does not currently support infer for computing standard errors and confidence intervals:

infer(fpm_vem)
ArgumentError: Asymptotic inference is not supported for the estimation method VEM (EXPERIMENTAL).
Stacktrace:
 [1] infer(fpm::Pumas.FittedPumasModel{PumasModel{(pop_CL = 1, pop_V = 1, pop_tabs = 1, pop_lag = 1, pk_Ω = 6, σ_prop = 1, σ_add = 1, pop_e0 = 1, pop_emax = 1, pop_c50 = 1, pop_tover = 1, pd_Ω = 6, σ_fx = 1), 6, (:Depot, :Central, :Turnover), (:CL, :Vc, :tabs, :Ka, :e0, :emax, :c50, :tover, :kout, :rin), ParamSet{@NamedTuple{pop_CL::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_V::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tabs::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_lag::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pk_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_prop::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, σ_add::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_e0::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_emax::RealDomain{TransformVariables.Infinity{false}, TransformVariables.Infinity{true}, Float64}, pop_c50::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tover::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pd_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_fx::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}, ParamSet{@NamedTuple{pk_η::Pumas.RandomEffect{false, var"#1#12"}, pd_η::Pumas.RandomEffect{false, var"#2#13"}}}, Pumas.PreObj{Pumas.TimeDispatcher{var"#3#14", var"#4#15"}}, Pumas.DCPObj{var"#6#17"}, var"#7#18", SciMLBase.ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, Nothing, SciMLBase.ODEFunction{false, SciMLBase.FullSpecialize, ModelingToolkit.GeneratedFunctionWrapper{(2, 3, false), var"#8#19", var"#9#20"}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ODESystem, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardODEProblem}, Pumas.DerivedObj{(:conc, :pca), false, var"#10#21"}, var"#11#22", ModelingToolkit.ODESystem, Pumas.PumasModelOptions}, Vector{Subject{@NamedTuple{conc::Vector{Union{Missing, Float64}}, pca::Vector{Union{Missing, Float64}}}, Pumas.ConstantCovar{@NamedTuple{FSZV::Float64, FSZCL::Float64}}, DosageRegimen{Float64, Union{Missing, Int64}, Float64, Float64, Float64, Float64}, Vector{Float64}}}, Pumas.Experimental.VEMResult{Pumas.Experimental.ThreadedVEMObj{true, false, true, true, PumasModel{(pop_CL = 1, pop_V = 1, pop_tabs = 1, pop_lag = 1, pk_Ω = 6, σ_prop = 1, σ_add = 1, pop_e0 = 1, pop_emax = 1, pop_c50 = 1, pop_tover = 1, pd_Ω = 6, σ_fx = 1), 6, (:Depot, :Central, :Turnover), (:CL, :Vc, :tabs, :Ka, :e0, :emax, :c50, :tover, :kout, :rin), ParamSet{@NamedTuple{pop_CL::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_V::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tabs::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_lag::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pk_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_prop::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, σ_add::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_e0::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_emax::RealDomain{TransformVariables.Infinity{false}, TransformVariables.Infinity{true}, Float64}, pop_c50::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tover::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pd_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_fx::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}, ParamSet{@NamedTuple{pk_η::Pumas.RandomEffect{false, var"#1#12"}, pd_η::Pumas.RandomEffect{false, var"#2#13"}}}, Pumas.PreObj{Pumas.TimeDispatcher{var"#3#14", var"#4#15"}}, Pumas.DCPObj{var"#6#17"}, var"#7#18", SciMLBase.ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, Nothing, SciMLBase.ODEFunction{false, SciMLBase.FullSpecialize, ModelingToolkit.GeneratedFunctionWrapper{(2, 3, false), var"#8#19", var"#9#20"}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ODESystem, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardODEProblem}, Pumas.DerivedObj{(:conc, :pca), false, var"#10#21"}, var"#11#22", ModelingToolkit.ODESystem, Pumas.PumasModelOptions}, Vector{Subject{@NamedTuple{conc::Vector{Union{Missing, Float64}}, pca::Vector{Union{Missing, Float64}}}, Pumas.ConstantCovar{@NamedTuple{FSZV::Float64, FSZCL::Float64}}, DosageRegimen{Float64, Union{Missing, Int64}, Float64, Float64, Float64, Float64}, Vector{Float64}}}, @NamedTuple{pop_CL::Float64, pop_V::Float64, pop_tabs::Float64, pop_lag::Float64, pk_Ω::PDMats.PDiagMat{Float64, Vector{Float64}}, σ_prop::Float64, σ_add::Float64, pop_e0::Float64, pop_emax::Float64, pop_c50::Float64, pop_tover::Float64, pd_Ω::PDMats.PDiagMat{Float64, Vector{Float64}}, σ_fx::Float64}, Random.TaskLocalRNG, Nothing, Nothing, @NamedTuple{ss_abstol::Float64, ss_reltol::Float64, ss_maxiters::Int64, abstol::Float64, reltol::Float64, alg::CompositeAlgorithm{0, Tuple{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}}, OrdinaryDiffEqCore.AutoSwitch{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}, Rational{Int64}, Int64}}, continuity::Symbol}, Vector{Array{Float64, 3}}, TransformVariables.TransformTuple{@NamedTuple{pop_CL::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_V::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tabs::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_lag::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pk_Ω::Pumas.PSDTransform, σ_prop::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, σ_add::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_e0::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_emax::TransformVariables.Identity, pop_c50::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tover::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pd_Ω::Pumas.PSDTransform, σ_fx::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}}}, TransformVariables.TransformTuple{@NamedTuple{pop_CL::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_V::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tabs::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_lag::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pk_Ω::TransformVariables.Constant{Nothing}, σ_prop::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, σ_add::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_e0::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_emax::TransformVariables.Identity, pop_c50::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tover::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pd_Ω::TransformVariables.Constant{Nothing}, σ_fx::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}}}, Vector{Vector{Float64}}, Vector{Vector{Float64}}, DifferentiationInterfaceForwardDiffExt.ForwardDiffGradientPrep{Tuple{Pumas._L_param{false, (:pop_CL, :pop_V, :pop_tabs, :pop_lag, :pk_Ω, :σ_prop, :σ_add, :pop_e0, :pop_emax, :pop_c50, :pop_tover, :pd_Ω, :σ_fx), TransformVariables.TransformTuple{@NamedTuple{pop_CL::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_V::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tabs::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_lag::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pk_Ω::Pumas.PSDTransform, σ_prop::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, σ_add::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_e0::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_emax::TransformVariables.Identity, pop_c50::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tover::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pd_Ω::Pumas.PSDTransform, σ_fx::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}}}, ParamSet{@NamedTuple{pop_CL::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_V::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tabs::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_lag::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pk_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_prop::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, σ_add::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_e0::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_emax::RealDomain{TransformVariables.Infinity{false}, TransformVariables.Infinity{true}, Float64}, pop_c50::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tover::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pd_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_fx::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}}, ADTypes.AutoForwardDiff{nothing, Nothing}, SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}, Tuple{}}, ForwardDiff.GradientConfig{ForwardDiff.Tag{Pumas._L_param{false, (:pop_CL, :pop_V, :pop_tabs, :pop_lag, :pk_Ω, :σ_prop, :σ_add, :pop_e0, :pop_emax, :pop_c50, :pop_tover, :pd_Ω, :σ_fx), TransformVariables.TransformTuple{@NamedTuple{pop_CL::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_V::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tabs::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_lag::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pk_Ω::Pumas.PSDTransform, σ_prop::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, σ_add::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_e0::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_emax::TransformVariables.Identity, pop_c50::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tover::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pd_Ω::Pumas.PSDTransform, σ_fx::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}}}, ParamSet{@NamedTuple{pop_CL::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_V::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tabs::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_lag::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pk_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_prop::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, σ_add::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_e0::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_emax::RealDomain{TransformVariables.Infinity{false}, TransformVariables.Infinity{true}, Float64}, pop_c50::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tover::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pd_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_fx::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}}, Float64}, Float64, 12, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Pumas._L_param{false, (:pop_CL, :pop_V, :pop_tabs, :pop_lag, :pk_Ω, :σ_prop, :σ_add, :pop_e0, :pop_emax, :pop_c50, :pop_tover, :pd_Ω, :σ_fx), TransformVariables.TransformTuple{@NamedTuple{pop_CL::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_V::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tabs::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_lag::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pk_Ω::Pumas.PSDTransform, σ_prop::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, σ_add::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_e0::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_emax::TransformVariables.Identity, pop_c50::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tover::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pd_Ω::Pumas.PSDTransform, σ_fx::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}}}, ParamSet{@NamedTuple{pop_CL::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_V::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tabs::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_lag::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pk_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_prop::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, σ_add::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_e0::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_emax::RealDomain{TransformVariables.Infinity{false}, TransformVariables.Infinity{true}, Float64}, pop_c50::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tover::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pd_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_fx::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}}, Float64}, Float64, 12}}}, Tuple{}}, Vector{DifferentiationInterfaceReverseDiffExt.ReverseDiffGradientPrep{Tuple{typeof(Pumas.Experimental.G1_prior), ADTypes.AutoReverseDiff{false}, Vector{Float64}, Tuple{DifferentiationInterface.Constant{Pumas.Experimental._L_rfx2{true, (:pop_CL, :pop_V, :pop_tabs, :pop_lag, :pk_Ω, :σ_prop, :σ_add, :pop_e0, :pop_emax, :pop_c50, :pop_tover, :pd_Ω, :σ_fx), TransformVariables.TransformTuple{@NamedTuple{pop_CL::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_V::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tabs::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_lag::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pk_Ω::Pumas.PSDTransform, σ_prop::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, σ_add::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_e0::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_emax::TransformVariables.Identity, pop_c50::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tover::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pd_Ω::Pumas.PSDTransform, σ_fx::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}}}, TransformVariables.TransformTuple{@NamedTuple{pop_CL::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_V::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tabs::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_lag::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pk_Ω::TransformVariables.Constant{Nothing}, σ_prop::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, σ_add::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_e0::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_emax::TransformVariables.Identity, pop_c50::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tover::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pd_Ω::TransformVariables.Constant{Nothing}, σ_fx::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}}}, PumasModel{(pop_CL = 1, pop_V = 1, pop_tabs = 1, pop_lag = 1, pk_Ω = 6, σ_prop = 1, σ_add = 1, pop_e0 = 1, pop_emax = 1, pop_c50 = 1, pop_tover = 1, pd_Ω = 6, σ_fx = 1), 6, (:Depot, :Central, :Turnover), (:CL, :Vc, :tabs, :Ka, :e0, :emax, :c50, :tover, :kout, :rin), ParamSet{@NamedTuple{pop_CL::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_V::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tabs::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_lag::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pk_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_prop::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, σ_add::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_e0::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_emax::RealDomain{TransformVariables.Infinity{false}, TransformVariables.Infinity{true}, Float64}, pop_c50::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tover::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pd_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_fx::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}, ParamSet{@NamedTuple{pk_η::Pumas.RandomEffect{false, var"#1#12"}, pd_η::Pumas.RandomEffect{false, var"#2#13"}}}, Pumas.PreObj{Pumas.TimeDispatcher{var"#3#14", var"#4#15"}}, Pumas.DCPObj{var"#6#17"}, var"#7#18", SciMLBase.ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, Nothing, SciMLBase.ODEFunction{false, SciMLBase.FullSpecialize, ModelingToolkit.GeneratedFunctionWrapper{(2, 3, false), var"#8#19", var"#9#20"}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ODESystem, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardODEProblem}, Pumas.DerivedObj{(:conc, :pca), false, var"#10#21"}, var"#11#22", ModelingToolkit.ODESystem, Pumas.PumasModelOptions}, Subject{@NamedTuple{conc::Vector{Union{Missing, Float64}}, pca::Vector{Union{Missing, Float64}}}, Pumas.ConstantCovar{@NamedTuple{FSZV::Float64, FSZCL::Float64}}, DosageRegimen{Float64, Union{Missing, Int64}, Float64, Float64, Float64, Float64}, Vector{Float64}}, @NamedTuple{ss_abstol::Float64, ss_reltol::Float64, ss_maxiters::Int64, abstol::Float64, reltol::Float64, alg::CompositeAlgorithm{0, Tuple{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}}, OrdinaryDiffEqCore.AutoSwitch{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}, Rational{Int64}, Int64}}, continuity::Symbol}, Matrix{Float64}, @NamedTuple{pop_CL::Float64, pop_V::Float64, pop_tabs::Float64, pop_lag::Float64, pk_Ω::PDMats.PDMat{Float64, Matrix{Float64}}, σ_prop::Float64, σ_add::Float64, pop_e0::Float64, pop_emax::Float64, pop_c50::Float64, pop_tover::Float64, pd_Ω::PDMats.PDMat{Float64, Matrix{Float64}}, σ_fx::Float64}}}}}, ReverseDiff.GradientConfig{ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}}, Nothing}}, Vector{DifferentiationInterfaceForwardDiffExt.ForwardDiffGradientPrep{Tuple{typeof(Pumas.Experimental.G1_likelihood), ADTypes.AutoForwardDiff{nothing, Nothing}, Vector{Float64}, Tuple{DifferentiationInterface.Constant{Pumas.Experimental._L_rfx2{true, (:pop_CL, :pop_V, :pop_tabs, :pop_lag, :pk_Ω, :σ_prop, :σ_add, :pop_e0, :pop_emax, :pop_c50, :pop_tover, :pd_Ω, :σ_fx), TransformVariables.TransformTuple{@NamedTuple{pop_CL::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_V::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tabs::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_lag::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pk_Ω::Pumas.PSDTransform, σ_prop::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, σ_add::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_e0::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_emax::TransformVariables.Identity, pop_c50::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tover::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pd_Ω::Pumas.PSDTransform, σ_fx::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}}}, TransformVariables.TransformTuple{@NamedTuple{pop_CL::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_V::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tabs::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_lag::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pk_Ω::TransformVariables.Constant{Nothing}, σ_prop::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, σ_add::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_e0::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_emax::TransformVariables.Identity, pop_c50::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tover::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pd_Ω::TransformVariables.Constant{Nothing}, σ_fx::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}}}, PumasModel{(pop_CL = 1, pop_V = 1, pop_tabs = 1, pop_lag = 1, pk_Ω = 6, σ_prop = 1, σ_add = 1, pop_e0 = 1, pop_emax = 1, pop_c50 = 1, pop_tover = 1, pd_Ω = 6, σ_fx = 1), 6, (:Depot, :Central, :Turnover), (:CL, :Vc, :tabs, :Ka, :e0, :emax, :c50, :tover, :kout, :rin), ParamSet{@NamedTuple{pop_CL::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_V::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tabs::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_lag::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pk_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_prop::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, σ_add::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_e0::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_emax::RealDomain{TransformVariables.Infinity{false}, TransformVariables.Infinity{true}, Float64}, pop_c50::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tover::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pd_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_fx::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}, ParamSet{@NamedTuple{pk_η::Pumas.RandomEffect{false, var"#1#12"}, pd_η::Pumas.RandomEffect{false, var"#2#13"}}}, Pumas.PreObj{Pumas.TimeDispatcher{var"#3#14", var"#4#15"}}, Pumas.DCPObj{var"#6#17"}, var"#7#18", SciMLBase.ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, Nothing, SciMLBase.ODEFunction{false, SciMLBase.FullSpecialize, ModelingToolkit.GeneratedFunctionWrapper{(2, 3, false), var"#8#19", var"#9#20"}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ODESystem, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardODEProblem}, Pumas.DerivedObj{(:conc, :pca), false, var"#10#21"}, var"#11#22", ModelingToolkit.ODESystem, Pumas.PumasModelOptions}, Subject{@NamedTuple{conc::Vector{Union{Missing, Float64}}, pca::Vector{Union{Missing, Float64}}}, Pumas.ConstantCovar{@NamedTuple{FSZV::Float64, FSZCL::Float64}}, DosageRegimen{Float64, Union{Missing, Int64}, Float64, Float64, Float64, Float64}, Vector{Float64}}, @NamedTuple{ss_abstol::Float64, ss_reltol::Float64, ss_maxiters::Int64, abstol::Float64, reltol::Float64, alg::CompositeAlgorithm{0, Tuple{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}}, OrdinaryDiffEqCore.AutoSwitch{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}, Rational{Int64}, Int64}}, continuity::Symbol}, Matrix{Float64}, @NamedTuple{pop_CL::Float64, pop_V::Float64, pop_tabs::Float64, pop_lag::Float64, pk_Ω::PDMats.PDMat{Float64, Matrix{Float64}}, σ_prop::Float64, σ_add::Float64, pop_e0::Float64, pop_emax::Float64, pop_c50::Float64, pop_tover::Float64, pd_Ω::PDMats.PDMat{Float64, Matrix{Float64}}, σ_fx::Float64}}}}}, ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(Pumas.Experimental.G1_likelihood), Float64}, Float64, 9, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(Pumas.Experimental.G1_likelihood), Float64}, Float64, 9}}}, Tuple{Nothing}}}, Vector{DifferentiationInterfaceReverseDiffExt.ReverseDiffTwoArgJacobianPrep{Tuple{typeof(Pumas.Experimental.get_randeffs!), SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}, ADTypes.AutoReverseDiff{false}, Vector{Float64}, Tuple{DifferentiationInterface.Constant{Pumas.Experimental._L_rfx2{true, (:pop_CL, :pop_V, :pop_tabs, :pop_lag, :pk_Ω, :σ_prop, :σ_add, :pop_e0, :pop_emax, :pop_c50, :pop_tover, :pd_Ω, :σ_fx), TransformVariables.TransformTuple{@NamedTuple{pop_CL::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_V::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tabs::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_lag::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pk_Ω::Pumas.PSDTransform, σ_prop::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, σ_add::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_e0::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_emax::TransformVariables.Identity, pop_c50::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tover::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pd_Ω::Pumas.PSDTransform, σ_fx::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}}}, TransformVariables.TransformTuple{@NamedTuple{pop_CL::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_V::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tabs::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_lag::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pk_Ω::TransformVariables.Constant{Nothing}, σ_prop::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, σ_add::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_e0::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_emax::TransformVariables.Identity, pop_c50::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tover::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pd_Ω::TransformVariables.Constant{Nothing}, σ_fx::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}}}, PumasModel{(pop_CL = 1, pop_V = 1, pop_tabs = 1, pop_lag = 1, pk_Ω = 6, σ_prop = 1, σ_add = 1, pop_e0 = 1, pop_emax = 1, pop_c50 = 1, pop_tover = 1, pd_Ω = 6, σ_fx = 1), 6, (:Depot, :Central, :Turnover), (:CL, :Vc, :tabs, :Ka, :e0, :emax, :c50, :tover, :kout, :rin), ParamSet{@NamedTuple{pop_CL::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_V::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tabs::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_lag::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pk_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_prop::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, σ_add::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_e0::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_emax::RealDomain{TransformVariables.Infinity{false}, TransformVariables.Infinity{true}, Float64}, pop_c50::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tover::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pd_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_fx::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}, ParamSet{@NamedTuple{pk_η::Pumas.RandomEffect{false, var"#1#12"}, pd_η::Pumas.RandomEffect{false, var"#2#13"}}}, Pumas.PreObj{Pumas.TimeDispatcher{var"#3#14", var"#4#15"}}, Pumas.DCPObj{var"#6#17"}, var"#7#18", SciMLBase.ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, Nothing, SciMLBase.ODEFunction{false, SciMLBase.FullSpecialize, ModelingToolkit.GeneratedFunctionWrapper{(2, 3, false), var"#8#19", var"#9#20"}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ODESystem, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardODEProblem}, Pumas.DerivedObj{(:conc, :pca), false, var"#10#21"}, var"#11#22", ModelingToolkit.ODESystem, Pumas.PumasModelOptions}, Subject{@NamedTuple{conc::Vector{Union{Missing, Float64}}, pca::Vector{Union{Missing, Float64}}}, Pumas.ConstantCovar{@NamedTuple{FSZV::Float64, FSZCL::Float64}}, DosageRegimen{Float64, Union{Missing, Int64}, Float64, Float64, Float64, Float64}, Vector{Float64}}, @NamedTuple{ss_abstol::Float64, ss_reltol::Float64, ss_maxiters::Int64, abstol::Float64, reltol::Float64, alg::CompositeAlgorithm{0, Tuple{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}}, OrdinaryDiffEqCore.AutoSwitch{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}, Rational{Int64}, Int64}}, continuity::Symbol}, Matrix{Float64}, @NamedTuple{pop_CL::Float64, pop_V::Float64, pop_tabs::Float64, pop_lag::Float64, pk_Ω::PDMats.PDMat{Float64, Matrix{Float64}}, σ_prop::Float64, σ_add::Float64, pop_e0::Float64, pop_emax::Float64, pop_c50::Float64, pop_tover::Float64, pd_Ω::PDMats.PDMat{Float64, Matrix{Float64}}, σ_fx::Float64}}}}}, ReverseDiff.JacobianConfig{ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}, Vector{ReverseDiff.TrackedReal{Float64, Float64, Nothing}}}, Nothing}}, ADTypes.AutoForwardDiff{nothing, Nothing}, ADTypes.AutoReverseDiff{false}, ADTypes.AutoForwardDiff{nothing, Nothing}, Vector{Matrix{Float64}}}, Optim.MultivariateOptimizationResults{Optim.LBFGS{Nothing, LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Returns{Nothing}}, Vector{Float64}, Float64, Float64, Vector{Optim.OptimizationState{Float64, Optim.LBFGS{Nothing, LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Returns{Nothing}}}}, @NamedTuple{x_converged::Bool, f_converged::Bool, g_converged::Bool, f_limit_reached::Bool, g_limit_reached::Bool, h_limit_reached::Bool, time_limit::Bool, callback::Bool, f_increased::Bool, ls_failed::Bool, iterations::Bool, small_trustregion_radius::Bool}}, Vector{Float64}}, VEM{Optim.LBFGS{Nothing, LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Returns{Nothing}}, @NamedTuple{show_every::Int64}, @NamedTuple{}, Random.TaskLocalRNG, EnsembleThreads, ADTypes.AutoForwardDiff{nothing, Nothing}}, Vector{SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}}, @NamedTuple{diffeq_options::@NamedTuple{ss_abstol::Float64, ss_reltol::Float64, ss_maxiters::Int64, abstol::Float64, reltol::Float64, alg::CompositeAlgorithm{0, Tuple{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}}, OrdinaryDiffEqCore.AutoSwitch{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}, Rational{Int64}, Int64}}, continuity::Symbol}, ensemblealg::EnsembleThreads, constantcoef::Tuple{}}, ParamSet{@NamedTuple{pop_CL::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_V::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tabs::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_lag::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pk_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_prop::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, σ_add::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_e0::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_emax::RealDomain{TransformVariables.Infinity{false}, TransformVariables.Infinity{true}, Float64}, pop_c50::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tover::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pd_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_fx::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}, Pumas.OptimState{NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Float64, Vector{Float64}}, Optim.LBFGSState{Vector{Float64}, Vector{Vector{Float64}}, Vector{Vector{Float64}}, Float64, Vector{Float64}}}}; level::Float64, rethrow_error::Bool, ensemblealg::EnsembleThreads, sandwich_estimator::Bool)
   @ Pumas.Experimental ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/GZeMg/src/estimation/vem.jl:2596
 [2] infer(fpm::Pumas.FittedPumasModel{PumasModel{(pop_CL = 1, pop_V = 1, pop_tabs = 1, pop_lag = 1, pk_Ω = 6, σ_prop = 1, σ_add = 1, pop_e0 = 1, pop_emax = 1, pop_c50 = 1, pop_tover = 1, pd_Ω = 6, σ_fx = 1), 6, (:Depot, :Central, :Turnover), (:CL, :Vc, :tabs, :Ka, :e0, :emax, :c50, :tover, :kout, :rin), ParamSet{@NamedTuple{pop_CL::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_V::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tabs::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_lag::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pk_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_prop::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, σ_add::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_e0::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_emax::RealDomain{TransformVariables.Infinity{false}, TransformVariables.Infinity{true}, Float64}, pop_c50::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tover::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pd_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_fx::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}, ParamSet{@NamedTuple{pk_η::Pumas.RandomEffect{false, var"#1#12"}, pd_η::Pumas.RandomEffect{false, var"#2#13"}}}, Pumas.PreObj{Pumas.TimeDispatcher{var"#3#14", var"#4#15"}}, Pumas.DCPObj{var"#6#17"}, var"#7#18", SciMLBase.ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, Nothing, SciMLBase.ODEFunction{false, SciMLBase.FullSpecialize, ModelingToolkit.GeneratedFunctionWrapper{(2, 3, false), var"#8#19", var"#9#20"}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ODESystem, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardODEProblem}, Pumas.DerivedObj{(:conc, :pca), false, var"#10#21"}, var"#11#22", ModelingToolkit.ODESystem, Pumas.PumasModelOptions}, Vector{Subject{@NamedTuple{conc::Vector{Union{Missing, Float64}}, pca::Vector{Union{Missing, Float64}}}, Pumas.ConstantCovar{@NamedTuple{FSZV::Float64, FSZCL::Float64}}, DosageRegimen{Float64, Union{Missing, Int64}, Float64, Float64, Float64, Float64}, Vector{Float64}}}, Pumas.Experimental.VEMResult{Pumas.Experimental.ThreadedVEMObj{true, false, true, true, PumasModel{(pop_CL = 1, pop_V = 1, pop_tabs = 1, pop_lag = 1, pk_Ω = 6, σ_prop = 1, σ_add = 1, pop_e0 = 1, pop_emax = 1, pop_c50 = 1, pop_tover = 1, pd_Ω = 6, σ_fx = 1), 6, (:Depot, :Central, :Turnover), (:CL, :Vc, :tabs, :Ka, :e0, :emax, :c50, :tover, :kout, :rin), ParamSet{@NamedTuple{pop_CL::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_V::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tabs::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_lag::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pk_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_prop::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, σ_add::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_e0::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_emax::RealDomain{TransformVariables.Infinity{false}, TransformVariables.Infinity{true}, Float64}, pop_c50::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tover::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pd_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_fx::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}, ParamSet{@NamedTuple{pk_η::Pumas.RandomEffect{false, var"#1#12"}, pd_η::Pumas.RandomEffect{false, var"#2#13"}}}, Pumas.PreObj{Pumas.TimeDispatcher{var"#3#14", var"#4#15"}}, Pumas.DCPObj{var"#6#17"}, var"#7#18", SciMLBase.ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, Nothing, SciMLBase.ODEFunction{false, SciMLBase.FullSpecialize, ModelingToolkit.GeneratedFunctionWrapper{(2, 3, false), var"#8#19", var"#9#20"}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ODESystem, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardODEProblem}, Pumas.DerivedObj{(:conc, :pca), false, var"#10#21"}, var"#11#22", ModelingToolkit.ODESystem, Pumas.PumasModelOptions}, Vector{Subject{@NamedTuple{conc::Vector{Union{Missing, Float64}}, pca::Vector{Union{Missing, Float64}}}, Pumas.ConstantCovar{@NamedTuple{FSZV::Float64, FSZCL::Float64}}, DosageRegimen{Float64, Union{Missing, Int64}, Float64, Float64, Float64, Float64}, Vector{Float64}}}, @NamedTuple{pop_CL::Float64, pop_V::Float64, pop_tabs::Float64, pop_lag::Float64, pk_Ω::PDMats.PDiagMat{Float64, Vector{Float64}}, σ_prop::Float64, σ_add::Float64, pop_e0::Float64, pop_emax::Float64, pop_c50::Float64, pop_tover::Float64, pd_Ω::PDMats.PDiagMat{Float64, Vector{Float64}}, σ_fx::Float64}, Random.TaskLocalRNG, Nothing, Nothing, @NamedTuple{ss_abstol::Float64, ss_reltol::Float64, ss_maxiters::Int64, abstol::Float64, reltol::Float64, alg::CompositeAlgorithm{0, Tuple{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}}, OrdinaryDiffEqCore.AutoSwitch{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}, Rational{Int64}, Int64}}, continuity::Symbol}, Vector{Array{Float64, 3}}, TransformVariables.TransformTuple{@NamedTuple{pop_CL::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_V::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tabs::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_lag::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pk_Ω::Pumas.PSDTransform, σ_prop::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, σ_add::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_e0::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_emax::TransformVariables.Identity, pop_c50::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tover::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pd_Ω::Pumas.PSDTransform, σ_fx::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}}}, TransformVariables.TransformTuple{@NamedTuple{pop_CL::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_V::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tabs::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_lag::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pk_Ω::TransformVariables.Constant{Nothing}, σ_prop::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, σ_add::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_e0::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_emax::TransformVariables.Identity, pop_c50::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tover::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pd_Ω::TransformVariables.Constant{Nothing}, σ_fx::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}}}, Vector{Vector{Float64}}, Vector{Vector{Float64}}, DifferentiationInterfaceForwardDiffExt.ForwardDiffGradientPrep{Tuple{Pumas._L_param{false, (:pop_CL, :pop_V, :pop_tabs, :pop_lag, :pk_Ω, :σ_prop, :σ_add, :pop_e0, :pop_emax, :pop_c50, :pop_tover, :pd_Ω, :σ_fx), TransformVariables.TransformTuple{@NamedTuple{pop_CL::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_V::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tabs::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_lag::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pk_Ω::Pumas.PSDTransform, σ_prop::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, σ_add::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_e0::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_emax::TransformVariables.Identity, pop_c50::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tover::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pd_Ω::Pumas.PSDTransform, σ_fx::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}}}, ParamSet{@NamedTuple{pop_CL::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_V::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tabs::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_lag::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pk_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_prop::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, σ_add::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_e0::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_emax::RealDomain{TransformVariables.Infinity{false}, TransformVariables.Infinity{true}, Float64}, pop_c50::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tover::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pd_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_fx::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}}, ADTypes.AutoForwardDiff{nothing, Nothing}, SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}, Tuple{}}, ForwardDiff.GradientConfig{ForwardDiff.Tag{Pumas._L_param{false, (:pop_CL, :pop_V, :pop_tabs, :pop_lag, :pk_Ω, :σ_prop, :σ_add, :pop_e0, :pop_emax, :pop_c50, :pop_tover, :pd_Ω, :σ_fx), TransformVariables.TransformTuple{@NamedTuple{pop_CL::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_V::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tabs::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_lag::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pk_Ω::Pumas.PSDTransform, σ_prop::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, σ_add::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_e0::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_emax::TransformVariables.Identity, pop_c50::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tover::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pd_Ω::Pumas.PSDTransform, σ_fx::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}}}, ParamSet{@NamedTuple{pop_CL::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_V::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tabs::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_lag::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pk_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_prop::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, σ_add::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_e0::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_emax::RealDomain{TransformVariables.Infinity{false}, TransformVariables.Infinity{true}, Float64}, pop_c50::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tover::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pd_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_fx::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}}, Float64}, Float64, 12, Vector{ForwardDiff.Dual{ForwardDiff.Tag{Pumas._L_param{false, (:pop_CL, :pop_V, :pop_tabs, :pop_lag, :pk_Ω, :σ_prop, :σ_add, :pop_e0, :pop_emax, :pop_c50, :pop_tover, :pd_Ω, :σ_fx), TransformVariables.TransformTuple{@NamedTuple{pop_CL::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_V::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tabs::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_lag::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pk_Ω::Pumas.PSDTransform, σ_prop::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, σ_add::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_e0::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_emax::TransformVariables.Identity, pop_c50::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tover::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pd_Ω::Pumas.PSDTransform, σ_fx::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}}}, ParamSet{@NamedTuple{pop_CL::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_V::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tabs::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_lag::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pk_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_prop::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, σ_add::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_e0::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_emax::RealDomain{TransformVariables.Infinity{false}, TransformVariables.Infinity{true}, Float64}, pop_c50::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tover::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pd_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_fx::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}}, Float64}, Float64, 12}}}, Tuple{}}, Vector{DifferentiationInterfaceReverseDiffExt.ReverseDiffGradientPrep{Tuple{typeof(Pumas.Experimental.G1_prior), ADTypes.AutoReverseDiff{false}, Vector{Float64}, Tuple{DifferentiationInterface.Constant{Pumas.Experimental._L_rfx2{true, (:pop_CL, :pop_V, :pop_tabs, :pop_lag, :pk_Ω, :σ_prop, :σ_add, :pop_e0, :pop_emax, :pop_c50, :pop_tover, :pd_Ω, :σ_fx), TransformVariables.TransformTuple{@NamedTuple{pop_CL::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_V::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tabs::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_lag::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pk_Ω::Pumas.PSDTransform, σ_prop::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, σ_add::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_e0::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_emax::TransformVariables.Identity, pop_c50::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tover::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pd_Ω::Pumas.PSDTransform, σ_fx::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}}}, TransformVariables.TransformTuple{@NamedTuple{pop_CL::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_V::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tabs::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_lag::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pk_Ω::TransformVariables.Constant{Nothing}, σ_prop::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, σ_add::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_e0::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_emax::TransformVariables.Identity, pop_c50::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tover::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pd_Ω::TransformVariables.Constant{Nothing}, σ_fx::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}}}, PumasModel{(pop_CL = 1, pop_V = 1, pop_tabs = 1, pop_lag = 1, pk_Ω = 6, σ_prop = 1, σ_add = 1, pop_e0 = 1, pop_emax = 1, pop_c50 = 1, pop_tover = 1, pd_Ω = 6, σ_fx = 1), 6, (:Depot, :Central, :Turnover), (:CL, :Vc, :tabs, :Ka, :e0, :emax, :c50, :tover, :kout, :rin), ParamSet{@NamedTuple{pop_CL::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_V::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tabs::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_lag::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pk_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_prop::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, σ_add::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_e0::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_emax::RealDomain{TransformVariables.Infinity{false}, TransformVariables.Infinity{true}, Float64}, pop_c50::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tover::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pd_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_fx::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}, ParamSet{@NamedTuple{pk_η::Pumas.RandomEffect{false, var"#1#12"}, pd_η::Pumas.RandomEffect{false, var"#2#13"}}}, Pumas.PreObj{Pumas.TimeDispatcher{var"#3#14", var"#4#15"}}, Pumas.DCPObj{var"#6#17"}, var"#7#18", SciMLBase.ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, Nothing, SciMLBase.ODEFunction{false, SciMLBase.FullSpecialize, ModelingToolkit.GeneratedFunctionWrapper{(2, 3, false), var"#8#19", var"#9#20"}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ODESystem, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardODEProblem}, Pumas.DerivedObj{(:conc, :pca), false, var"#10#21"}, var"#11#22", ModelingToolkit.ODESystem, Pumas.PumasModelOptions}, Subject{@NamedTuple{conc::Vector{Union{Missing, Float64}}, pca::Vector{Union{Missing, Float64}}}, Pumas.ConstantCovar{@NamedTuple{FSZV::Float64, FSZCL::Float64}}, DosageRegimen{Float64, Union{Missing, Int64}, Float64, Float64, Float64, Float64}, Vector{Float64}}, @NamedTuple{ss_abstol::Float64, ss_reltol::Float64, ss_maxiters::Int64, abstol::Float64, reltol::Float64, alg::CompositeAlgorithm{0, Tuple{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}}, OrdinaryDiffEqCore.AutoSwitch{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}, Rational{Int64}, Int64}}, continuity::Symbol}, Matrix{Float64}, @NamedTuple{pop_CL::Float64, pop_V::Float64, pop_tabs::Float64, pop_lag::Float64, pk_Ω::PDMats.PDMat{Float64, Matrix{Float64}}, σ_prop::Float64, σ_add::Float64, pop_e0::Float64, pop_emax::Float64, pop_c50::Float64, pop_tover::Float64, pd_Ω::PDMats.PDMat{Float64, Matrix{Float64}}, σ_fx::Float64}}}}}, ReverseDiff.GradientConfig{ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}}, Nothing}}, Vector{DifferentiationInterfaceForwardDiffExt.ForwardDiffGradientPrep{Tuple{typeof(Pumas.Experimental.G1_likelihood), ADTypes.AutoForwardDiff{nothing, Nothing}, Vector{Float64}, Tuple{DifferentiationInterface.Constant{Pumas.Experimental._L_rfx2{true, (:pop_CL, :pop_V, :pop_tabs, :pop_lag, :pk_Ω, :σ_prop, :σ_add, :pop_e0, :pop_emax, :pop_c50, :pop_tover, :pd_Ω, :σ_fx), TransformVariables.TransformTuple{@NamedTuple{pop_CL::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_V::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tabs::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_lag::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pk_Ω::Pumas.PSDTransform, σ_prop::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, σ_add::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_e0::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_emax::TransformVariables.Identity, pop_c50::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tover::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pd_Ω::Pumas.PSDTransform, σ_fx::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}}}, TransformVariables.TransformTuple{@NamedTuple{pop_CL::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_V::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tabs::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_lag::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pk_Ω::TransformVariables.Constant{Nothing}, σ_prop::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, σ_add::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_e0::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_emax::TransformVariables.Identity, pop_c50::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tover::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pd_Ω::TransformVariables.Constant{Nothing}, σ_fx::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}}}, PumasModel{(pop_CL = 1, pop_V = 1, pop_tabs = 1, pop_lag = 1, pk_Ω = 6, σ_prop = 1, σ_add = 1, pop_e0 = 1, pop_emax = 1, pop_c50 = 1, pop_tover = 1, pd_Ω = 6, σ_fx = 1), 6, (:Depot, :Central, :Turnover), (:CL, :Vc, :tabs, :Ka, :e0, :emax, :c50, :tover, :kout, :rin), ParamSet{@NamedTuple{pop_CL::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_V::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tabs::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_lag::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pk_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_prop::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, σ_add::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_e0::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_emax::RealDomain{TransformVariables.Infinity{false}, TransformVariables.Infinity{true}, Float64}, pop_c50::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tover::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pd_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_fx::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}, ParamSet{@NamedTuple{pk_η::Pumas.RandomEffect{false, var"#1#12"}, pd_η::Pumas.RandomEffect{false, var"#2#13"}}}, Pumas.PreObj{Pumas.TimeDispatcher{var"#3#14", var"#4#15"}}, Pumas.DCPObj{var"#6#17"}, var"#7#18", SciMLBase.ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, Nothing, SciMLBase.ODEFunction{false, SciMLBase.FullSpecialize, ModelingToolkit.GeneratedFunctionWrapper{(2, 3, false), var"#8#19", var"#9#20"}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ODESystem, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardODEProblem}, Pumas.DerivedObj{(:conc, :pca), false, var"#10#21"}, var"#11#22", ModelingToolkit.ODESystem, Pumas.PumasModelOptions}, Subject{@NamedTuple{conc::Vector{Union{Missing, Float64}}, pca::Vector{Union{Missing, Float64}}}, Pumas.ConstantCovar{@NamedTuple{FSZV::Float64, FSZCL::Float64}}, DosageRegimen{Float64, Union{Missing, Int64}, Float64, Float64, Float64, Float64}, Vector{Float64}}, @NamedTuple{ss_abstol::Float64, ss_reltol::Float64, ss_maxiters::Int64, abstol::Float64, reltol::Float64, alg::CompositeAlgorithm{0, Tuple{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}}, OrdinaryDiffEqCore.AutoSwitch{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}, Rational{Int64}, Int64}}, continuity::Symbol}, Matrix{Float64}, @NamedTuple{pop_CL::Float64, pop_V::Float64, pop_tabs::Float64, pop_lag::Float64, pk_Ω::PDMats.PDMat{Float64, Matrix{Float64}}, σ_prop::Float64, σ_add::Float64, pop_e0::Float64, pop_emax::Float64, pop_c50::Float64, pop_tover::Float64, pd_Ω::PDMats.PDMat{Float64, Matrix{Float64}}, σ_fx::Float64}}}}}, ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(Pumas.Experimental.G1_likelihood), Float64}, Float64, 9, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(Pumas.Experimental.G1_likelihood), Float64}, Float64, 9}}}, Tuple{Nothing}}}, Vector{DifferentiationInterfaceReverseDiffExt.ReverseDiffTwoArgJacobianPrep{Tuple{typeof(Pumas.Experimental.get_randeffs!), SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}, ADTypes.AutoReverseDiff{false}, Vector{Float64}, Tuple{DifferentiationInterface.Constant{Pumas.Experimental._L_rfx2{true, (:pop_CL, :pop_V, :pop_tabs, :pop_lag, :pk_Ω, :σ_prop, :σ_add, :pop_e0, :pop_emax, :pop_c50, :pop_tover, :pd_Ω, :σ_fx), TransformVariables.TransformTuple{@NamedTuple{pop_CL::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_V::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tabs::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_lag::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pk_Ω::Pumas.PSDTransform, σ_prop::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, σ_add::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_e0::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_emax::TransformVariables.Identity, pop_c50::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tover::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pd_Ω::Pumas.PSDTransform, σ_fx::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}}}, TransformVariables.TransformTuple{@NamedTuple{pop_CL::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_V::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tabs::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_lag::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pk_Ω::TransformVariables.Constant{Nothing}, σ_prop::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, σ_add::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_e0::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_emax::TransformVariables.Identity, pop_c50::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pop_tover::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}, pd_Ω::TransformVariables.Constant{Nothing}, σ_fx::TransformVariables.CompositeScalarTransform{Tuple{TransformVariables.TVShift{Float64}, TransformVariables.TVExp}}}}, PumasModel{(pop_CL = 1, pop_V = 1, pop_tabs = 1, pop_lag = 1, pk_Ω = 6, σ_prop = 1, σ_add = 1, pop_e0 = 1, pop_emax = 1, pop_c50 = 1, pop_tover = 1, pd_Ω = 6, σ_fx = 1), 6, (:Depot, :Central, :Turnover), (:CL, :Vc, :tabs, :Ka, :e0, :emax, :c50, :tover, :kout, :rin), ParamSet{@NamedTuple{pop_CL::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_V::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tabs::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_lag::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pk_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_prop::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, σ_add::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_e0::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_emax::RealDomain{TransformVariables.Infinity{false}, TransformVariables.Infinity{true}, Float64}, pop_c50::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tover::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pd_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_fx::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}, ParamSet{@NamedTuple{pk_η::Pumas.RandomEffect{false, var"#1#12"}, pd_η::Pumas.RandomEffect{false, var"#2#13"}}}, Pumas.PreObj{Pumas.TimeDispatcher{var"#3#14", var"#4#15"}}, Pumas.DCPObj{var"#6#17"}, var"#7#18", SciMLBase.ODEProblem{Nothing, Tuple{Nothing, Nothing}, false, Nothing, SciMLBase.ODEFunction{false, SciMLBase.FullSpecialize, ModelingToolkit.GeneratedFunctionWrapper{(2, 3, false), var"#8#19", var"#9#20"}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ODESystem, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardODEProblem}, Pumas.DerivedObj{(:conc, :pca), false, var"#10#21"}, var"#11#22", ModelingToolkit.ODESystem, Pumas.PumasModelOptions}, Subject{@NamedTuple{conc::Vector{Union{Missing, Float64}}, pca::Vector{Union{Missing, Float64}}}, Pumas.ConstantCovar{@NamedTuple{FSZV::Float64, FSZCL::Float64}}, DosageRegimen{Float64, Union{Missing, Int64}, Float64, Float64, Float64, Float64}, Vector{Float64}}, @NamedTuple{ss_abstol::Float64, ss_reltol::Float64, ss_maxiters::Int64, abstol::Float64, reltol::Float64, alg::CompositeAlgorithm{0, Tuple{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}}, OrdinaryDiffEqCore.AutoSwitch{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}, Rational{Int64}, Int64}}, continuity::Symbol}, Matrix{Float64}, @NamedTuple{pop_CL::Float64, pop_V::Float64, pop_tabs::Float64, pop_lag::Float64, pk_Ω::PDMats.PDMat{Float64, Matrix{Float64}}, σ_prop::Float64, σ_add::Float64, pop_e0::Float64, pop_emax::Float64, pop_c50::Float64, pop_tover::Float64, pd_Ω::PDMats.PDMat{Float64, Matrix{Float64}}, σ_fx::Float64}}}}}, ReverseDiff.JacobianConfig{ReverseDiff.TrackedArray{Float64, Float64, 1, Vector{Float64}, Vector{Float64}}, Vector{ReverseDiff.TrackedReal{Float64, Float64, Nothing}}}, Nothing}}, ADTypes.AutoForwardDiff{nothing, Nothing}, ADTypes.AutoReverseDiff{false}, ADTypes.AutoForwardDiff{nothing, Nothing}, Vector{Matrix{Float64}}}, Optim.MultivariateOptimizationResults{Optim.LBFGS{Nothing, LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Returns{Nothing}}, Vector{Float64}, Float64, Float64, Vector{Optim.OptimizationState{Float64, Optim.LBFGS{Nothing, LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Returns{Nothing}}}}, @NamedTuple{x_converged::Bool, f_converged::Bool, g_converged::Bool, f_limit_reached::Bool, g_limit_reached::Bool, h_limit_reached::Bool, time_limit::Bool, callback::Bool, f_increased::Bool, ls_failed::Bool, iterations::Bool, small_trustregion_radius::Bool}}, Vector{Float64}}, VEM{Optim.LBFGS{Nothing, LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Returns{Nothing}}, @NamedTuple{show_every::Int64}, @NamedTuple{}, Random.TaskLocalRNG, EnsembleThreads, ADTypes.AutoForwardDiff{nothing, Nothing}}, Vector{SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}}, @NamedTuple{diffeq_options::@NamedTuple{ss_abstol::Float64, ss_reltol::Float64, ss_maxiters::Int64, abstol::Float64, reltol::Float64, alg::CompositeAlgorithm{0, Tuple{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}}, OrdinaryDiffEqCore.AutoSwitch{Vern7{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Rodas5P{1, ADTypes.AutoForwardDiff{1, Nothing}, LinearSolve.GenericLUFactorization{LinearAlgebra.RowMaximum}, typeof(OrdinaryDiffEqCore.DEFAULT_PRECS), Val{:forward}(), true, nothing, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!)}, Rational{Int64}, Int64}}, continuity::Symbol}, ensemblealg::EnsembleThreads, constantcoef::Tuple{}}, ParamSet{@NamedTuple{pop_CL::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_V::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tabs::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_lag::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pk_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_prop::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, σ_add::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_e0::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_emax::RealDomain{TransformVariables.Infinity{false}, TransformVariables.Infinity{true}, Float64}, pop_c50::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pop_tover::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}, pd_Ω::PSDDomain{PDMats.PDiagMat{Float64, Vector{Float64}}}, σ_fx::RealDomain{Float64, TransformVariables.Infinity{true}, Float64}}}, Pumas.OptimState{NLSolversBase.OnceDifferentiable{Float64, Vector{Float64}, Float64, Vector{Float64}}, Optim.LBFGSState{Vector{Float64}, Vector{Vector{Float64}}, Vector{Vector{Float64}}, Float64, Vector{Float64}}}})
   @ Pumas.Experimental ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/GZeMg/src/estimation/vem.jl:2589
 [3] top-level scope
   @ ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/tutorials/estimation/01-vem-introduction.qmd:405

The recommended workaround is to use the VEM estimates as a starting point for a FOCE/LaplaceI fit, then run inference on the FOCE/LaplaceI result. This typically converges very quickly since VEM provides excellent initial values:

@time fpm_foce_from_vem = fit(
    warfarin_model,
    pop,
    coef(fpm_vem),
    FOCE();
    init_randeffs = empirical_bayes(fpm_vem),
    optim_options = (; show_every = 20),
)
Iter     Function value   Gradient norm 
     0     1.072198e+03     1.238181e+01
 * time: 7.104873657226562e-5
    20     1.069323e+03     3.312721e+01
 * time: 80.10134983062744
    40     1.061185e+03     3.274055e+01
 * time: 115.6318678855896
135.722583 seconds (120.42 M allocations: 30.401 GiB, 4.19% gc time, 0.22% compilation time)
FittedPumasModel

Dynamical system type:               Nonlinear ODE
Solver(s): (OrdinaryDiffEqVerner.Vern7,OrdinaryDiffEqRosenbrock.Rodas5P)

Number of subjects:                             32

Observation records:         Active        Missing
    conc:                       251             47
    pca:                        232             66
    Total:                      483            113

Number of parameters:      Constant      Optimized
                                  0             23

Likelihood approximation:                     FOCE
Likelihood optimizer:                         BFGS

Termination Reason:                   GradientNorm
Log-likelihood value:                   -1060.6836

------------------------
            Estimate
------------------------
pop_CL       0.13498
pop_V        8.0347
pop_tabs     0.54029
pop_lag      0.87451
pk_Ω₁,₁      0.06856
pk_Ω₂,₁      0.007987
pk_Ω₃,₁      0.032911
pk_Ω₂,₂      0.020884
pk_Ω₃,₂     -0.0036756
pk_Ω₃,₃      0.83843
σ_prop       0.088881
σ_add        0.40867
pop_e0      96.361
pop_emax    -1.0878
pop_c50      1.6319
pop_tover   14.423
pd_Ω₁,₁      0.0028601
pd_Ω₂,₁      0.0017294
pd_Ω₃,₁      0.00052742
pd_Ω₂,₂      0.14354
pd_Ω₃,₂     -0.036232
pd_Ω₃,₃      0.017757
σ_fx         3.5709
------------------------
loglikelihood(fpm_foce_from_vem)
-1060.6835570268515
infer(fpm_foce_from_vem)
[ Info: Calculating: variance-covariance matrix.
[ Info: Done.

Asymptotic inference results using sandwich estimator

Dynamical system type:               Nonlinear ODE
Solver(s): (OrdinaryDiffEqVerner.Vern7,OrdinaryDiffEqRosenbrock.Rodas5P)

Number of subjects:                             32

Observation records:         Active        Missing
    conc:                       251             47
    pca:                        232             66
    Total:                      483            113

Number of parameters:      Constant      Optimized
                                  0             23

Likelihood approximation:                     FOCE
Likelihood optimizer:                         BFGS

Termination Reason:                   GradientNorm
Log-likelihood value:                   -1060.6836

-----------------------------------------------------------------
            Estimate      SE          95.0% C.I.
-----------------------------------------------------------------
pop_CL       0.13498      0.0066379   [  0.12197   ;  0.14799  ]
pop_V        8.0347       0.21543     [  7.6125    ;  8.4569   ]
pop_tabs     0.54029      0.1776      [  0.19221   ;  0.88837  ]
pop_lag      0.87451      0.05685     [  0.76309   ;  0.98594  ]
pk_Ω₁,₁      0.06856      0.024778    [  0.019995  ;  0.11712  ]
pk_Ω₂,₁      0.007987     0.0072022   [ -0.006129  ;  0.022103 ]
pk_Ω₃,₁      0.032911     0.055627    [ -0.076116  ;  0.14194  ]
pk_Ω₂,₂      0.020884     0.0063164   [  0.0085036 ;  0.033264 ]
pk_Ω₃,₂     -0.0036756    0.059614    [ -0.12052   ;  0.11317  ]
pk_Ω₃,₃      0.83843      0.38633     [  0.081234  ;  1.5956   ]
σ_prop       0.088881     0.01387     [  0.061696  ;  0.11607  ]
σ_add        0.40867      0.069893    [  0.27169   ;  0.54566  ]
pop_e0      96.361        1.1784      [ 94.051     ; 98.67     ]
pop_emax    -1.0878       0.046168    [ -1.1783    ; -0.99734  ]
pop_c50      1.6319       0.33117     [  0.98282   ;  2.281    ]
pop_tover   14.423        0.76978     [ 12.914     ; 15.932    ]
pd_Ω₁,₁      0.0028601    0.0010039   [  0.00089247;  0.0048277]
pd_Ω₂,₁      0.0017294    0.0037375   [ -0.0055959 ;  0.0090548]
pd_Ω₃,₁      0.00052742   0.0022402   [ -0.0038633 ;  0.0049182]
pd_Ω₂,₂      0.14354      0.04335     [  0.058575  ;  0.2285   ]
pd_Ω₃,₂     -0.036232     0.0143      [ -0.06426   ; -0.0082043]
pd_Ω₃,₃      0.017757     0.0092593   [ -0.00039054;  0.035905 ]
σ_fx         3.5709       0.30611     [  2.9709    ;  4.1709   ]
-----------------------------------------------------------------
NoteTrue Posterior Mode vs Posterior Approximation Mode

Initializing a FOCE fit with VEM estimates and running for 0 or more iterations is also a way to get a fitted Pumas model with standard EBEs (true posterior mode, instead of the VEM posterior approximation’s mode) that can be used for downstream functions, such as predict, inspect and loglikelihood.

TipVEM as an Initializer

Using VEM to initialize FOCE is a powerful workflow. VEM provides good population parameter estimates and per-subject random effect posteriors, which together allow FOCE to converge quickly and reliably.

3 Comparing VEM and FOCE

3.1 Direct FOCE Fit

Let us also fit the model directly with FOCE from default initial values for comparison:

@time fpm_foce_direct = fit(
    warfarin_model,
    pop,
    init_params(warfarin_model),
    FOCE();
    optim_options = (; show_every = 20),
)
Iter     Function value   Gradient norm 
     0     5.225547e+06     1.026248e+07
 * time: 9.918212890625e-5
    20     4.290829e+03     1.936591e+03
 * time: 32.39234113693237
    40     1.622225e+03     5.112982e+02
 * time: 87.54792618751526
    60     1.380720e+03     3.611768e+02
 * time: 117.313725233078
    80     1.271772e+03     1.229157e+02
 * time: 145.53272700309753
   100     1.264818e+03     3.745549e+01
 * time: 172.6462550163269
   120     1.255091e+03     9.895804e+01
 * time: 198.96994304656982
   140     1.195682e+03     4.017472e+01
 * time: 230.9579951763153
   160     1.189284e+03     5.575148e+01
 * time: 261.5078411102295
   180     1.183513e+03     4.676253e+01
 * time: 290.3748950958252
   200     1.169818e+03     1.753033e+01
 * time: 322.53380513191223
   220     1.141433e+03     8.662493e+01
 * time: 357.9451320171356
   240     1.125162e+03     3.124826e+01
 * time: 405.4076101779938
   260     1.082971e+03     9.731992e+00
 * time: 458.8164231777191
   280     1.080982e+03     4.115973e+00
 * time: 505.6427872180939
   300     1.079409e+03     4.518367e+00
 * time: 554.656131029129
   320     1.071216e+03     7.824287e+00
 * time: 609.5648930072784
   340     1.064598e+03     6.447437e+00
 * time: 652.2475261688232
   360     1.060925e+03     2.996939e+01
 * time: 691.7175710201263
712.015951 seconds (363.20 M allocations: 139.020 GiB, 4.37% gc time, 0.03% compilation time)
FittedPumasModel

Dynamical system type:               Nonlinear ODE
Solver(s): (OrdinaryDiffEqVerner.Vern7,OrdinaryDiffEqRosenbrock.Rodas5P)

Number of subjects:                             32

Observation records:         Active        Missing
    conc:                       251             47
    pca:                        232             66
    Total:                      483            113

Number of parameters:      Constant      Optimized
                                  0             23

Likelihood approximation:                     FOCE
Likelihood optimizer:                         BFGS

Termination Reason:                      NoXChange
Log-likelihood value:                   -1060.6958

------------------------
            Estimate
------------------------
pop_CL       0.13506
pop_V        8.0358
pop_tabs     0.55913
pop_lag      0.87346
pk_Ω₁,₁      0.068711
pk_Ω₂,₁      0.0080307
pk_Ω₃,₁      0.034242
pk_Ω₂,₂      0.020889
pk_Ω₃,₂      0.0010522
pk_Ω₃,₃      0.83181
σ_prop       0.088868
σ_add        0.40866
pop_e0      96.364
pop_emax    -1.0873
pop_c50      1.6283
pop_tover   14.417
pd_Ω₁,₁      0.0028484
pd_Ω₂,₁      0.0017427
pd_Ω₃,₁      0.00054849
pd_Ω₂,₂      0.1444
pd_Ω₃,₂     -0.036343
pd_Ω₃,₃      0.017419
σ_fx         3.5752
------------------------
loglikelihood(fpm_foce_direct)
-1060.6957822778843

3.2 Compare Estimates

compare_estimates(; VEM = fpm_vem, VEM_FOCE = fpm_foce_from_vem, FOCE = fpm_foce_direct)
23×4 DataFrame
Row parameter VEM VEM_FOCE FOCE
String Float64? Float64? Float64?
1 pop_CL 0.134396 0.134975 0.135059
2 pop_V 8.01812 8.0347 8.03583
3 pop_tabs 0.407714 0.540286 0.559127
4 pop_lag 0.881609 0.874513 0.873459
5 pk_Ω₁,₁ 0.0658831 0.0685598 0.0687108
6 pk_Ω₂,₁ 0.00899818 0.00798703 0.00803068
7 pk_Ω₃,₁ 0.0346546 0.0329108 0.0342423
8 pk_Ω₂,₂ 0.0205447 0.0208836 0.0208894
9 pk_Ω₃,₂ -0.0269033 -0.00367559 0.00105221
10 pk_Ω₃,₃ 0.928079 0.838426 0.831814
11 σ_prop 0.0884536 0.0888807 0.0888683
12 σ_add 0.388876 0.408673 0.408663
13 pop_e0 96.4755 96.3606 96.3638
14 pop_emax -1.07136 -1.08782 -1.08733
15 pop_c50 1.53539 1.6319 1.62827
16 pop_tover 14.1607 14.4232 14.4168
17 pd_Ω₁,₁ 0.000319028 0.00286007 0.00284845
18 pd_Ω₂,₁ 0.00639603 0.00172944 0.00174267
19 pd_Ω₃,₁ -0.00138542 0.000527422 0.000548494
20 pd_Ω₂,₂ 0.129104 0.14354 0.144404
21 pd_Ω₃,₂ -0.0309213 -0.0362321 -0.036343
22 pd_Ω₃,₃ 0.0173689 0.0177573 0.0174188
23 σ_fx 4.33198 3.57089 3.57517

3.3 Compare Log-Likelihoods

DataFrame(
    Method = [
        "Deterministic VEM",
        "FOCE initialized from deterministic VEM",
        "FOCE direct",
    ],
    LogLikelihood = [
        loglikelihood(fpm_vem, FOCE()),
        loglikelihood(fpm_foce_from_vem),
        loglikelihood(fpm_foce_direct),
    ],
) |> simple_table
Method LogLikelihood
Deterministic VEM -1072.0
FOCE initialized from deterministic VEM -1061.0
FOCE direct -1061.0

4 VEM Options Reference

This section provides a comprehensive reference for all VEM keyword arguments. All options below are keyword arguments to VEM() unless otherwise noted.

4.1 Presampling vs. Resampling (presample)

The presample keyword controls which VEM variant is used. Many other defaults change depending on this choice:

Option Presampling (presample=true, default) Resampling (presample=false)
Optimizer optim_alg L-BFGS (deterministic) Adam (stochastic)
Phases nphases 1 10
Initial samples init_nsamples 15 3
Sample increment nsamples_increment half the initial samples (7 by default) the initial samples (3 by default)
Max samples max_nsamples Unlimited 25

Presampling draws samples once per phase and reuses them across iterations, which allows a deterministic optimizer (L-BFGS) to be used. Resampling draws fresh samples at every iteration, which requires a stochastic optimizer (Adam) but can be more appropriate for DeepPumas models with neural networks. However, generally the stochastic optimizer needs many more iterations to converge than the deterministic one, so presampling is a good default for most models.

Note that all of the above options can be overridden regardless of the presample setting.

# Resampling VEM
@time fpm_vem_resample = fit(
    warfarin_model,
    pop,
    init_params(warfarin_model),
    VEM(;
        presample = false,
        nphases = 5,
        optim_options = (; iterations = 1000, show_every = 20),
    ),
)
[ Info:           phase             iter       obj      ∇obj     s∇obj  
[ Info:           1 / 5         1 / 1000   1.4e+08   4.9e+08   1.0e+00
[ Info:           1 / 5        20 / 1000   4.6e+06   3.8e+07   7.8e-02
[ Info:           1 / 5        40 / 1000   1.2e+06   1.4e+07   2.9e-02
[ Info:           1 / 5        50 / 1000   7.1e+05   3.4e+06   7.0e-03
[ Info: Terminating phase due to sgtol.
[ Info:           2 / 5         1 / 1000   7.0e+05   3.5e+06   1.0e+00
[ Info:           2 / 5        20 / 1000   2.7e+05   1.9e+06   5.4e-01
[ Info:           2 / 5        40 / 1000   1.4e+05   3.5e+05   1.0e-01
[ Info:           2 / 5        60 / 1000   9.4e+04   2.1e+05   6.1e-02
[ Info:           2 / 5        75 / 1000   7.2e+04   1.5e+05   4.3e-02
[ Info: Terminating phase due to sgtol.
[ Info:           3 / 5         1 / 1000   7.3e+04   1.6e+05   1.0e+00
[ Info:           3 / 5        20 / 1000   6.1e+04   1.2e+05   7.7e-01
[ Info:           3 / 5        40 / 1000   5.2e+04   1.2e+05   7.3e-01
[ Info:           3 / 5        60 / 1000   4.4e+04   9.3e+04   5.8e-01
[ Info:           3 / 5        80 / 1000   3.9e+04   8.1e+04   5.1e-01
[ Info:           3 / 5       100 / 1000   3.4e+04   6.7e+04   4.2e-01
[ Info:           3 / 5       120 / 1000   3.0e+04   6.3e+04   3.9e-01
[ Info:           3 / 5       140 / 1000   2.7e+04   5.6e+04   3.5e-01
[ Info:           3 / 5       160 / 1000   2.4e+04   6.8e+04   4.3e-01
[ Info:           3 / 5       180 / 1000   2.2e+04   6.0e+04   3.7e-01
[ Info:           3 / 5       200 / 1000   2.1e+04   4.1e+04   2.6e-01
[ Info:           3 / 5       220 / 1000   1.9e+04   3.8e+04   2.4e-01
[ Info:           3 / 5       240 / 1000   1.7e+04   3.4e+04   2.1e-01
[ Info:           3 / 5       260 / 1000   1.6e+04   4.0e+04   2.5e-01
[ Info:           3 / 5       280 / 1000   1.5e+04   3.0e+04   1.9e-01
[ Info:           3 / 5       300 / 1000   1.4e+04   3.4e+04   2.1e-01
[ Info:           3 / 5       320 / 1000   1.3e+04   3.3e+04   2.1e-01
[ Info:           3 / 5       340 / 1000   1.3e+04   2.4e+04   1.5e-01
[ Info:           3 / 5       360 / 1000   1.2e+04   2.6e+04   1.6e-01
[ Info:           3 / 5       380 / 1000   1.1e+04   2.3e+04   1.5e-01
[ Info:           3 / 5       400 / 1000   1.1e+04   2.4e+04   1.5e-01
[ Info:           3 / 5       420 / 1000   1.0e+04   1.9e+04   1.2e-01
[ Info:           3 / 5       440 / 1000   9.8e+03   2.0e+04   1.3e-01
[ Info:           3 / 5       460 / 1000   9.2e+03   1.9e+04   1.2e-01
[ Info:           3 / 5       480 / 1000   8.9e+03   1.7e+04   1.1e-01
[ Info:           3 / 5       500 / 1000   8.5e+03   2.0e+04   1.2e-01
[ Info:           3 / 5       520 / 1000   8.1e+03   1.7e+04   1.0e-01
[ Info:           3 / 5       540 / 1000   7.9e+03   1.7e+04   1.1e-01
[ Info:           3 / 5       560 / 1000   7.5e+03   1.5e+04   9.3e-02
[ Info:           3 / 5       580 / 1000   7.3e+03   1.3e+04   8.3e-02
[ Info:           3 / 5       600 / 1000   7.0e+03   1.4e+04   8.7e-02
[ Info:           3 / 5       620 / 1000   6.7e+03   1.2e+04   7.7e-02
[ Info:           3 / 5       640 / 1000   6.5e+03   1.2e+04   7.2e-02
[ Info:           3 / 5       660 / 1000   6.3e+03   1.1e+04   6.9e-02
[ Info:           3 / 5       680 / 1000   6.1e+03   1.1e+04   6.9e-02
[ Info:           3 / 5       700 / 1000   6.0e+03   1.2e+04   7.7e-02
[ Info:           3 / 5       720 / 1000   5.7e+03   9.7e+03   6.1e-02
[ Info:           3 / 5       740 / 1000   5.6e+03   1.1e+04   6.6e-02
[ Info:           3 / 5       760 / 1000   5.4e+03   9.1e+03   5.7e-02
[ Info:           3 / 5       780 / 1000   5.3e+03   1.0e+04   6.3e-02
[ Info:           3 / 5       800 / 1000   5.1e+03   8.6e+03   5.4e-02
[ Info:           3 / 5       820 / 1000   5.0e+03   9.7e+03   6.1e-02
[ Info:           3 / 5       840 / 1000   4.9e+03   8.6e+03   5.4e-02
[ Info:           3 / 5       860 / 1000   4.8e+03   8.4e+03   5.3e-02
[ Info:           3 / 5       880 / 1000   4.6e+03   8.0e+03   5.0e-02
[ Info:           3 / 5       900 / 1000   4.5e+03   9.5e+03   6.0e-02
[ Info:           3 / 5       920 / 1000   4.4e+03   8.2e+03   5.2e-02
[ Info:           3 / 5       940 / 1000   4.3e+03   7.2e+03   4.5e-02
[ Info:           3 / 5       960 / 1000   4.2e+03   6.8e+03   4.3e-02
[ Info:           3 / 5       977 / 1000   4.1e+03   7.0e+03   4.4e-02
[ Info: Terminating phase due to sgtol.
[ Info:           4 / 5         1 / 1000   4.2e+03   6.5e+03   1.0e+00
[ Info:           4 / 5        20 / 1000   4.1e+03   6.6e+03   1.0e+00
[ Info:           4 / 5        40 / 1000   4.0e+03   7.5e+03   1.2e+00
[ Info:           4 / 5        60 / 1000   3.9e+03   6.0e+03   9.2e-01
[ Info:           4 / 5        80 / 1000   3.8e+03   5.8e+03   9.0e-01
[ Info:           4 / 5       100 / 1000   3.8e+03   5.8e+03   9.0e-01
[ Info:           4 / 5       120 / 1000   3.7e+03   6.3e+03   9.8e-01
[ Info:           4 / 5       140 / 1000   3.6e+03   5.5e+03   8.5e-01
[ Info:           4 / 5       160 / 1000   3.6e+03   5.3e+03   8.3e-01
[ Info:           4 / 5       180 / 1000   3.5e+03   5.3e+03   8.3e-01
[ Info:           4 / 5       200 / 1000   3.5e+03   5.3e+03   8.3e-01
[ Info:           4 / 5       220 / 1000   3.4e+03   5.3e+03   8.2e-01
[ Info:           4 / 5       240 / 1000   3.3e+03   5.3e+03   8.3e-01
[ Info:           4 / 5       260 / 1000   3.3e+03   4.9e+03   7.5e-01
[ Info:           4 / 5       280 / 1000   3.2e+03   4.5e+03   7.0e-01
[ Info:           4 / 5       300 / 1000   3.2e+03   4.8e+03   7.4e-01
[ Info:           4 / 5       320 / 1000   3.2e+03   4.5e+03   6.9e-01
[ Info:           4 / 5       340 / 1000   3.1e+03   4.5e+03   7.0e-01
[ Info:           4 / 5       360 / 1000   3.1e+03   4.9e+03   7.6e-01
[ Info:           4 / 5       380 / 1000   3.0e+03   4.1e+03   6.3e-01
[ Info:           4 / 5       400 / 1000   3.0e+03   4.0e+03   6.1e-01
[ Info:           4 / 5       420 / 1000   2.9e+03   4.0e+03   6.3e-01
[ Info:           4 / 5       440 / 1000   2.9e+03   4.1e+03   6.4e-01
[ Info:           4 / 5       460 / 1000   2.9e+03   4.1e+03   6.3e-01
[ Info:           4 / 5       480 / 1000   2.8e+03   5.2e+03   8.0e-01
[ Info:           4 / 5       500 / 1000   2.8e+03   3.6e+03   5.6e-01
[ Info:           4 / 5       520 / 1000   2.8e+03   3.8e+03   5.8e-01
[ Info:           4 / 5       540 / 1000   2.7e+03   3.4e+03   5.3e-01
[ Info:           4 / 5       560 / 1000   2.7e+03   3.5e+03   5.4e-01
[ Info:           4 / 5       580 / 1000   2.7e+03   3.3e+03   5.1e-01
[ Info:           4 / 5       600 / 1000   2.6e+03   3.3e+03   5.2e-01
[ Info:           4 / 5       620 / 1000   2.6e+03   3.9e+03   6.0e-01
[ Info:           4 / 5       640 / 1000   2.6e+03   3.1e+03   4.9e-01
[ Info:           4 / 5       660 / 1000   2.5e+03   3.1e+03   4.8e-01
[ Info:           4 / 5       680 / 1000   2.5e+03   3.2e+03   4.9e-01
[ Info:           4 / 5       700 / 1000   2.5e+03   3.2e+03   4.9e-01
[ Info:           4 / 5       720 / 1000   2.5e+03   3.1e+03   4.8e-01
[ Info:           4 / 5       740 / 1000   2.4e+03   2.9e+03   4.4e-01
[ Info:           4 / 5       760 / 1000   2.4e+03   2.9e+03   4.5e-01
[ Info:           4 / 5       780 / 1000   2.4e+03   2.7e+03   4.2e-01
[ Info:           4 / 5       800 / 1000   2.4e+03   2.8e+03   4.4e-01
[ Info:           4 / 5       820 / 1000   2.3e+03   2.7e+03   4.1e-01
[ Info:           4 / 5       840 / 1000   2.3e+03   3.2e+03   5.0e-01
[ Info:           4 / 5       860 / 1000   2.3e+03   2.9e+03   4.5e-01
[ Info:           4 / 5       880 / 1000   2.3e+03   2.6e+03   4.1e-01
[ Info:           4 / 5       900 / 1000   2.3e+03   3.2e+03   5.0e-01
[ Info:           4 / 5       920 / 1000   2.2e+03   3.2e+03   4.9e-01
[ Info:           4 / 5       940 / 1000   2.2e+03   2.6e+03   4.0e-01
[ Info:           4 / 5       960 / 1000   2.2e+03   3.2e+03   5.0e-01
[ Info:           4 / 5       980 / 1000   2.2e+03   2.4e+03   3.8e-01
[ Info:           4 / 5      1000 / 1000   2.2e+03   2.3e+03   3.5e-01
[ Info:           5 / 5         1 / 1000   2.2e+03   2.4e+03   1.0e+00
[ Info:           5 / 5        20 / 1000   2.1e+03   2.3e+03   9.6e-01
[ Info:           5 / 5        40 / 1000   2.2e+03   2.5e+03   1.0e+00
[ Info:           5 / 5        60 / 1000   2.1e+03   2.4e+03   1.0e+00
[ Info:           5 / 5        80 / 1000   2.1e+03   2.4e+03   9.9e-01
[ Info:           5 / 5       100 / 1000   2.1e+03   2.5e+03   1.1e+00
[ Info:           5 / 5       120 / 1000   2.1e+03   2.2e+03   9.2e-01
[ Info:           5 / 5       140 / 1000   2.1e+03   2.2e+03   9.4e-01
[ Info:           5 / 5       160 / 1000   2.1e+03   2.2e+03   9.2e-01
[ Info:           5 / 5       180 / 1000   2.1e+03   2.3e+03   9.5e-01
[ Info:           5 / 5       200 / 1000   2.1e+03   2.2e+03   9.2e-01
[ Info:           5 / 5       220 / 1000   2.1e+03   2.2e+03   9.1e-01
[ Info:           5 / 5       240 / 1000   2.1e+03   2.2e+03   9.0e-01
[ Info:           5 / 5       260 / 1000   2.1e+03   2.2e+03   9.0e-01
[ Info:           5 / 5       280 / 1000   2.1e+03   2.1e+03   8.9e-01
[ Info:           5 / 5       300 / 1000   2.1e+03   2.1e+03   8.9e-01
[ Info:           5 / 5       320 / 1000   2.1e+03   2.6e+03   1.1e+00
[ Info:           5 / 5       340 / 1000   2.1e+03   2.1e+03   9.0e-01
[ Info:           5 / 5       360 / 1000   2.1e+03   2.4e+03   9.9e-01
[ Info:           5 / 5       380 / 1000   2.1e+03   2.6e+03   1.1e+00
[ Info:           5 / 5       400 / 1000   2.1e+03   2.1e+03   8.6e-01
[ Info:           5 / 5       420 / 1000   2.1e+03   2.1e+03   9.0e-01
[ Info:           5 / 5       440 / 1000   2.0e+03   2.1e+03   8.7e-01
[ Info:           5 / 5       460 / 1000   2.1e+03   2.1e+03   8.6e-01
[ Info:           5 / 5       480 / 1000   2.0e+03   2.1e+03   8.7e-01
[ Info:           5 / 5       500 / 1000   2.0e+03   2.2e+03   9.0e-01
[ Info:           5 / 5       520 / 1000   2.0e+03   2.2e+03   9.2e-01
[ Info:           5 / 5       540 / 1000   2.0e+03   2.7e+03   1.1e+00
[ Info:           5 / 5       560 / 1000   2.0e+03   2.1e+03   8.9e-01
[ Info:           5 / 5       580 / 1000   2.0e+03   2.0e+03   8.3e-01
[ Info:           5 / 5       600 / 1000   2.0e+03   2.0e+03   8.3e-01
[ Info:           5 / 5       620 / 1000   2.0e+03   2.0e+03   8.2e-01
[ Info:           5 / 5       640 / 1000   2.0e+03   2.0e+03   8.1e-01
[ Info:           5 / 5       660 / 1000   2.0e+03   2.0e+03   8.5e-01
[ Info:           5 / 5       680 / 1000   2.0e+03   1.9e+03   8.1e-01
[ Info:           5 / 5       700 / 1000   2.0e+03   2.0e+03   8.3e-01
[ Info:           5 / 5       720 / 1000   2.0e+03   1.9e+03   8.0e-01
[ Info:           5 / 5       740 / 1000   2.0e+03   2.1e+03   8.9e-01
[ Info:           5 / 5       760 / 1000   2.0e+03   1.9e+03   7.9e-01
[ Info:           5 / 5       780 / 1000   2.0e+03   1.9e+03   8.0e-01
[ Info:           5 / 5       800 / 1000   2.0e+03   2.0e+03   8.4e-01
[ Info:           5 / 5       820 / 1000   2.0e+03   1.9e+03   8.1e-01
[ Info:           5 / 5       840 / 1000   2.0e+03   1.9e+03   7.8e-01
[ Info:           5 / 5       860 / 1000   2.0e+03   1.9e+03   7.8e-01
[ Info:           5 / 5       880 / 1000   2.0e+03   2.1e+03   8.7e-01
[ Info:           5 / 5       900 / 1000   1.9e+03   1.8e+03   7.7e-01
[ Info:           5 / 5       920 / 1000   1.9e+03   1.8e+03   7.7e-01
[ Info:           5 / 5       940 / 1000   1.9e+03   1.9e+03   7.7e-01
[ Info:           5 / 5       960 / 1000   1.9e+03   1.8e+03   7.6e-01
[ Info:           5 / 5       980 / 1000   1.9e+03   2.2e+03   9.1e-01
[ Info:           5 / 5      1000 / 1000   1.9e+03   1.9e+03   7.8e-01
1227.239509 seconds (1.95 G allocations: 424.287 GiB, 6.66% gc time, 0.38% compilation time)

FittedPumasModel

Dynamical system type:               Nonlinear ODE
Solver(s): (OrdinaryDiffEqVerner.Vern7,OrdinaryDiffEqRosenbrock.Rodas5P)

Number of subjects:                             32

Observation records:         Active        Missing
    conc:                       251             47
    pca:                        232             66
    Total:                      483            113

Number of parameters:      Constant      Optimized
                                  0             23

Estimation method:              VEM (EXPERIMENTAL)

-----------------------
            Estimate
-----------------------
pop_CL       0.15187
pop_V        8.4261
pop_tabs     0.68308
pop_lag      1.6572
pk_Ω₁,₁      0.042539
pk_Ω₂,₁     -0.068331
pk_Ω₃,₁      0.12408
pk_Ω₂,₂      0.12059
pk_Ω₃,₂     -0.084261
pk_Ω₃,₃      1.9341
σ_prop       0.13437
σ_add        2.8207
pop_e0      95.435
pop_emax    -1.0988
pop_c50      1.7177
pop_tover   13.207
pd_Ω₁,₁      0.0088784
pd_Ω₂,₁      0.012993
pd_Ω₃,₁     -0.0074459
pd_Ω₂,₂      0.17391
pd_Ω₃,₂     -0.054657
pd_Ω₃,₃      0.04374
σ_fx         0.75077
-----------------------
loglikelihood(fpm_vem_resample, FOCE())
-1986.6312948851134

The number of iterations and phases have been restricted in the above example for demonstration purposes.

NoteReading the VEM Convergence Trace

When presampling is enabled, the trace is deterministic and should show steady convergence. The trace printed during fitting has three columns:

Column Meaning
Iter Current iteration within the phase.
Function value The negative ELBO — the objective being minimized. Lower values indicate a better fit.
Gradient norm The gradient norm — how steep the objective landscape is. Decreasing values indicate the optimizer is converging.

When resampling is enabled, the trace is stochastic and may show more fluctuation. In the stochastic case, the trace printed during fitting has five columns:

Column Meaning
phase Current phase / total phases. Each phase increases the sample count.
iter Current iteration / max iterations within this phase.
obj The negative ELBO — the objective being minimized. Lower values indicate a better fit.
∇obj The gradient norm — how steep the objective landscape is. Decreasing values indicate the optimizer is converging.
s∇obj The scaled gradient norm — gradient norm relative to the first iteration of the current phase (normalized to 1.0). Shows relative convergence progress within each phase.

A phase may terminate early with “Terminating phase due to sgtol” when the scaled gradient norm drops below the stopping tolerance (0.05), triggering the next phase, typically with an increase in the sample size and a decrease in the learning rate. For stochastic optimizers like Adam (resampling VEM), some fluctuation in ∇obj is normal because the gradients are estimated from random samples.

compare_estimates(;
    VEM_presample = fpm_vem,
    VEM_resample = fpm_vem_resample,
    VEM_FOCE = fpm_foce_from_vem,
    FOCE = fpm_foce_direct,
)
23×5 DataFrame
Row parameter VEM_presample VEM_resample VEM_FOCE FOCE
String Float64? Float64? Float64? Float64?
1 pop_CL 0.134396 0.151872 0.134975 0.135059
2 pop_V 8.01812 8.42612 8.0347 8.03583
3 pop_tabs 0.407714 0.683085 0.540286 0.559127
4 pop_lag 0.881609 1.65724 0.874513 0.873459
5 pk_Ω₁,₁ 0.0658831 0.0425391 0.0685598 0.0687108
6 pk_Ω₂,₁ 0.00899818 -0.0683308 0.00798703 0.00803068
7 pk_Ω₃,₁ 0.0346546 0.124083 0.0329108 0.0342423
8 pk_Ω₂,₂ 0.0205447 0.120592 0.0208836 0.0208894
9 pk_Ω₃,₂ -0.0269033 -0.0842606 -0.00367559 0.00105221
10 pk_Ω₃,₃ 0.928079 1.93415 0.838426 0.831814
11 σ_prop 0.0884536 0.134371 0.0888807 0.0888683
12 σ_add 0.388876 2.82072 0.408673 0.408663
13 pop_e0 96.4755 95.435 96.3606 96.3638
14 pop_emax -1.07136 -1.09879 -1.08782 -1.08733
15 pop_c50 1.53539 1.71773 1.6319 1.62827
16 pop_tover 14.1607 13.2072 14.4232 14.4168
17 pd_Ω₁,₁ 0.000319028 0.00887835 0.00286007 0.00284845
18 pd_Ω₂,₁ 0.00639603 0.0129925 0.00172944 0.00174267
19 pd_Ω₃,₁ -0.00138542 -0.00744588 0.000527422 0.000548494
20 pd_Ω₂,₂ 0.129104 0.173908 0.14354 0.144404
21 pd_Ω₃,₂ -0.0309213 -0.0546573 -0.0362321 -0.036343
22 pd_Ω₃,₃ 0.0173689 0.0437396 0.0177573 0.0174188
23 σ_fx 4.33198 0.75077 3.57089 3.57517 
DataFrame(
    Method = [
        "Deterministic VEM",
        "Stochastic VEM",
        "FOCE initialized from deterministic VEM",
        "FOCE direct",
    ],
    LogLikelihood = [
        loglikelihood(fpm_vem, FOCE()),
        loglikelihood(fpm_vem_resample, FOCE()),
        loglikelihood(fpm_foce_from_vem),
        loglikelihood(fpm_foce_direct),
    ],
) |> simple_table
Method LogLikelihood
Deterministic VEM -1072.0
Stochastic VEM -1987.0
FOCE initialized from deterministic VEM -1061.0
FOCE direct -1061.0

Note that stochastic VEM took much longer to converge and it only got close to the best parameters. Deterministic VEM will generally converge faster. Stochastic VEM is mostly useful for models with neural networks and large populations where subject_ratio < 1 (to be demonstrated below) can be used to get quick estimates by running the algorithm for a few iterations only.

4.2 Sample Size Control

VEM draws samples from a base distribution to estimate the ELBO gradient. More samples improve accuracy but increase computation time.

When nsamples is set to a nonzero value, it overrides init_nsamples and sets nsamples_increment to 0, so the sample count stays fixed across all phases:

# Fixed 30 samples per subject across all phases
@time fpm_vem_nsamples = fit(
    warfarin_model,
    pop,
    init_params(warfarin_model),
    VEM(; nsamples = 30, optim_options = (; show_every = 20)),
)
[ 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     1.401185e+08     3.465709e+08
 * time: 2.2172927856445312e-5
    20     1.779390e+03     1.044671e+03
 * time: 19.64454698562622
    40     1.169506e+03     1.804089e+02
 * time: 34.248207092285156
    60     1.077426e+03     2.209827e+02
 * time: 49.936097145080566
    80     9.980722e+02     5.705262e+01
 * time: 66.4448812007904
   100     9.408512e+02     2.606154e+01
 * time: 85.47910809516907
   120     9.093491e+02     3.011094e+01
 * time: 107.69267416000366
   140     8.646411e+02     7.028322e+01
 * time: 129.8445611000061
   160     8.107487e+02     1.264954e+02
 * time: 152.70870804786682
   180     7.761239e+02     5.635161e+01
 * time: 176.03288507461548
   200     7.493039e+02     4.154542e+01
 * time: 201.0869140625
   220     7.427083e+02     9.517225e+00
 * time: 228.38062620162964
   240     7.337475e+02     2.446213e+01
 * time: 260.50611305236816
   260     7.282774e+02     8.923319e+00
 * time: 303.6597511768341
   280     7.253234e+02     1.270644e+01
 * time: 339.3359830379486
   300     7.236454e+02     3.827494e+00
 * time: 378.8135449886322
   320     7.228151e+02     1.233090e+00
 * time: 418.01746106147766
441.361976 seconds (630.95 M allocations: 129.036 GiB, 6.38% gc time)

FittedPumasModel

Dynamical system type:               Nonlinear ODE
Solver(s): (OrdinaryDiffEqVerner.Vern7,OrdinaryDiffEqRosenbrock.Rodas5P)

Number of subjects:                             32

Observation records:         Active        Missing
    conc:                       251             47
    pca:                        232             66
    Total:                      483            113

Number of parameters:      Constant      Optimized
                                  0             23

Estimation method:              VEM (EXPERIMENTAL)

-----------------------
            Estimate
-----------------------
pop_CL       0.13413
pop_V        8.0016
pop_tabs     0.38799
pop_lag      0.88115
pk_Ω₁,₁      0.067509
pk_Ω₂,₁      0.010518
pk_Ω₃,₁      0.030996
pk_Ω₂,₂      0.022964
pk_Ω₃,₂     -0.034492
pk_Ω₃,₃      1.0017
σ_prop       0.087658
σ_add        0.38644
pop_e0      96.233
pop_emax    -1.0717
pop_c50      1.5647
pop_tover   14.102
pd_Ω₁,₁      0.0031265
pd_Ω₂,₁     -0.004371
pd_Ω₃,₁      0.0025402
pd_Ω₂,₂      0.15374
pd_Ω₃,₂     -0.036933
pd_Ω₃,₃      0.011822
σ_fx         3.8995
-----------------------
loglikelihood(fpm_vem_nsamples, FOCE())
-1064.3370951585323

Notice that increasing the number of samples improved the quality of the ELBO however that does not guarantee an improvement in the FOCE-approximated log likelihood of the solution found. In this case, it improved it. Heuristically, setting nsamples to 15-30 seems to be a good starting point for most models, but the best value can vary based on the number of random effects in the model.

For multi-phase runs, the number of samples in phase p is: min(init_nsamples + (p - 1) * nsamples_increment, max_nsamples). More samples in later phases refine the estimate as the parameters approach their final values.

# Resampling with more phases and custom sample progression
fit(
    model,
    pop,
    params,
    VEM(;
        presample = false,
        nphases = 15,
        init_nsamples = 5,
        nsamples_increment = 5,
        max_nsamples = 50,
    ),
)

4.3 Optimizer Options (optim_options)

Fine-grained control over the optimizer is passed through the optim_options named tuple. The available options depend on which variant is used.

For presampling (L-BFGS):

Option Default Description
iterations 3000 Maximum iterations per phase
show_trace true Print convergence information
show_every 5 Print every N iterations
f_reltol 1e-5 Relative tolerance on the objective
g_tol 1e-3 Gradient norm tolerance 
fit(model, pop, params, VEM(; optim_options = (; iterations = 5000, show_every = 10)))

For resampling (Adam):

Option Default Description
iterations 20000 Maximum iterations per phase
show_trace true Print convergence information
show_every 10 Print every N iterations
subject_ratio 1.0 Fraction of subjects per iteration (mini-batching)

Setting subject_ratio < 1.0 enables mini-batching, where only a random subset of subjects is used per iteration. This can speed up each iteration for large datasets at the cost of noisier gradients:

fit(
    model,
    pop,
    params,
    VEM(; presample = false, optim_options = (; iterations = 30000, subject_ratio = 0.5)),
)

4.4 Prior on Population Parameters (prior)

By default (prior=true), VEM adds the log prior on the population parameters, if defined in the @param block, to the ELBO. This allows for maximum a-posteriori (MAP) estimation of the population parameters, after marginalizing the random effects. Setting prior=false ignores the prior contribution, which means the ELBO depends only on the likelihood:

fit(model, pop, params, VEM(; prior = false))

If the @param block has no prior distributions and only includes domains, this option is not used.

4.5 Base Distribution Standard Deviation (base_dist_std)

VEM’s variational family uses a reparameterization trick: samples are drawn from a base distribution and transformed to approximate the posterior. The base_dist_std parameter controls the standard deviation of this base distribution (default: 0.1):

fit(model, pop, params, VEM(; base_dist_std = 0.05))

Smaller values start the optimization with a tighter initial approximation to the posterior, which can help stability for sensitive models. Larger values allow broader initial exploration. The default of 0.1 works well for most models.

4.6 Posterior Covariance Structure (diagonal)

VEM approximates each subject’s random effect posterior with a multivariate normal distribution. The diagonal keyword controls the covariance structure:

  • diagonal=true (default): assumes independent random effects in the posterior. Faster and requires fewer variational parameters.
  • diagonal=false: captures correlations between random effects in the posterior. Slower but may be more accurate when random effects are correlated.
fpm_full_cov = fit(
    warfarin_model,
    pop,
    init_params(warfarin_model),
    VEM(; diagonal = false, optim_options = (; show_every = 20)),
)
[ 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     1.389956e+08     3.440554e+08
 * time: 2.3126602172851562e-5
    20     1.796437e+03     1.123577e+03
 * time: 13.070409059524536
    40     1.158877e+03     1.060873e+02
 * time: 23.59917712211609
    60     1.049436e+03     3.510573e+02
 * time: 34.05941820144653
    80     9.942815e+02     1.373982e+02
 * time: 48.534911155700684
   100     8.809855e+02     2.043411e+02
 * time: 65.42767715454102
   120     8.049193e+02     4.325779e+01
 * time: 80.48091411590576
   140     7.605175e+02     7.079459e+01
 * time: 97.34817409515381
   160     7.444436e+02     9.263546e+01
 * time: 112.41037011146545
   180     7.371362e+02     1.138209e+01
 * time: 126.38635015487671
   200     7.323646e+02     7.557712e+00
 * time: 140.4344391822815
   220     7.307324e+02     2.023366e+01
 * time: 154.487722158432
   240     7.287573e+02     1.147916e+01
 * time: 170.6592299938202
   260     7.267072e+02     1.428650e+01
 * time: 187.3003170490265
   280     7.249038e+02     2.100036e+00
 * time: 202.41970014572144
   300     7.230817e+02     4.044257e+00
 * time: 217.6405110359192
   320     7.213803e+02     1.051188e+00
 * time: 232.4256660938263
   340     7.200559e+02     3.963884e+00
 * time: 247.0711271762848
   360     7.193628e+02     6.851190e+00
 * time: 262.0750300884247
   380     7.184348e+02     5.934083e+00
 * time: 276.96819615364075
   400     7.170894e+02     1.430359e+01
 * time: 291.7761721611023
   420     7.160352e+02     4.594379e+00
 * time: 307.0602240562439
   440     7.150667e+02     2.376845e+00
 * time: 321.99050116539
   460     7.142590e+02     2.580791e+00
 * time: 337.0013391971588
   480     7.136455e+02     1.245333e+01
 * time: 352.19328117370605

FittedPumasModel

Dynamical system type:               Nonlinear ODE
Solver(s): (OrdinaryDiffEqVerner.Vern7,OrdinaryDiffEqRosenbrock.Rodas5P)

Number of subjects:                             32

Observation records:         Active        Missing
    conc:                       251             47
    pca:                        232             66
    Total:                      483            113

Number of parameters:      Constant      Optimized
                                  0             23

Estimation method:              VEM (EXPERIMENTAL)

------------------------
            Estimate
------------------------
pop_CL       0.13499
pop_V        8.0007
pop_tabs     0.65201
pop_lag      0.81868
pk_Ω₁,₁      0.06716
pk_Ω₂,₁      0.0066393
pk_Ω₃,₁      0.023315
pk_Ω₂,₂      0.009388
pk_Ω₃,₂      0.092654
pk_Ω₃,₃      0.94296
σ_prop       0.11638
σ_add        0.38621
pop_e0      96.513
pop_emax    -1.0774
pop_c50      1.5551
pop_tover   14.236
pd_Ω₁,₁      0.00037626
pd_Ω₂,₁      0.0070386
pd_Ω₃,₁     -0.0016503
pd_Ω₂,₂      0.13238
pd_Ω₃,₂     -0.033564
pd_Ω₃,₃      0.017918
σ_fx         4.2117
------------------------
loglikelihood(fpm_full_cov, FOCE())
-1088.6036716469694

Note that setting diagonal = false does not always lead to a better fit. It is generally recommended to start with diagonal = true and only switch to diagonal = false if there is reason to believe that the random effects are highly correlated.

4.7 Holding Parameters Constant (constantcoef)

You can hold specific parameters constant during estimation using the constantcoef keyword. This is a keyword to fit, not to VEM():

fit(model, pop, params, VEM(); constantcoef = (:σ_prop,))
Note

constantcoef is a keyword argument to fit, not to VEM(). This is consistent with how it works for FOCE and other algorithms.

4.8 Reproducibility (rng)

VEM uses random sampling internally, so results vary across runs by default. Set a specific random number generator for reproducible results:

rng = default_rng()
seed!(rng, 123)
fit(model, pop, params, VEM(; rng))
Note

Even with presampling, different seeds produce different results because the pre-drawn samples differ. For exact reproducibility across runs, always set the rng.

4.9 Threading (ensemblealg)

Control parallelism across subjects:

  • EnsembleThreads() (default): uses multiple threads for gradient computation.
  • EnsembleSerial(): single-threaded, useful for debugging.
fit(model, pop, params, VEM(; ensemblealg = EnsembleSerial()))

4.10 ODE Solver Options (diffeq_options)

Pass differential equation solver options through diffeq_options:

fit(model, pop, params, VEM(; diffeq_options = (; abstol = 1e-8, reltol = 1e-8)))

This is useful when the default ODE tolerances are insufficient for a stiff or sensitive model.

5 When to Use VEM vs FOCE/LaplaceI

Choose VEM when:

  • You have a response model that FOCE does not support and LaplaceI is too slow, e.g. truncated, censored or time-to-event data.
  • The model has many random effects and FOCE is too slow or struggles to converge, and you want a faster alternative, e.g. to initialize FOCE or as a standalone method.
  • You have a very large dataset and want a fast initial fit to use as a starting point for FOCE/LaplaceI, then stochastic VEM with mini-batching can help.

Choose FOCE when:

  • You have a model without too many random effects and the response model is supported by FOCE.
  • You need a well-established (non-experimental) fitting method, diagnostics and inference tools.

Choose LaplaceI when:

  • Your response model is not supported by FOCE.
  • You need a well-established (non-experimental) fitting method, diagnostics and inference tools.

6 Conclusion

In this tutorial, we introduced VEM as a flexible estimation method in Pumas. We covered its conceptual foundations (ELBO maximization, presampling vs. resampling, and posterior covariance structure), demonstrated its use with a warfarin PK model, showed how to compare results with FOCE, and provided a reference for the main VEM options.