#' @title Robust random-effect meta-regression with dependent effect sizes
#'
#' @description
#' \code{rrma} returns robust variance estimates from random effect
#' meta-regression with dependent effect-size estimates.
#'
#' @details
#' This is an implementation based on the \code{robust.se} function from
#' Hedges et al. (2010).
#'
#' @references
#' Hedges, L.V., Tipton, E. & Johnson, M.C. (2010). Robust variance
#' estimation in meta-regression with dependent effect_size estimates.
#' Res. Synth. Method., 1, 39-65.
#'
#' @author Main body of code is from the appendix of Hedges et al.
#' (2010) and written by the paper authors.
#' Code adapted to an R package with a formula interface and various
#' convenience functions by Sean Anderson and Jarrett Byrnes.
#'
#' @param formula The meta-regression formula of the form \code{y ~ x1 +
#' x2...}, where \code{y} is the effect size. An intercept only meta-analysis
#' can be performed with \code{y ~ 1}. The intercept can be excluded with
#' \code{y ~ -1 + x1 + x2 ...}
#' @param data The input data frame.
#' @param study_id The study IDs. Can be in any form (character,
#' numeric, factor). Will be converted to factor and then numeric form
#' internally.
#' @param var_eff_size The variance on each effect size.
#' @param rho The assumed correlation between effect sizes of the same
#' study including correlation induced by the random effects. Hedges
#' et al. (2010) suggest conducting a sensitivity analysis on
#' \code{rho}.
#' @param level Tolerance/confidence interval
#'
#' @export
#' @return A list with the robust covariance matrix, the Qe statistic,
#' the tau_sq estimate (variance between studies), a data frame of
#' the robust coefficient estimates and standard errors, and a number
#' of diagnostics and input values.
#'
#' @examples
#' data(broad)
#' (m <- rrma(formula = lnorReg ~ d18OresidualMean.cent, data =
#' broad, study_id = study.ID, var_eff = vlnorReg, rho = 0.5))
#'
rrma <- function(formula, data, study_id, var_eff_size, rho, level = 0.05) {
# [r]obust [r]andom effect [m]eta-[a]nalysis regression Hedges et al. 2010
# code to parse the formula and create the required data structure:
require(plyr) # using plyr to save time because I already coded it this way
cl <- match.call()
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "data", "study_id", "var_eff_size"),
names(mf), 0L)
mf <- mf[c(1L, m)]
mf$drop.unused.levels <- TRUE
mf[[1L]] <- as.name("model.frame")
mf <- eval(mf, parent.frame())
mf <- data.frame(effect_size = mf[,1], model.matrix(formula, mf),
study_id=mf[["(study_id)"]],
var_eff_size=mf[["(var_eff_size)"]])
names(mf) <- gsub("\\(", "", names(mf))
names(mf) <- gsub("\\)", "", names(mf))
names(mf) <- gsub("X.Intercept.", "intercept", names(mf)) #for later use
ind_var_names <- names(mf)[-c(1, match(c("study_id", "var_eff_size"), names(mf)))]
mf <- ddply(mf, "study_id", transform, k = length(study_id),
mean_v = mean(var_eff_size), s = sqrt(var_eff_size))
mf <- ddply(mf, "study_id", transform, weights = 1/(k * mean_v))
#mf$intercept <- 1
names(mf)[1] <- "effect_size"
mf$study <- as.numeric(factor(mf$study_id))
names_design_matrix <- c("study", ind_var_names)
names_data_matrix <- c("study", "k", "mean_v", "weights",
"effect_size", "var_eff_size", "s")
# now create 2 data frames: input_data and X.full
input_data <- mf[,names_data_matrix]
X_full <- mf[,names_design_matrix]
# and proceed with the function from Hedges et al. 2010:
p <- ncol(X_full) - 2 # number of non-intercept betas
N <- max(input_data$study) # number of studies
sumXWX <- 0
sumXWy <- 0
sumXWJWX <- 0
sumXWVWX <- 0
sumXW.sig.m.v.WX <- 0
for (i in (1:N)) { # loop through each study
tab <- input_data[input_data$study == i, ] # get the meta-analysis data for this study
W <- diag(tab$weights, tab$k[1]) # weights are 1 / (k * mean_v)
tab2 <- X_full[X_full$study == i, ] # design matrix for this study
tab3 <- cbind(tab2[-c(1)]) # remove study column
X <- data.matrix(tab3) # convert to matrix
dimnames(X) <- NULL # clean column and row names
y <- cbind(tab$effect_size) # create effect_size vector
one <- cbind(rep(1, tab$k[1])) # a matrix of ones k times
J <- one %*% t(one) # a k x k matrix of ones
sigma <- (tab$s %*% t(tab$s)) # variance-covariance
vee <- diag(tab$s^2, tab$k[1]) # variance diagonals
SigmV <- sigma - vee # the covariance elements
sumXWX <- sumXWX + t(X) %*% W %*% X # first half of equation 3
sumXWy <- sumXWy + t(X) %*% W %*% y # second half of equation 3
sumXWJWX <- sumXWJWX + t(X) %*% W %*% J %*% W %*% X # equation 4; needed for denominator of equation 15
sumXWVWX <- sumXWVWX + t(X) %*% W %*% vee %*% W %*% X # needed for equation 15
sumXW.sig.m.v.WX <- sumXW.sig.m.v.WX + t(X) %*% W %*% SigmV %*% W %*% X # needed for equation 15
}
b <- solve(sumXWX) %*% sumXWy # solve for the (not-robust) beta means
X <- data.matrix(X_full[-c(1)])
dimnames(X) <- NULL
input_data$pred <- X %*% b # predictions
input_data$e <- input_data$effect_size - input_data$pred # e = residuals
W <- diag(input_data$weights) # matrix with weights on diagonal
sumW <- sum(input_data$weights) # sum of all weights across all effects
Qe <- t(input_data$e) %*% W %*% input_data$e # equation 14, the Qe statistic
denom <- sumW - sum(diag(solve(sumXWX) %*% sumXWJWX)) # denominator of equation 15
term1 <- (Qe - N + sum(diag(solve(sumXWX) %*% sumXWVWX)))/denom # part of equation 15
term2 <- (sum(diag(solve(sumXWX) %*% sumXW.sig.m.v.WX)))/denom # part of equation 15
tau.sq1 <- term1 + rho * term2 # the full equation 15
tau.sq <- ifelse(tau.sq1 < 0, 0, tau.sq1) # don't allow negative tau
input_data$r.weights <- 1/(input_data$k * (input_data$mean_v + tau.sq)) # weights adjusted for tau.sq; [r]obust weights; equation 17
sumXWX.r <- 0
sumXWy.r <- 0
for (i in (1:N)) { # now repeat beta estimation with robust weights
tab <- input_data[input_data$study == i, ]
W <- diag(tab$r.weights, tab$k[1])
tab2 <- X_full[X_full$study == i, ]
tab3 <- cbind(tab2[-c(1)])
X <- data.matrix(tab3)
dimnames(X) <- NULL
y <- cbind(tab$effect_size)
sumXWX.r <- sumXWX.r + t(X) %*% W %*% X
sumXWy.r <- sumXWy.r + t(X) %*% W %*% y
}
b.r <- solve(sumXWX.r) %*% sumXWy.r # solve for robust beta means
X <- data.matrix(X_full[-c(1)])
dimnames(X) <- NULL
input_data$pred.r <- X %*% b.r # predictions from robust betas
input_data$e.r <- cbind(input_data$effect_size) - input_data$pred.r # residuals from robust betas
sumXWeeWX.r <- 0
for (i in (1:N)) { # now get estimated variance-covariance matrix using the robust residuals
tab <- input_data[input_data$study == i, ]
sigma.hat.r <- tab$e.r %*% t(tab$e.r)
W <- diag(tab$r.weights, tab$k[1])
tab2 <- X_full[X_full$study == i, ]
tab3 <- cbind(tab2[-c(1)])
X <- data.matrix(tab3)
dimnames(X) <- NULL
sumXWeeWX.r <- sumXWeeWX.r + t(X) %*% W %*% sigma.hat.r %*%
W %*% X
}
VR.r <- solve(sumXWX.r) %*% sumXWeeWX.r %*% solve(sumXWX.r) # the robust (approximated) variance-covariance matrix
SE <- c(rep(0, p + 1))
for (i in (1:(p + 1))) { # now convert the variances (the diagonal) to standard errors
SE[i] <- sqrt(VR.r[i, i]) * sqrt(N/(N - (p + 1)))
}
labels <- names(X_full)[-1]
output <- data.frame(beta = labels, estimate = b.r, SE = SE)
#fill out the table with z-scores and CI information to match rma
crit <- qt(level/2, df=(N-(p+1)), lower.tail=FALSE )
output <- cbind(output, with(output, data.frame(
tval = estimate/SE,
pval = pt(abs(estimate/SE), df = (N - (p + 1)), lower.tail=F)*2,
ci.lb = estimate - crit*SE,
ci.ub = estimate + crit*SE
)))
rownames(VR.r) <- colnames(VR.r) <- output[,2]
input_data$study_id <- mf$study_id
ret <- list(VR_r = VR.r, Qe = Qe, `tau_sq_est` = tau.sq, `est` =
output, W=W, X_full = X_full, input_data=input_data,
rho=rho, mf=mf, call=cl, n=nrow(input_data), pred =
input_data$pred.r, res = input_data$e.r, k = N, yi =
input_data$effect_size, vi = input_data$var_eff_size,
df = N-(p+1), crit = crit, level = level, individual_weights = 1/input_data$var_eff_size, weights = input_data$weights, robust_weights = input_data$r.weights)
class(ret) <- "rrma"
ret
}
data(broad)
m <- rrma(formula = lnorReg ~ d18OresidualMean.cent, data =
broad, study_id = study.ID, var_eff = vlnorReg, rho = 0.5)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.