R/mcreport.R

In RTMBdist: Distributions Compatible with Automatic Differentiation by 'RTMB'

# helper functions to mcreport()
transpose_samples <- function(x) {
  nm <- unique(unlist(lapply(x, names)))
  out <- stats::setNames(vector("list", length(nm)), nm)
  for(k in nm) {
    out[[k]] <- lapply(x, function(el) el[[k]])
  }
  out
}
append_report_suffix <- function(report_names, par_names, suffix = ".report") {
  ind <- report_names %in% par_names
  report_names[ind] <- paste0(report_names[ind], suffix)
  report_names
}


#' 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 \code{\link[RTMB]{MakeADFun}} call (potentially including additional \code{REPORT}ed quantities). Each entry is a list with \code{nSamples} entries.
#'
#' @export
#'
#' @importFrom stats optimHess
#' @importFrom Matrix Matrix
#' @importFrom RTMB sdreport
#'
#' @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)
#' # parameters
#' mean(unlist(samples$loc)); sd(unlist(samples$loc))
#' # reported transformed parameters
#' mean(unlist(samples$scale)); sd(unlist(samples$scale))
#'
#' # compare to delta method sd from sdreport
#' sdr <- sdreport(obj)
#' summary(sdr)
#'
#'
#'
#' ### 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)
#'
#' # run sdreport
#' sdr <- sdreport(obj)
#'
#' # standard deviation using delta method
#' as.list(sdr, "Std", report = TRUE)$probs
#'
#' # standard deviation from samples
#' apply(simplify2array(samples$probs), 1, sd)
mcreport <- function(obj,
                     nSamples = 1000,
                     include_random_pars = TRUE,
                     report = FALSE,
                     Q = NULL,
                     ...) {
  if(nSamples < 1) stop("At least one sample needs to be drawn")

  # Function to map parameter vector to list
  relist_par <- obj$env$parList

  # MLE or c(MLE, predicted random effects)
  p_hat <- obj$env$last.par.best

  # Are random effects present?
  random_ind <- obj$env$random # index of random effects
  random <- !is.null(random_ind)
  random_names <- unique(names(p_hat[random_ind]))

  # Quantities to exclude when REPORT sampling
  excl <- c("allprobs", "type", "trackID")
  # Should allprobs be included?
  dots <- list(...)
  include_allprobs <- isTRUE(dots$include_allprobs)
  if(include_allprobs) excl <- setdiff(excl, "allprobs")

  isRTMB <- ifelse(obj$env$DLL == "RTMB", TRUE, FALSE)

  if(random) {
    if(is.null(Q)) {
      # model with random effects -> use joint precision approx from TMB
      message("Computing joint precision...")
      if(isRTMB) {
        Q <- RTMB::sdreport(obj, getJointPrecision = TRUE, skip.delta.method = TRUE,
                            getReportCovariance = FALSE)$jointPrecision
      } else {
        Q <- TMB::sdreport(obj, getJointPrecision = TRUE, skip.delta.method = TRUE,
                           getReportCovariance = FALSE)$jointPrecision
      }

    } else {
      # check if Q is sparse matrix
      if(!inherits(Q, "sparseMatrix")) {
        stop("Q does not seem to be generated by sdreport(..., getJointPrecision = TRUE)")
      }
      # check if Q has correct dimensions
      if(any(dim(Q) != length(p_hat))) {
        stop("Q has incorrect dimensions. Expected ", length(p_hat), " but got ", nrow(Q))
      }
    }
  } else {
    # model without random effects -> use Hessian
    if(!is.null(Q)) {
      warning("Argument Q is ignored for models without random effects")
    }
    message("Evaluating Hessian...")
    Hessian <- obj$he(p_hat)
    Q <- Matrix(Hessian, sparse = TRUE)
  }

  # Sampling parameter vectors
  samples <- rgmrf(nSamples, p_hat, Q)

  # Relisting samples
  par_samples <- lapply(seq_len(nSamples), function(i) {
    p <- samples[i, ]
    relist_par(par = p)
  })
  # nicer shape: named list outside
  par_samples <- transpose_samples(par_samples)

  if(!include_random_pars && random) {
    par_samples <- par_samples[!names(par_samples) %in% random_names]
  }

  # if report = TRUE, call obj$report() on each sampled parameter vector
  if(report) {
    message("Computing reported quantities...")
    report_samples <- lapply(seq_len(nSamples), function(i) {
      p <- samples[i, ]
      r <- obj$report(par = p)
      r[!names(r) %in% excl] # exclude large, unnecessary quantities
    })
    # nicer shape: named list outside
    report_samples <- transpose_samples(report_samples)
    # adding .report to names when par and report share names to avoid collision when merging
    names(report_samples) <- append_report_suffix(names(report_samples),names(par_samples))
    # merging
    par_samples <- c(par_samples, report_samples)
  }

  return(par_samples)
}