Sample parameters from approximate Gaussian posterior distribution

Efficient Monte Carlo sampling of parameters (and REPORTed quantities) from the approximate posterior of an (R)TMB model. See sdreport for details on posterior variance-covariance in random effects models.

Usage

MCreport(
  obj,
  nSamples = 1000,
  include_random_pars = TRUE,
  report = FALSE,
  Q = NULL,
  joint = FALSE,
  ...
)

Arguments

obj: Optimised RTMB object generated by MakeADFun, or a model object from lama or qreml
nSamples: Number of samples to draw
include_random_pars: Logical; Should random parameters be included in the output?
report: Logical; Should REPORTed quantities be sampled as well? Defaults to FALSE because this may be slow depending on your model.
Q: Optional precalculated precision matrix, e.g. returned by sdreport(..., getJointPrecision = TRUE)$jointPrecision for models with random effects, or by joint_uncertainty() of a model fitted by qreml. If not provided, computed internally using sdreport for models with random effects and the Hessian otherwise.
joint: Logical; only for models fitted by qreml. If TRUE, the samples account for the uncertainty in the estimated penalty strengths, by sampling from the joint precision of the parameters and the log penalty strengths returned by joint_uncertainty() of the model, and dropping the latter. This requires 2L additional model fits for L penalty strengths, see qreml. If the joint precision has already been computed, it can instead be passed as Q.
...: For internal use only

Value

A list structured like the original parameter list used in the MakeADFun call (potentially including additional REPORTed quantities). Each entry is a list with nSamples entries.

Details

Caution: This has nothing to do with Bayesian posterior sampling. It is simple Monte Carlo sampling from a Gaussian distribution with mean at the MLE and covariance given by the inverse of the Hessian (for fixed effects models) or the joint precision matrix (for random effects models). This is a common approach to get an approximate idea of parameter uncertainty around the MLE, but it relies on large-sample asymptotics.

Sampling random effects from their posterior as compared to calculating marginal standard deviations like sdreport is particularly useful for multidimensional random effects (e.g. for locations $x$ and $y$) where pointwise confidence intervals (e.g. along a path) based on standard deviations are not possible.

For models fitted by qreml, the samples are centred at the final estimate. By default, they are conditional on the estimated penalty strengths. With joint = TRUE, or a Q from joint_uncertainty(), they also account for their uncertainty.

Examples

step <- trex$step[1:2000] # subsetting trex data
N <- 2                    # 2 states

# custom likelihood
nll <- function(par) {
  getAll(par)
  Gamma <- tpm(eta)
  delta <- stationary(Gamma)
  mu <- exp(log_mu); REPORT(mu)
  sigma <- exp(log_sigma); REPORT(sigma)
  allprobs <- matrix(1, length(step), N)
  for(j in 1:N) allprobs[,j] <- dgamma2(step, mu[j], sigma[j])
  -forward(delta, Gamma, allprobs)
}

# initial parameters in named list
par0 <- list(eta = rep(-1.5,2), 
             log_mu = log(c(0.3, 2.5)), 
             log_sigma = log(c(0.3, 1.5)))
    
# constructing AD object        
obj <- MakeADFun(nll, par0, silent = TRUE)

# optimising
opt <- nlminb(obj$par, obj$fn, obj$gr)

# sampling from distribution of the MLE
samples <- MCreport(obj, nSamples = 10, report = TRUE)
#> Evaluating Hessian...
#> Computing reported quantities...