Bias-reduced linearization covariance estimator and degrees of freedom for multi-variate meta-analysis models

metafor-BRL.R

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

#-----------------------------------------------
# Identify outer-most clustering variable
#-----------------------------------------------

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 method for rma.mv
#-----------------------------------------------

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

#-----------------------------------------------
# estfun method for rma.mv
#-----------------------------------------------

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
}

#--------------------------------------------------
# Bias-reduced linearization covariance estimator
# for rma.mv objects
#--------------------------------------------------

meatBRL <- function(obj, correction = "BRL-S") {
  cluster <- droplevels(as.factor(findCluster(obj)))
  reslist <- tapply(residuals(obj), cluster, function(x) x, simplify = FALSE)
  Vlist <- lapply(levels(cluster), function(c) obj$M[cluster==c, cluster==c, drop=FALSE])
  Wlist <- lapply(Vlist, function(v) chol2inv(chol(v)))
  Vchol <- lapply(Vlist, chol)
  Xlist <- lapply(levels(cluster), function(c) obj$X[cluster==c,,drop=FALSE])
  XQX <- lapply(Xlist, function(x) x %*% obj$vb %*% t(x))
  
  if (correction == "BRL-A") { 
    U_XQX_cholinv <- mapply(function(v,x) solve(chol(v - x)), v = Vlist, x = XQX)
    A_list <- mapply(function(u,m) t(u) %*% t(m), u = Vchol, m = U_XQX_cholinv)  
  } else {
    G_eig <- mapply(function(a,b,c) eigen(a %*% (b - c) %*% t(a)), 
                    a = Vchol, b = Vlist, c = XQX, SIMPLIFY = FALSE)
    G_half <- lapply(G_eig, function(g) g$vectors %*% diag(1 / sqrt(g$values), nrow=length(g$values)) %*% t(g$vectors))
    A_list <- mapply(function(g,v) t(v) %*% g %*% v, g = G_half, v = Vchol)
  }
  
  m <- nlevels(cluster)
  psi <- mapply(function(x,w,a,r) t(x) %*% w %*% a %*% r, 
                   x = Xlist, w = Wlist, a = A_list, r = reslist, SIMPLIFY = FALSE)
  psi <- matrix(unlist(psi), nrow = obj$p, ncol = m, dimnames = list(colnames(obj$X)))
  rownames(psi)
  rval <- tcrossprod(psi) / m
  rval
}

#--------------------------------------------------
# Satterthwaite degrees of freedom
# for BRL covariance with rma.mv objects
#--------------------------------------------------

dfBRL <- function(obj, correction = "BRL-S") {
  cluster <- droplevels(as.factor(findCluster(obj)))
  reslist <- tapply(residuals(obj), cluster, function(x) x, simplify = FALSE)
  Vlist <- lapply(levels(cluster), function(c) obj$M[cluster==c, cluster==c, drop=FALSE])
  Wlist <- lapply(Vlist, function(v) chol2inv(chol(v)))
  Vchol <- lapply(Vlist, chol)
  Xlist <- lapply(levels(cluster), function(c) obj$X[cluster==c,,drop=FALSE])
  XQX <- lapply(Xlist, function(x) x %*% obj$vb %*% t(x))
  
  if (correction == "BRL-A") { 
    U_XQX_cholinv <- mapply(function(v,x) solve(chol(v - x)), v = Vlist, x = XQX)
    A_list <- mapply(function(u,m) t(u) %*% t(m), u = Vchol, m = U_XQX_cholinv)  
  } else {
    G_eig <- mapply(function(a,b,c) eigen(a %*% (b - c) %*% t(a)), 
                    a = Vchol, b = Vlist, c = XQX, SIMPLIFY = FALSE)
    G_half <- lapply(G_eig, function(g) g$vectors %*% diag(1 / sqrt(g$values), nrow=length(g$values)) %*% t(g$vectors))
    A_list <- mapply(function(g,v) t(v) %*% g %*% v, g = G_half, v = Vchol)
  }
  
  IH <- diag(1, nrow=obj$k) - obj$X %*% obj$vb %*% t(obj$X) %*% chol2inv(chol(obj$M))
  IH_list <- lapply(levels(cluster), function(c) IH[cluster==c,,drop=FALSE])

  g_list <- mapply(function(ih, a, w, x) t(ih) %*% t(a) %*% w %*% x %*% obj$vb,
                   ih = IH_list, a = A_list, w = Wlist, x = Xlist, SIMPLIFY = FALSE)
  g_array <- array(unlist(g_list), dim = c(obj$k, obj$p, length(g_list)))
  V_eig <- eigen(obj$M)
  V_half <- V_eig$vectors %*% diag(sqrt(V_eig$values), nrow=length(V_eig$values)) %*% t(V_eig$vectors)
  df_eigen <- apply(g_array, 2, function(g) eigen(V_half %*% tcrossprod(g) %*% V_half, symmetric = TRUE, only.values=TRUE)$values)
  df <- apply(df_eigen, 2, function(x) sum(x)^2 / (sum(x^2)))
  names(df) <- rownames(obj$b)
  df
}

#-----------------------------------------------------
# robust t-tests for fixed effects in rma.mv objects
#-----------------------------------------------------
  
RobustResults <- function(obj, correction = "BRL-S") {
  if (correction == "HTJ") {
    vcov. <- sandwich(obj, adjust = TRUE)
    cluster <- findCluster(obj)  
    df. <- length(unique(cluster)) - obj$p
  } else {
    vcov. <- sandwich(obj, meat. = meatBRL, correction=correction)
    df. <- dfBRL(obj, correction = correction)
  }
  results <- coeftest(obj, vcov. = vcov., df = df.)
  results <- cbind(results[,1:3, drop=FALSE], df = df., results[,4,drop=FALSE])
  
  results
}