Calculate sandwich covariance estimators for multi-variate meta-analysis models

metafor-sandwich.R

require(Formula)
require(metafor)
require(sandwich)
require(zoo)
require(lmtest)

#-----------------------------------------------
# Functions for making sandwich standard errors
#-----------------------------------------------

findCluster <- function(obj) {
  if (is.null(obj$cluster)) {
    if (obj$withS) {
      r <- which.min(obj$s.nlevels)
      cluster <- obj$mf.r[[r]][[obj$s.names[r]]]
    } else if (obj$withG) {
      cluster <- obj$mf.r[[1]][[obj$g.names[2]]]
    } else {
      stop("No clustering variable specified.")
    }
  } else {
    cluster <- obj$cluster
  }
  cluster
}

bread.rma.mv <- function(obj) {
  cluster <- findCluster(obj)
  length(unique(cluster)) * obj$vb  
}

estfun.rma.mv <- function(obj) {
  cluster <- droplevels(as.factor(findCluster(obj)))
  res <- residuals(obj)
  WX <- chol2inv(chol(obj$M)) %*% obj$X
  rval <- by(cbind(res, WX), cluster, 
             function(x) colSums(x[,1] * x[,-1, drop = FALSE]))
  rval <- matrix(unlist(rval), length(unique(cluster)), obj$p, byrow=TRUE)
  colnames(rval) <- colnames(obj$X)
  rval
}

RobustResults <- function(obj, adjust = TRUE) {
  cluster <- findCluster(obj)  
  vcov. <- sandwich(obj, adjust = adjust)
  df. <- length(unique(cluster)) - obj$p
  coeftest(obj, vcov. = vcov., df = df.)
}