Skip to content

Instantly share code, notes, and snippets.

@jepusto
Created April 21, 2014 14:16
Show Gist options
  • Save jepusto/11144005 to your computer and use it in GitHub Desktop.
Save jepusto/11144005 to your computer and use it in GitHub Desktop.
Calculate sandwich covariance estimators for multi-variate meta-analysis models
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.)
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment