Skip to content

Instantly share code, notes, and snippets.

@bwiernik
Created February 16, 2021 14:14
Show Gist options
  • Save bwiernik/be04c8817046ea0ffea9e9838928ea96 to your computer and use it in GitHub Desktop.
Save bwiernik/be04c8817046ea0ffea9e9838928ea96 to your computer and use it in GitHub Desktop.
Posterior (predictive) simulation for metafor rma models
# 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)
}
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment