Variational Expectation Maximization (VEM)

Authors

Mohamed Tarek

Vijay Ivaturi

Experimental 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
Same 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.

Random 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: 5.507469177246094e-5
    20     1.782428e+03     1.102616e+03
 * time: 19.335196018218994
    40     1.158299e+03     1.183852e+02
 * time: 32.52602815628052
    60     1.082972e+03     5.027322e+02
 * time: 44.819297075271606
    80     1.033477e+03     1.661968e+02
 * time: 59.093116998672485
   100     9.683433e+02     9.345296e+01
 * time: 72.76539421081543
   120     9.206464e+02     2.761997e+01
 * time: 90.99577903747559
   140     8.862514e+02     9.321403e+01
 * time: 107.25813698768616
   160     8.339402e+02     2.689859e+01
 * time: 122.82978820800781
   180     7.830751e+02     5.439420e+00
 * time: 141.59120416641235
   200     7.463957e+02     4.798360e+01
 * time: 159.35408806800842
   220     7.338703e+02     2.505895e+01
 * time: 179.7639501094818
   240     7.269872e+02     6.267607e+00
 * time: 200.60550117492676
   260     7.245847e+02     3.532142e+00
 * time: 229.56864309310913
   280     7.237545e+02     6.477516e+00
 * time: 248.1946120262146
271.837089 seconds (285.32 M allocations: 57.372 GiB, 5.95% gc time, 1.54% 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
Log-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
Weighted 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])

wres_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
VEM 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, 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, 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, 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/lK187/src/estimation/vem.jl:2594
 [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, 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, 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, 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/lK187/src/estimation/vem.jl:2587
 [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: 9.393692016601562e-5
    20     1.069323e+03     3.312721e+01
 * time: 97.22760677337646
    40     1.061185e+03     3.274055e+01
 * time: 142.4888298511505
167.069051 seconds (120.41 M allocations: 30.402 GiB, 3.96% gc time, 0.23% 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   ]
-----------------------------------------------------------------
True 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.

VEM 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: 0.0001289844512939453
    20     4.290829e+03     1.936591e+03
 * time: 40.58150291442871
    40     1.622225e+03     5.112982e+02
 * time: 110.57725405693054
    60     1.380720e+03     3.611768e+02
 * time: 149.61574792861938
    80     1.271772e+03     1.229157e+02
 * time: 201.49630498886108
   100     1.264818e+03     3.745549e+01
 * time: 251.8838930130005
   120     1.255091e+03     9.895804e+01
 * time: 289.92402601242065
   140     1.195682e+03     4.017472e+01
 * time: 326.14822483062744
   160     1.189284e+03     5.575148e+01
 * time: 361.68308901786804
   180     1.183513e+03     4.676253e+01
 * time: 401.9361319541931
   200     1.169818e+03     1.753033e+01
 * time: 476.3280129432678
   220     1.141433e+03     8.662493e+01
 * time: 564.7921049594879
   240     1.125162e+03     3.124826e+01
 * time: 661.1453309059143
   260     1.082971e+03     9.731992e+00
 * time: 737.7062220573425
   280     1.080982e+03     4.115973e+00
 * time: 792.7242250442505
   300     1.079409e+03     4.518367e+00
 * time: 848.3918039798737
   320     1.071216e+03     7.824287e+00
 * time: 931.0139198303223
   340     1.064598e+03     6.447437e+00
 * time: 990.1208908557892
   360     1.060925e+03     2.996939e+01
 * time: 1048.3557970523834
1072.991300 seconds (363.16 M allocations: 139.020 GiB, 3.47% 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
1703.812826 seconds (1.95 G allocations: 424.332 GiB, 6.15% gc time, 0.35% 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.

Reading 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: 4.100799560546875e-5
    20     1.779390e+03     1.044671e+03
 * time: 24.929980039596558
    40     1.169506e+03     1.804089e+02
 * time: 43.72418785095215
    60     1.077426e+03     2.209827e+02
 * time: 64.7934000492096
    80     9.980722e+02     5.705262e+01
 * time: 88.03809404373169
   100     9.408512e+02     2.606154e+01
 * time: 114.60482501983643
   120     9.093491e+02     3.011094e+01
 * time: 147.11614084243774
   140     8.646411e+02     7.028322e+01
 * time: 182.01240801811218
   160     8.107487e+02     1.264954e+02
 * time: 214.99422597885132
   180     7.761239e+02     5.635161e+01
 * time: 245.26629304885864
   200     7.493039e+02     4.154542e+01
 * time: 277.0907199382782
   220     7.427083e+02     9.517225e+00
 * time: 311.88948798179626
   240     7.337475e+02     2.446213e+01
 * time: 348.2466230392456
   260     7.282774e+02     8.923319e+00
 * time: 387.7636330127716
   280     7.253234e+02     1.270644e+01
 * time: 426.74086689949036
   300     7.236454e+02     3.827494e+00
 * time: 465.3489818572998
   320     7.228151e+02     1.233090e+00
 * time: 506.0413439273834
529.183229 seconds (630.95 M allocations: 129.061 GiB, 6.13% 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.6941299438476562e-5
    20     1.796437e+03     1.123577e+03
 * time: 13.35403299331665
    40     1.158877e+03     1.060873e+02
 * time: 23.03213095664978
    60     1.049436e+03     3.510573e+02
 * time: 33.88619112968445
    80     9.942815e+02     1.373982e+02
 * time: 46.71410894393921
   100     8.809855e+02     2.043411e+02
 * time: 59.53821396827698
   120     8.049193e+02     4.325779e+01
 * time: 73.61116695404053
   140     7.605175e+02     7.079459e+01
 * time: 88.43140697479248
   160     7.444436e+02     9.263546e+01
 * time: 104.4025809764862
   180     7.371362e+02     1.138209e+01
 * time: 119.92857098579407
   200     7.323646e+02     7.557712e+00
 * time: 135.90189409255981
   220     7.307324e+02     2.023366e+01
 * time: 152.31331205368042
   240     7.287573e+02     1.147916e+01
 * time: 167.19561910629272
   260     7.267072e+02     1.428650e+01
 * time: 182.9012279510498
   280     7.249038e+02     2.100036e+00
 * time: 199.25247406959534
   300     7.230817e+02     4.044257e+00
 * time: 214.75546503067017
   320     7.213803e+02     1.051188e+00
 * time: 230.2975459098816
   340     7.200559e+02     3.963884e+00
 * time: 245.46332097053528
   360     7.193628e+02     6.851190e+00
 * time: 260.6239900588989
   380     7.184348e+02     5.934083e+00
 * time: 276.0297830104828
   400     7.170894e+02     1.430359e+01
 * time: 292.14627504348755
   420     7.160352e+02     4.594379e+00
 * time: 307.47974705696106
   440     7.150667e+02     2.376845e+00
 * time: 324.04764199256897
   460     7.142590e+02     2.580791e+00
 * time: 340.25500106811523
   480     7.136455e+02     1.245333e+01
 * time: 355.6853630542755

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.