Estimating Parameters of Pumas Models

Estimating Parameters of Pumas Models

Pumas can use the observational data of a Subject or Population to estimate the parameters of an NLME model. This is done by two classes of methods. Maximum likelihood methods find the parameters such that the observational data has the highest probability of occurring according to the chosen error distributions. Bayesian methods find a posterior probability distribution for the parameters to describe the chance that a parameter has a given value given the data. The following section describes how to fit an NLME model in Pumas via the two methods.

Defining Data for Estimation

Estimation is done by looking at the likelihood of likewise names. Thus, for example if subject.observations is a NamedTuple with names dv and resp, the likelihood calculation will be between values from derived named dv and resp. If dv is a scalar in the observation data, then dv from derived should also be a scalar. Likewise, if dv is an array like a time series, then dv should be a size-matching time series when returned from derived. Note that likelihoods are calculated between the probability distribution from derived and the matching observation from subject.observations. If no likelihood is associated with a derived value, then the value has an implicit standard normal interpretation, which amounts to having the L2 Euclidian distance taken as the likelihood during fitting procedures.

Note: Currently only BayesMCMC and LaplaceI support multiple output series. all other likelihood approximations must have the derived and observation data under the name dv

Note: Currently the estimation procedures require that there only exists a single random effect (vector), and this vector must be named η (\eta)

Maximum Likelihood Estimation

Maximum Likelihood Estimation (MLE) is performed using the fit function. This function's signature is:

Distributions.fit(m::PumasModel,
                  data::Population,
                  param::NamedTuple,
                  approx::LikelihoodApproximation,
                  args...;
                  optimize_fn = DEFAULT_OPTIMIZE_FN,
                  kwargs...)

The arguments are:

function DEFAULT_OPTIMIZE_FN(cost,p)
  Optim.optimize(cost,p,BFGS(linesearch=Optim.LineSearches.BackTracking()),
                 Optim.Options(show_trace=verbose, # Print progress
                               store_trace=true,
                               extended_trace=true,
                               g_tol=1e-3),
                  autodiff=:finite)
end

The return type of fit is a FittedPumasModel.

Marginal Likelihood Approximations

The following choices are available for the likelihood approximations:

FittedPumasModel

The relevant fields of a FittedPumasModel are:

Additionally, the following functions help interpret the fit:

Additionally the function:

predict(fpm::FittedPumasModel, approx=fpm.approx;
        nsim=nothing, timegrid=false, newdata=false, useEBEs=true)

Returns a FittedPumasPrediction which contains the solution of all population diagnostics in the field population and all individual diagnostics in the field individual. For more information on the diagnostics, please see the Simulation and Estimation Diagnostics

Bayesian Estimation

Bayesian parameter estimation is performed by using the fit function as follows:

fit(model::PumasModel, data::Population, ::BayesMCMC,
                       args...; nsamples=5000, kwargs...)

The arguments are:

The result is a BayesMCMCResults type. Notice that initial parameter values are not utilized in Bayesian estimation.

BayesMCMCResults

The MCMC chain is stored in the chain field of the returned BayesMCMCResults. The following functions help with querying common results on the Bayesian posterior: