Correlated effects example from Tanner-Smith & Tipton (2013)

RVE-correlated.R

# robumeta calculations
library(grid)
library(robumeta)
data(corrdat)
rho <- 0.8

HTJ <- robu(effectsize ~ males + college + binge,
            data = corrdat, 
            modelweights = "CORR", rho = rho,
            studynum = studyid,
            var.eff.size = var, small = FALSE)
HTJ

# reproduce using metafor
source("R/metafor-sandwich.R") # See https://gist.github.com/jepusto/11144005

corrdat <- within(corrdat, {
  var_mean <- tapply(var, studyid, mean)[studyid]
  k <- table(studyid)[studyid]
  var_HTJ <- as.numeric(k * (var_mean + as.numeric(HTJ$tau.sq)))
})

meta1 <- rma.mv(effectsize ~ males + college + binge, 
                V = var_HTJ, 
                data = corrdat, method = "FE")
meta1$cluster <- corrdat$studyid
RobustResults(meta1)


# fit model using metafor
rho = 0.0
library(Matrix)
equicorr <- function(x, rho) {
  corr <- rho + (1 - rho) * diag(nrow = length(x))
  tcrossprod(x) * corr 
} 
covMat <- as.matrix(bdiag(with(corrdat, tapply(var_mean, studyid, equicorr, rho = rho, simplify = FALSE))))
meta2 <- rma.mv(yi = effectsize ~ males + college + binge, 
                V = covMat, random = ~ 1 | studyid, 
                data = corrdat,
                method = "REML")
meta2$sigma2
RobustResults(meta2)

# sensitivity analysis for between-study heterogeneity
sigma2 <- function(rho) {
  covMat <- as.matrix(bdiag(with(corrdat, tapply(var_mean, studyid, equicorr, rho = rho, simplify = FALSE))))
  rma.mv(yi = effectsize ~ males + college + binge, 
                  V = covMat, random = ~ 1 | studyid, 
                  data = corrdat,
                  method = "REML")$sigma2
}
rho_sens <- seq(0,0.9,0.1)
sigma2_sens <- sapply(rho_sens, sigma2)
round(sigma2_sens, 4)