Posterior (predictive) simulation for metafor rma models

sim.rma.uni.R

# Copyright (2021) Brenton M. Wiernik
# Licensed under GPL-3.0
#
sim.rma.uni <- function(object, n.sims = object$k, sim.type = c("coef", "sample"), newdata = NULL, control = NULL, ...) {
  sim.type <- match.arg(sim.type)
  if (!inherits(object, "rma.uni")) 
    stop(mstyle$stop("Argument 'object' must be an object of class \"rma.uni\"."))
  if (inherits(object, "robust.rma")) 
    stop(mstyle$stop("Method not available for objects of class \"robust.rma\"."))
  if (inherits(object, "rma.ls")) 
    stop(mstyle$stop("Method not available for objects of class \"rma.ls\"."))
  if (inherits(object, "rma.uni.selmodel")) 
    stop(mstyle$stop("Method not available for objects of class \"rma.uni.selmodel\"."))
  tau2.min <- ifelse(is.null(object$control$tau2.min), 
                     min(0, object$tau2), 
                     object$control$tau2.min)
  tau2.max <- ifelse(is.null(object$control$tau2.max), 
                     max(100, object$tau2 * 10, tau2.min * 10), 
                     object$control$tau2.max)
  con <- list(tol = .Machine$double.eps^0.25, 
              maxiter = 1000, 
              tau2.min = tau2.min, 
              tau2.max = tau2.max, 
              verbose = FALSE)
  con.pos <- pmatch(names(control), names(con))
  con[c(na.omit(con.pos))] <- control[!is.na(con.pos)]
  
  .sim.coef <- function (object, n.sims = 100, con = NULL) {
    k <- object$k
    p <- object$p
    yi <- object$yi
    vi <- object$vi
    X <- object$X
    Y <- cbind(yi)
    weights <- object$weights
    
    stXX <- metafor:::.invcalc(X = X, W = diag(k), k = k)
    P <- diag(k) - X %*% Matrix::tcrossprod(stXX, X)
    V <- diag(vi, nrow = k, ncol = k)
    PV <- P %*% V
    trPV <- metafor:::.tr(PV)
    
    summ <- summary(object)
    coef <- cbind(summ$b, summ$se)
    dimnames(coef)[[2]] <- c("coef.est", "coef.sd")
    tau2.hat <- m$tau2
    beta.hat <- coef[, 1, drop = FALSE]
    V.beta <- vcov(object) / object$s2w
    tau2 <- rep(NA, n.sims)
    s2w <- rep(NA, n.sims)
    beta <- array(NA, c(n.sims, p))
    dimnames(beta) <- list(NULL, rownames(beta.hat))
    
    for (s in 1:n.sims) {
      crit.q <- rchisq(1, k - p)
      res <- try(uniroot(metafor:::.QE.func, 
                         interval = c(con$tau2.min, con$tau2.max), 
                         tol = con$tol, 
                         maxiter = con$maxiter, 
                         Y = Y, vi = vi, X = X, k = k, 
                         objective = crit.q, 
                         digits = digits)$root, 
                 silent = TRUE)
      if (!inherits(res, "try-error")) {
        tau2[s] <- res
      } else {
        tau2[s] <- NA
      }
      
      s2w[s] <- (tau2[s] * (k - p) + trPV ) / (k - p)
      beta[s, ] <- MASS::mvrnorm(1, beta.hat, V.beta * s2w[s])
    }
    ans <- list(coef = beta, tau = sqrt(tau2))
    return(ans)
  }
  res.coef <- .sim.coef(object, n.sims = n.sims, con = con, ...)
  if (sim.type == "coef") {
    return(data.frame(res.coef$coef, tau = res.coef$tau))
  } else {
    if (is.null(newdata)) {
      yhat <- rowSums(object$X * res.coef$coef)
      yrep <- rnorm(n.sims, yhat, res.coef$tau)
      attributes(yrep) <- list(coef = data.frame(res.coef$coef, tau = res.coef$tau))
      return(yrep)
    } else {
      ypred <- newdata %*% res.coef$coef
      ypred <- rnorm(n.sims, yhat, res.coef$tau)
      attributes(ypred) <- list(coef = data.frame(res.coef$coef, tau = res.coef$tau))
      return(ypred)
    }
  }
}