Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Correlated effects example from Tanner-Smith & Tipton (2013)
# 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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment