R/mcreport.R

In LaMa: Fast Numerical Maximum Likelihood Estimation for Latent Markov Models

#' Sample parameters from approximate Gaussian posterior distribution
#' 
#' Efficient Monte Carlo sampling of parameters (and \code{REPORT}ed quantities) from the approximate posterior of an \code{(R)TMB} model.
#' See \code{\link[TMB]{sdreport}} for details on posterior variance-covariance in random effects models.
#'
#' @details
#' \strong{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 \code{\link[TMB]{sdreport}} is particularly useful for multidimensional random effects (e.g. for locations \eqn{x} and \eqn{y})
#' where pointwise confidence intervals (e.g. along a path) based on standard deviations are not possible.
#'
#' @param obj 
#' Optimised \code{RTMB} object generated by \code{\link[RTMB]{MakeADFun}}
#' 
#' @param nSamples 
#' Number of samples to draw
#' 
#' @param include_random_pars 
#' Logical; Should random parameters be included in the output?
#' 
#' @param report 
#' Logical; Should \code{REPORT}ed quantities be sampled as well? 
#' Defaults to \code{FALSE} because this may be slow depending on your model.
#' 
#' @param Q 
#' Optional precalculated sparse precision matrix returned by \code{sdreport(..., getJointPrecision = TRUE)$jointPrecision}.
#' If not provided, computed internally using \code{sdreport}. Only used for models with random effects.
#' 
#' @param ... For internal use only
#'
#' @returns A list structured like the original parameter list used in the \link[RTMB]{MakeADFun} call (potentially including additional \code{REPORT}ed quantities). Each entry is a list with \code{nSamples} entries.
#' 
#' @export
#' 
#' @importFrom RTMBdist mcreport
#'
#' @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)
MCreport <- function(obj,
                     nSamples = 1000,
                     include_random_pars = TRUE,
                     report = FALSE,
                     Q = NULL,
                     ...) {
  
  RTMBdist::mcreport(obj = obj,
                     nSamples = nSamples,
                     include_random_pars = include_random_pars,
                     report = report,
                     Q = Q,
                     ...)
}