12.1 Robust variance estimation

The package clubSandwich can perform robust variance estimation for a variety of models - including the rma models created by metafor. This is convenient, because it means that we can simply specify our model in the same way as before. Then, we take the output of our rma analysis, and use it as input for the robust variance estimation function clubSandwich::coef_test.

Robust variance estimation is another way to account for dependent data (multiple outcomes from one study/sample). Thus, we need to indicate a ‘clustering’ variable again. When we pass our rma output to coef_test, the function automatically uses the outermost level of our multilevel model.

We first re-fit the three-level multilevel meta-analysis from before, which uses REML to estimate the variance components:

m_multi <- rma.mv(yi = d, V = vi, random = list(~1 | effect_id, ~1 | study_id), data = df)
m_multi

## 
## Multivariate Meta-Analysis Model (k = 56; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed     factor 
## sigma^2.1  0.0000  0.0000     56     no  effect_id 
## sigma^2.2  0.0838  0.2895     27     no   study_id 
## 
## Test for Heterogeneity:
## Q(df = 55) = 156.9109, p-val < .0001
## 
## Model Results:
## 
## estimate      se    zval    pval   ci.lb   ci.ub 
##   0.2836  0.0650  4.3638  <.0001  0.1562  0.4109  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The standard errors, t-tests, and confidence intervals from rma are all parameteric and model-based. The clubSandwich::coef_test function will calculate robust standard errors and robust t-tests for each of the coefficients:

# Install if necessary
install.packages("clubSandwich")
library(clubSandwich)

clubSandwich::coef_test(m_multi, vcov = "CR2")

##     Coef. Estimate    SE t-stat d.f. p-val (Satt) Sig.
## 1 intrcpt    0.284 0.065   4.37 24.6       <0.001  ***

Note that coef_test assumed that it should cluster based on studyID, which is the outer-most random effect in the metafor model. This can be specified explicitly by including the option cluster = df$effect_id in the call.

The vcov = "CR2" option means that the standard errors will be corrected using the bias-reduced linearization estimator described in Tipton and Pustejovsky (2015).

---