Estimating Parameters using SAEM

PumasEMModels can optionally be fitted using SAEM. Here is an example PumasEMModel definition:

using Pumas
using PumasUtilities
using Random

covariate_saem_cov = @emmodel begin
    @random begin
        CL ~ 1 + logwt | LogNormal
        v ~ 1 | LogNormal
    end
    @covariance 2
    @covariates wt
    @pre begin
        Vc = wt * v
    end
    @dynamics Central1
    @post begin
        cp = Central / Vc
    end
    @error begin
        dv ~ ProportionalNormal(cp)
    end
end
PumasEMModel
  Parameters with random effects:
    CL ~ (1, :logwt) | LogNormal
    v ~ (1,) | LogNormal
  Covariates: wt
  Pre-dynamical variables: Vc
  Dynamical system variables: Central
  Post-dynamical variables: cp

See the documentation on the @emmodel macro interface for an explanation of the syntax. To fit this model, we'll simulate data using simobs.

sim_params_covariate_cov = (;
    CL = (4.0, 0.75),
    v = 70.0,
    Ω = (Pumas.@SMatrix([0.1 0.05; 0.05 0.1]),),
    σ = ((0.2,),),
)

obstimes = 0.0:100.0

dose = DosageRegimen(1_000; addl = 2, ii = 24)

function choose_covariates()
    wt = (55 + 25rand()) / 70
    return (; wt, logwt = log(wt))
end

pop = [Subject(; id = i, events = dose, covariates = choose_covariates()) for i = 1:72]
sims = simobs(
    covariate_saem_cov,
    pop,
    sim_params_covariate_cov;
    obstimes,
    ensemblealg = EnsembleSerial(),
)
reread_df = DataFrame(sims);

pop_covariate = read_pumas(reread_df; observations = [:dv], covariates = [:logwt, :wt])
Population
  Subjects: 72
  Covariates: logwt, wt
  Observations: dv

When specifying inits, it is not necessary to specify any variance parameters for the random effects or error models.

init_covariate = (; CL = (2.0, 2 / 3), v = 50.0)
(CL = (2.0, 0.6666666666666666),
 v = 50.0,)

If unspecified, they will be initialized to 1.0 or the identity matrix for covariances. SAEM's stochastic exploration phase is most effective when these variance parameters are much larger than the true variances. Thus, if the true variances are believed to be around 1.0 or larger, it is recommended to specify larger initial values manually.

It is possible to pass a vector of random number generates as an argument to fit. If so, the fit will use one thread per RNG. By specifying the seeds of each RNG, SAEM() can be fully reproducible:

rngv = [MersenneTwister(1941964947i + 1) for i ∈ 1:Threads.nthreads()];

fit_covariate_cov1 = fit(
    covariate_saem_cov,
    pop_covariate,
    init_covariate,
    SAEM();
    ensemblealg = EnsembleThreads(),
    rng = rngv,
)
rngv = [MersenneTwister(1941964947i + 1) for i ∈ 1:Threads.nthreads()];

fit_covariate_cov2 = fit(
    covariate_saem_cov,
    pop_covariate,
    init_covariate,
    SAEM();
    ensemblealg = EnsembleThreads(),
    rng = rngv,
)
coef(fit_covariate_cov1) == coef(fit_covariate_cov2) # true
true

One can also specify the number of iterations for each of the three phases (rapid exploration, convergence, smoothing):

fit_covariate_cov =
    fit(covariate_saem_cov, pop_covariate, init_covariate, SAEM(; iters = (1000, 500, 500)))
FittedPumasEMModel

Dynamical system type:                 Closed form

Number of subjects:                             72

Observation records:         Active        Missing
    dv:                        7272              0
    Total:                     7272              0

Number of parameters:      Constant      Optimized
                                  0              7

Likelihood approximation:                     SAEM

Log-likelihood value:                   -12563.264

-------------------
         Estimate
-------------------
CL₁       3.5701
CL₂       0.51811
v        70.444
Ω₁,₁      0.1203
Ω₂,₁      0.043111
Ω₂,₂      0.088156
σ         0.20033
-------------------

The results of a fit can be analyzed normally:

infer_cov = infer(fit_covariate_cov)
Asymptotic inference results using sandwich estimator

Dynamical system type:                 Closed form

Number of subjects:                             72

Observation records:         Active        Missing
    dv:                        7272              0
    Total:                     7272              0

Number of parameters:      Constant      Optimized
                                  0              7

Likelihood approximation:                     SAEM

Log-likelihood value:                   -12563.264

---------------------------------------------------------
           Estimate   SE          95.0% C.I.
---------------------------------------------------------
CL_base     3.5701    0.14974     [  3.2766  ;  3.8635 ]
CL_logwt    0.51811   0.32111     [ -0.11125 ;  1.1475 ]
v_base     70.444     2.4839      [ 65.576   ; 75.312  ]
ω_1₁,₁      0.34685   0.028935    [  0.29014 ;  0.40356]
ω_1₂,₁      0.12429   0.031662    [  0.062237;  0.18635]
ω_1₂,₂      0.26964   0.019835    [  0.23077 ;  0.30852]
σ_0         0.20033   0.0014327   [  0.19752 ;  0.20314]
---------------------------------------------------------
coeftable(infer_cov)
7×7 DataFrame
Rowparameterconstantestimateserelative_seci_lowerci_upper
StringBoolFloat64Float64Float64Float64Float64
1CL_basefalse3.570050.1497360.04194223.276573.86353
2CL_logwtfalse0.5181090.3211070.619767-0.1112491.14747
3v_basefalse70.44412.483890.035260565.575875.3124
4ω_1₁,₁false0.3468490.02893510.08342270.2901370.403561
5ω_1₂,₁false0.1242940.03166230.2547370.0622370.186351
6ω_1₂,₂false0.2696430.01983510.07356060.2307670.308519
7σ_0false0.200330.001432720.007151820.1975220.203138
inspect_cov = inspect(fit_covariate_cov)
[ Info: Calculating predictions.
[ Info: Calculating weighted residuals.
[ Info: Calculating empirical bayes.
[ Info: Evaluating dose control parameters.
[ Info: Evaluating individual parameters.
[ Info: Done.