[R-meta] Questions about Robust Variance Models using Meta-for

[James Pustejovsky]

Emily,

There are several differences between the models and estimation techniques
used in robumeta versus in metafor:

1. With the default option of modelweights = "CORR", robumeta estimates a
between-cluster variance component, assuming that effect sizes are nested
within clusters. If Yij is the ith effect size from cluster j, the model is

Yij = mu + uj + eij,

where Var(eij) = Vij, the sampling variance of the ith effect size from
cluster j; ej is a between-cluster random effect with Var(ej) = tau^2; and
effect size estimates from the same cluster are assumed to be correlated,
so corr(ehj, eij) = rho. So this model is a bit different than the model
you fit with rma.uni, which (again with the default options) assumes that
there is a between-study variance component but that every effect size is
from a different study. So generally speaking, the tau^2 estimate from
robumeta will not be the same as the tau^2 estimate from rma.uni. It will
typically be closer to a model fit using rma.mv.

2. robumeta uses a special formula to determine the weights and also
averages the sampling variances across effect sizes from a given cluster.
Say that cluster j includes a total of kj effect sizes and let V*j be the
average of the variances from cluster j, so V*j = (1 / kj) sum_{i=1}^{kj}
Vijj. Then robumeta uses weights

Wij = 1 / [kj *(V*j + tau^2)].

In contrast, rma.uni uses weights given by Wij = 1 / (Vij + tau^2).

3. robumeta uses a generally accurate small-sample correction (called
"bias-reduced linearization") when calculating cluster-robust standard
errors, whereas the robust() function in metafor uses a more approximate
small-sample correction. The small-sample adjustment in metafor should be
in about the right ballpark for basic meta-analysis models without any
covariates, but degrades in meta-regression models that include covariates.
The clubSandwich package calculates robust standard errors with a variety
of different small-sample corrections, including the one used in robumeta,
and can be used in combination with metafor.

All that is to say that it's possible to make robumeta and rma.uni agree
exactly, but it takes a bit of work. I give an example below.

Cheers,
James


library(robumeta)
library(metafor)
library(clubSandwich)

# fit a simple meta-analysis model using robumeta

data(corrdat)
corrdat$studyid <- factor(corrdat$studyid)

robu_fit <- robu(effectsize ~ 1, data = corrdat,
                 modelweights = "CORR", studynum = studyid,
                 var.eff.size = var)
robu_fit

# note value of tau^2
tau_sq_robu <- as.numeric(robu_fit$mod_info$tau.sq)

# fit the same specification using rma.uni
rma_uni_fit <- rma.uni(yi = effectsize, vi = var, data = corrdat)

# we get a different estimate of tau^2
rma_uni_fit$tau2

# robust standard error from metafor
robust(rma_uni_fit, cluster = corrdat$studyid, adjust = TRUE)

# we can get the same robust SE from clubSandwich using vcov = "CR1"
coef_test(rma_uni_fit, cluster = corrdat$studyid, vcov = "CR1")

# Using the more accurate small-sample approximation vcov = "CR2" leads to
slightly larger SE
coef_test(rma_uni_fit, cluster = corrdat$studyid, vcov = "CR2")

# calculate the inverse of the weight used by robumeta
corrdat$k <- with(corrdat, table(studyid)[studyid])
corrdat$Vbar <- with(corrdat, tapply(var, studyid, mean)[studyid])
corrdat$Vij <- with(corrdat, as.numeric(k * (Vbar + tau_sq_robu)))

# use the robumeta weights in rma.uni
rma_uni_robu <- rma.uni(yi = effectsize, vi = Vij, data = corrdat, method =
"FE")

# average effect size estimates now agree
rma_uni_robu$beta
robu_fit$reg_table$b.r

# we have to use "CR2" to get the standard errors to agree
coef_test(rma_uni_robu, cluster = corrdat$studyid, vcov = "CR2")
sqrt(vcovCR(rma_uni_robu, cluster = corrdat$studyid, type = "CR2"))
robu_fit$reg_table$SE

# the robumeta model is closer to the following multivariate model
V_mat <- with(corrdat, impute_covariance_matrix(vi = var, cluster =
studyid, r = 0.7))
rma_mv_fit <- rma.mv(yi = effectsize, V = V_mat, data = corrdat,
                     random = ~ studyid)
coef_test(rma_mv_fit, cluster = corrdat$studyid, vcov = "CR2")


On Tue, Sep 26, 2017 at 4:36 PM, Emily Quinn <quinnem at ohsu.edu> wrote:

> Dear Dr. Wolfgang Viechtbauer and other meta-analysists,
>
> I hope that this email finds you well.
> I was writing for your advice re: an issue related to the metafor R
> package.
>
> Currently, I am estimating a robust variance model using metafor rma.uni
> and robust functions and I noticed that the estimations were different than
> the estimates when I ran the same model in STATA using the robumeta package
> and in R using Robumeta See:
>
> STATA Robumeta:  hedges g= 0.3694, SE=0.1326, UL CI= 0.0939 - 0.6450
> R Robumeta: hedges g=0.369, SE=0.1333, UL CI = 0.0939- 0.645
> R metafor: hedges g = 0.307, SE 0.1149, UL 0.0736- 0.5479.
>
> Specifically, I was wondering if there was any arguments that I can
> specify that would make the metafor estimates closer to/ more consistent
> with the STATA and R Robumeta estimates:
>
> This is the syntax that I'm using to estimate the mode in metafor:
> library(metafor)
> Emily.t2 = rma.uni(yi=E$g,vi=E$var_g)
> Emily.t3 = robust(Emily.t2,E$author,adjust=T)
> Emily.t3
>
> Thank you so much for your time and efforts in advance.
>
> Sincerely,
> Emily D. Quinn
>
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-meta-analysis mailing list
> R-sig-meta-analysis at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-meta-analysis
>

	[[alternative HTML version deleted]]