In a previous post, I provided some code to do robust variance estimation with metafor and sandwich. Here’s another example, replicating some more of the calculations from Tanner-Smith & Tipton (2013). (See here for the complete code.)

As a starting point, here are the results produced by the robumeta package:

Code

library(grid)library(robumeta)data(corrdat)rho <-0.8HTJ <-robu(effectsize ~ males + college + binge,data = corrdat, modelweights ="CORR", rho = rho,studynum = studyid,var.eff.size = var, small =FALSE)HTJ

RVE: Correlated Effects Model
Model: effectsize ~ males + college + binge
Number of studies = 39
Number of outcomes = 172 (min = 1 , mean = 4.41 , median = 4 , max = 18 )
Rho = 0.8
I.sq = 75.08352
Tau.sq = 0.1557714
Estimate StdErr t-value dfs P(|t|>) 95% CI.L 95% CI.U Sig
1 X.Intercept. 0.31936 0.27784 1.149 35 0.258 -0.2447 0.88340
2 males -0.00331 0.00376 -0.882 35 0.384 -0.0109 0.00431
3 college 0.41226 0.18685 2.206 35 0.034 0.0329 0.79159 **
4 binge 0.13774 0.12586 1.094 35 0.281 -0.1178 0.39326
---
Signif. codes: < .01 *** < .05 ** < .10 *
---

To exactly re-produce the results with metafor, I’ll need to use the weights proposed by HTJ. In their approach to the correlated effects case, effect size \(i\) from study \(j\) receives weight equal to \(\left[\left(v_{\cdot j} + \hat\tau^2\right)(1 + (k_j - 1) \rho)\right]^{-1}\), where \(v_{\cdot j}\) is the average sampling variance of the effect sizes from study \(j\), \(\hat\tau^2\) is an estimate of the between-study variance, \(k_j\) is the number of correlated effects in study \(j\), and \(\rho\) is a user-specified value of the intra-study correlation. However, it appears that robumeta actually uses a slightly different set weights, which are equivalent to taking \(\rho = 1\). I calculate the latter weights, fit the model in metafor, and output the robust standard errors and \(t\)-tests:

Code

devtools::source_gist(id ="11144005", filename ="metafor-sandwich.R")corrdat <-within(corrdat, { var_mean <-tapply(var, studyid, mean)[studyid] k <-table(studyid)[studyid] var_HTJ <-as.numeric(k * (var_mean +as.numeric(HTJ$mod_info$tau.sq)))})meta1 <-rma.mv(effectsize ~ males + college + binge, V = var_HTJ, data = corrdat, method ="FE")meta1$cluster <- corrdat$studyidRobustResults(meta1)

t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
intrcpt 0.3193586 0.2778360 1.1494 0.25816
males -0.0033143 0.0037573 -0.8821 0.38374
college 0.4122631 0.1868489 2.2064 0.03401 *
binge 0.1377393 0.1258637 1.0944 0.28127
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

One could specify a similar (though not exactly identical model) in metafor as follows. In the HTJ approach, \(\rho\) represents the total correlation induced by both the within-study sampling error and intra-study correlation in true effects. In contrast, the metafor approach would take \(\rho\) to be correlation due to within-study sampling error alone. I’ll first need to create a block-diagonal covariance matrix given a user-specified value of \(\rho\).

where \(Var(\nu_i) = \sigma^2\), \(Var(e_{ij}) = v_{ij}\), and \(Cor(e_{ij}, e_{ik}) = \rho\). Note that \(\sigma^2\) is now estimated via REML.

Code

meta2 <-rma.mv(yi = effectsize ~ males + college + binge, V = covMat, random =~1| studyid, data = corrdat,method ="REML")c(sigma.sq = meta2$sigma2)

sigma.sq
0.2477825

The between-study heterogeneity estimate is considerably larger than the moment estimate from robumeta. Together with the difference in weighting, this leads to some changes in the coefficient estimates and their estimated precision:

Code

RobustResults(meta2)

t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
intrcpt -0.8907096 0.4148219 -2.1472 0.038783 *
males 0.0163074 0.0055805 2.9222 0.006052 **
college 0.3180139 0.2273396 1.3988 0.170658
binge -0.0984026 0.0897269 -1.0967 0.280265
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

It is important to keep in mind that the estimate of between-study heterogeneity depends on the posited model for the covariance structure, including the assumed value of \(\rho\). HTJ recommend conducting sensitivity analysis across a range of values for the within-study effect correlation. Re-calculating the value of \(\sigma^2\) for \(\rho\) between 0.0 and 0.9 yields the following:

Code

sigma2 <-function(rho) { covMat <-as.matrix(bdiag(with(corrdat, tapply(var_mean, studyid, equicorr, rho = rho, simplify =FALSE))))rma.mv(yi = effectsize ~ males + college + binge, V = covMat, random =~1| studyid, data = corrdat,method ="REML")$sigma2}rho_sens <-seq(0,0.9,0.1)sigma2_sens <-sapply(rho_sens, sigma2)cbind(rho = rho_sens, sigma2 =round(sigma2_sens, 4))

The between-study heterogeneity is quite insensitive to the assumed value of \(\rho\).

The difference between the results based on metafor versus on robumeta appears to be due to the subtle difference in the weighting approach: metafor uses block-diagonal weights that contain off-diagonal terms for effects drawn from a common study, whereas robumeta uses entirely diagonal weights.

---title: Another meta-sandwichdate: '2014-04-23'categories:- meta-analysis- sandwiches- Rstats- robust variance estimationcode-tools: true---In [a previous post](/posts/Robust-meta-analysis-1/), I provided some code to do robust variance estimation with `metafor` and `sandwich`. Here's another example, replicating some more of the calculations from [Tanner-Smith & Tipton (2013)](http://doi.org/10.1002/jrsm.1091). ([See here](https://gist.github.com/jepusto/11147304) for the complete code.)As a starting point, here are the results produced by the `robumeta` package:```{r}library(grid)library(robumeta)data(corrdat)rho <-0.8HTJ <-robu(effectsize ~ males + college + binge,data = corrdat, modelweights ="CORR", rho = rho,studynum = studyid,var.eff.size = var, small =FALSE)HTJ```To exactly re-produce the results with `metafor`, I'll need to use the weights proposed by HTJ. In their approach to the correlated effects case, effect size $i$ from study $j$ receives weight equal to $\left[\left(v_{\cdot j} + \hat\tau^2\right)(1 + (k_j - 1) \rho)\right]^{-1}$, where $v_{\cdot j}$ is the average sampling variance of the effect sizes from study $j$, $\hat\tau^2$ is an estimate of the between-study variance, $k_j$ is the number of correlated effects in study $j$, and $\rho$ is a user-specified value of the intra-study correlation. However, it appears that `robumeta` actually uses a slightly different set weights, which are equivalent to taking $\rho = 1$. I calculate the latter weights, fit the model in `metafor`, and output the robust standard errors and $t$-tests:```{r, message=FALSE}devtools::source_gist(id = "11144005", filename = "metafor-sandwich.R")corrdat <- within(corrdat, { var_mean <- tapply(var, studyid, mean)[studyid] k <- table(studyid)[studyid] var_HTJ <- as.numeric(k * (var_mean + as.numeric(HTJ$mod_info$tau.sq)))})meta1 <- rma.mv(effectsize ~ males + college + binge, V = var_HTJ, data = corrdat, method = "FE")meta1$cluster <- corrdat$studyidRobustResults(meta1)```One could specify a similar (though not exactly identical model) in `metafor` as follows. In the HTJ approach, $\rho$ represents the total correlation induced by both the within-study sampling error and intra-study correlation in true effects. In contrast, the `metafor` approach would take $\rho$ to be correlation due to within-study sampling error alone. I'll first need to create a block-diagonal covariance matrix given a user-specified value of $\rho$. ```{r}library(Matrix)equicorr <-function(x, rho) { corr <- rho + (1- rho) *diag(nrow =length(x))tcrossprod(x) * corr } covMat <-as.matrix(bdiag(with(corrdat, tapply(var_mean, studyid, equicorr, rho =0.8, simplify =FALSE))))```Passing this block-diagonal covariance matrix to `rma.mv`, I now estimate the model $$T_{ij} = \mathbf{X}_{ij} \beta + \nu_i + e_{ij},$$where $Var(\nu_i) = \sigma^2$, $Var(e_{ij}) = v_{ij}$, and $Cor(e_{ij}, e_{ik}) = \rho$. Note that $\sigma^2$ is now estimated via REML.```{r}meta2 <-rma.mv(yi = effectsize ~ males + college + binge, V = covMat, random =~1| studyid, data = corrdat,method ="REML")c(sigma.sq = meta2$sigma2)```The between-study heterogeneity estimate is considerably larger than the moment estimate from `robumeta`. Together with the difference in weighting, this leads to some changes in the coefficient estimates and their estimated precision:```{r}RobustResults(meta2)```It is important to keep in mind that the estimate of between-study heterogeneity depends on the posited model for the covariance structure, including the assumed value of $\rho$. HTJ recommend conducting sensitivity analysis across a range of values for the within-study effect correlation. Re-calculating the value of $\sigma^2$ for $\rho$ between 0.0 and 0.9 yields the following:```{r}sigma2 <-function(rho) { covMat <-as.matrix(bdiag(with(corrdat, tapply(var_mean, studyid, equicorr, rho = rho, simplify =FALSE))))rma.mv(yi = effectsize ~ males + college + binge, V = covMat, random =~1| studyid, data = corrdat,method ="REML")$sigma2}rho_sens <-seq(0,0.9,0.1)sigma2_sens <-sapply(rho_sens, sigma2)cbind(rho = rho_sens, sigma2 =round(sigma2_sens, 4))```The between-study heterogeneity is quite insensitive to the assumed value of $\rho$. The difference between the results based on `metafor` versus on `robumeta` appears to be due to the subtle difference in the weighting approach: `metafor` uses block-diagonal weights that contain off-diagonal terms for effects drawn from a common study, whereas `robumeta` uses entirely diagonal weights.