Skip to content

Instantly share code, notes, and snippets.

@Yefeng0920
Forked from jepusto/metafor-BRL.R
Created February 25, 2024 14:22
Show Gist options
  • Save Yefeng0920/4204f1ea0a2feaa1deeb709e91146ba4 to your computer and use it in GitHub Desktop.
Save Yefeng0920/4204f1ea0a2feaa1deeb709e91146ba4 to your computer and use it in GitHub Desktop.
Bias-reduced linearization covariance estimator and degrees of freedom for multi-variate meta-analysis models
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
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment