Estimating Parameters using SAEM
PumasEMModel
s 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 init
s, 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
Likelihood approximation: SAEM
Dynamical system type: Closed form
Log-likelihood value: -11437.538
Number of subjects: 72
Number of parameters: Fixed Optimized
0 7
Observation records: Active Missing
dv: 7272 0
Total: 7272 0
-------------------
Estimate
-------------------
CL₁ 4.126
CL₂ 1.2434
v 72.166
Ω₁,₁ 0.10776
Ω₂,₁ 0.039826
Ω₂,₂ 0.093114
σ 0.20126
-------------------
The results of a fit can be analyzed normally:
infer_cov = infer(fit_covariate_cov)
Asymptotic inference results using sandwich estimator
Likelihood approximation: SAEM
Dynamical system type: Closed form
Log-likelihood value: -11437.538
Number of subjects: 72
Number of parameters: Fixed Optimized
0 7
Observation records: Active Missing
dv: 7272 0
Total: 7272 0
-------------------------------------------------------------------
Estimate SE 95.0% C.I.
-------------------------------------------------------------------
CL_base 4.126 0.16585 [ 3.801 ; 4.4511 ]
CL_logwt 1.2434 0.30879 [ 0.6382 ; 1.8486 ]
v_base 72.166 2.6111 [67.048 ; 77.283 ]
ω_1₁,₁ 0.32826 0.02201 [ 0.28512 ; 0.3714 ]
ω_1₂,₁ 0.12132 0.031493 [ 0.059598; 0.18305]
ω_1₂,₂ 0.27999 0.031657 [ 0.21794 ; 0.34204]
σ_0 0.20126 0.0017092 [ 0.19791 ; 0.20461]
-------------------------------------------------------------------
coeftable(infer_cov)
Row | parameter | estimate | se | relative_se | ci_lower | ci_upper |
---|---|---|---|---|---|---|
String | Float64 | Float64 | Float64 | Float64 | Float64 | |
1 | CL_base | 4.126 | 0.165845 | 0.0401952 | 3.80095 | 4.45106 |
2 | CL_logwt | 1.24341 | 0.308789 | 0.24834 | 0.638195 | 1.84863 |
3 | v_base | 72.1658 | 2.61111 | 0.0361822 | 67.0481 | 77.2835 |
4 | ω_1₁,₁ | 0.328261 | 0.0220096 | 0.0670492 | 0.285123 | 0.371399 |
5 | ω_1₂,₁ | 0.121324 | 0.0314935 | 0.259582 | 0.0595977 | 0.18305 |
6 | ω_1₂,₂ | 0.27999 | 0.0316571 | 0.113065 | 0.217943 | 0.342037 |
7 | σ_0 | 0.201257 | 0.00170923 | 0.00849278 | 0.197907 | 0.204607 |
inspect_cov = inspect(fit_covariate_cov)
[ Info: Calculating predictions.
[ Info: Calculating weighted residuals.
[ Info: Calculating empirical bayes.
[ Info: Evaluating individual parameters.
[ Info: Evaluating dose control parameters.
[ Info: Done.