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,
  par = NULL,
  ...
)

Arguments

obj: Optimised RTMB object generated by MakeADFun
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 precision matrix, dense or sparse. For models with random effects, this is typically the sparse matrix returned by sdreport(..., getJointPrecision = TRUE)$jointPrecision. If not provided, it is computed internally, using sdreport for models with random effects and the Hessian obj$he() otherwise. Q may have more rows than there are model parameters, see Details.
par: Optional mean of the samples. Defaults to obj$env$last.par.best. Must have length nrow(Q) if Q is supplied.
...: 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 devitions 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.

A user-supplied Q may be the joint precision of the model parameters and additional hyperparameters, such as log smoothing parameters, whose coordinates come after the model parameters. The mean par then needs to be given for all coordinates. Samples are drawn in the full space and the hyperparameters are dropped, so they are integrated out and the samples of the model parameters account for their uncertainty.

Examples

set.seed(123)

### Example 1: Without random effects ###

# simulate data
x <- rgumbel(100, location = 5, scale = 2)

# negative log-likelihood function
nll <- function(par) {
  getAll(par)
  x <- OBS(x) # mark x as the response
  scale <- exp(log_scale)
  ADREPORT(scale); REPORT(scale)
  -sum(dgumbel(x, loc, scale, log = TRUE))
}

# RTMB AD object
par <- list(loc = 5, log_scale = log(2))
obj <- MakeADFun(nll, par, silent = TRUE)

# model fitting using AD gradient
opt <- nlminb(obj$par, obj$fn, obj$gr)

# sample from asymptotic distribution of the MLE
samples <- mcreport(obj, report = TRUE)
#> Evaluating Hessian...
#> Computing reported quantities...
# parameters
mean(unlist(samples$loc)); sd(unlist(samples$loc))
#> [1] 5.00983
#> [1] 0.1991763
# reported transformed parameters
mean(unlist(samples$scale)); sd(unlist(samples$scale))
#> [1] 1.970023
#> [1] 0.1544266

# compare to delta method sd from sdreport
sdr <- sdreport(obj)
summary(sdr)
#>            Estimate Std. Error
#> loc       5.0015427 0.20659355
#> log_scale 0.6732893 0.07663174
#> scale     1.9606760 0.15025002



### Example 2: With random effects ###

# simulate random effects
id <- rep(1:10, each = 100) # 10 groups with 100 observations each
true_gamma <- rnorm(10, 0, 2) # random coefficients
true_probs <- plogis(1 + true_gamma) # linear predictor

# simulate Bernoulli data conditional on random coefficients
x <- rbinom(1000, 1, true_probs[id])

# negative log-likelihood function
nll <- function(par) {
  getAll(par)
  x <- OBS(x) # mark x as the response

  # random effect likelihood
  sd_gamma <- exp(log_sd_gamma)
  ADREPORT(sd_gamma); REPORT(sd_gamma)
  nll <- -sum(dnorm(gamma, 0, sd_gamma, log = TRUE))

  # data likelihood
  probs <- plogis(beta0 + gamma) # linear predictor
  ADREPORT(probs); REPORT(probs)
  nll <- nll - sum(dbinom(x, 1, probs[id], log = TRUE))
  return(nll)
}

# RTMB Laplace approximation
par <- list(beta0 = 1, log_sd_gamma = log(2), gamma = rep(0,10))
obj <- MakeADFun(nll, par, random = "gamma", silent = TRUE)

# model fitting using AD gradient
opt <- nlminb(obj$par, obj$fn, obj$gr)

# draw samples
samples <- mcreport(obj, report = TRUE)
#> Computing joint precision...
#> Computing reported quantities...

# run sdreport
sdr <- sdreport(obj)

# standard deviation using delta method
as.list(sdr, "Std", report = TRUE)$probs
#>  [1] 0.04426637 0.04799042 0.02115872 0.01608627 0.04376788 0.02817631
#>  [7] 0.04561685 0.02938416 0.03879701 0.04151316

# standard deviation from samples
apply(simplify2array(samples$probs), 1, sd)
#>  [1] 0.04405297 0.04911051 0.02291450 0.02046111 0.04423394 0.03029676
#>  [7] 0.04530734 0.03058053 0.03698335 0.04142802