Estimating Parameters of Pumas Models

fit(
  model::PumasModel,
  population::Population,
  param::NamedTuple,
  approx::Union{LikelihoodApproximation, MAP};
  optim_alg = Optim.BFGS(linesearch = Optim.LineSearches.BackTracking(), initial_stepnorm = 1.0),
  optim_options::NamedTuple = NamedTuple(),
  optimize_fn = nothing,
  constantcoef::Tuple = (),
  ensemblealg::SciMLBase.EnsembleAlgorithm = EnsembleThreads(),
  checkidentification = true,
  diffeq_options = NamedTuple(),
  verbose = true,
  ignore_numerical_error = true,
  )

Fit the Pumas model model to the dataset population with starting values param using the estimation method approx. Currently supported values for the approx argument are FO, FOCE, LaplaceI, NaivePooled, and BayesMCMC. See the online documentation for more details about the different methods.

To control the optimization procedure for finding population parameters (fixed effects), use the optim_alg and optim_options keyword arguments. In previous versions of Pumas the argument optimize_fn was used, but is now discouraged and will be removed in a later version of Pumas. These options control the optimization all approx methods except BayesMCMC. The default optimization function uses the quasi-Newton routine BFGS method from the Optim package. It can be changed by setting the optim_alg to an algorithm implemented in Optim.jl as long as it does not use second order derivatives. Optimization specific options can be passed in using the optim_options keyword and has to be a NamedTuple with names and values that match the Optim.Options type. For example, the optimization trace can be disabled and the algorithm can be changed to L-BFGS by passing optim_alg=Optim.LBFGS(), optim_options = ;(show_trace=false) to fit. See Optim for more defails.

It is possible to fix one or more parameters of the fit by passing a Tuple of Symbols as the constantcoef argument with elements corresponding to the names of the fixed parameters, e.g. constantcoef=(:σ,).

When models include an @random block and fitting with NaivePooled is requested, it is required that the user sets all variability parameters to zero with constantcoef such that these can be ignored in the optimization, e.g. constantcoef = (:Ω,) while overwriting the corresponding values in param with (; init_params..., Ω = zeros(3, 3)).

Parallelization of the optimization is supported for most estimation methods via the ensemble interface of DifferentialEquations.jl. Currently, the only supported options are:

• EnsembleThreads(): the default. Accelerate fits by using multiple threads.
• EnsembleSerial(): fit using a single thread.
• EnsembleDistributed(): fit by using multiple worker processes.

The fit function will check if any gradients and throw an exception if any of the elements are exactly zero unless checkidentification is set to false.

Further keyword arguments can be passed via the diffeq_options argument. This allows for passing arguments to the differential equations solver such as alg, abstol, and reltol. The default values for these are AutoVern7(Rodas5P(autodiff=true)), 1e-12, and 1e-8 respectively. See the DifferentialEquations.jl documentation for more details.

The keyword verbose controls if info statements about initial evaluations of the loglikelihood function and gradient should be printed or not. Defaults to true.

Since the numerical optimization routine can sometimes visit extreme regions of the parameter space during it's exploration we automatically handle the situation where the input parameters result in errors when solving the dynamical system. This allows the algorithm to recover and continue in many situations that would have otherwise stalled early. Sometimes, it is useful to turn this error handling off when debugging a model fit, and this can be done by setting ignore_numerical_error = false.

fit(
  model::PumasEMModel,
  population::Population,
  param::NamedTuple,
  approx::Union{SAEM,LaplaceI};
  ensemblealg = EnsembleThreads(),
  rng = default_rng()
)

Fit the PumasEMModel to the dataset population using the initial values specified in param. The supported methods for the approx argument are currently SAEM and LaplaceI. See the online documentation for more details about these two methods.

When fitting with SAEM the variance terms of the random effects (Ω) and the dispersion parameters for the error model (σ) are initialized to the identity matrix or 1 as appropriate. They may also be specified. With SAEM, it is recommended to choose an init larger than the true values to facilitate exploration of the parameter space and avoid getting trapped in local optima early. Currently LaplaceI fits with a PumasEMModel require Ω to be diagonal.

Options for ensemblealg are EnsembleSerial() and EnsembleThreads() (the default). The fit will be parallel if ensemblealg == EnsembleThreads(). SAEM imposes the additional requirement for threaded fits that either the rng used supports Future.randjump or that rng be a vector of rngs, one per thread. The default_rng() supports randjump().

Distributions.fit(
  model::PumasModel,
  data::Population,
  param::NamedTuple,
  alg::BayesMCMC,
)

Samples from the posterior distribution of the model's population and subject-specific parameters given the observations in data using the Hamiltonian Monte Carlo (HMC) based No-U-Turn Sampler (NUTS) algorithm. The MCMC sampler is initialized starting from param for the population parameters. The subject-specific parameters are also apprioriately initialized. The Bayesian inference algorithm's options are passed in the alg struct. Run ?BayesMCMC for more details on the options that can be set.