Aggregating or not?

With complex data structures in meta-analysis sometimes we want to aggregate effect size to simplify the analysis. As reported by James Pustejovsky, “sometimes, aggregating effect sizes is fine”. The blog post show formally why aggregating or not, in specific condition, brings the exactly same results with some additional pros in terms of computational speed. Crucially, this is true only in specific conditions.

In this post I’ll highlight one such circumstance, where aggregating effect size estimates is not only reasonable but leads to exactly the same results as a multivariate model. This occurs when two conditions are met: - We are not interested in within-study heterogeneity of effects - Any predictors included in the model vary between studies but not within a given study (i.e., effect sizes from the same study all have the same values of the predictors). In short, if all we care about is understanding between-study variation in effect sizes, then it is fine to aggregate them up to the study level.

The example refers to a situation where I have multiple effects collected on the same pool of subjects nested within papers. Thus if I have predictors at the level of the effect or I’m interested in modelling variability within papers, aggregating is not good.

In case where I do not have moderators or I’m interested only in between-papers variability, aggregating is totally fine.

We can expand this even to situations with more nesting level, outcomes nested within experiments nested within papers. If we are not interested in the effect sizes level (i.e., outcome level) aggregating is fine.

Example

Here I simulate a dataset with \(K\) papers. Each papers can have a maximum of \(J\) experiments and each experiments can have 1 (e.g., ACC or RT) or 2 (e.g., ACC and RT) outcomes colected on the same pool of subjects.

library(metafor)

Loading required package: Matrix

Loading required package: metadat

Loading required package: numDeriv


Loading the 'metafor' package (version 4.0-0). For an
introduction to the package please type: help(metafor)

seqw <- function(x){
  unlist(sapply(x, function(i) 1:i))
}
set.seed(2023)

We start by simulating the data structure. The number of experiments and outcomes is random:

K <- 100 # number of papers
J <- sample(1:3, K, replace = TRUE) # number of studies, within each paper
Z <- sample(1:2, sum(J), replace = TRUE) # number of outcomes per study/paper

dat <- data.frame(paper = rep(rep(1:K, J), Z), 
                  exp = rep(seqw(J), Z), 
                  effect = seqw(Z))

head(dat)

  paper exp effect
1     1   1      1
2     1   1      2
3     2   1      1
4     2   1      2
5     2   2      1
6     2   3      1

Then we set our simulation parameters. We have 3 heterogeneity components: \(\tau^2\), \(\omega^2\) and \(\zeta^2\). I simulate a meta-regression model because I simulate a difference between the effect on the two outcomes but I’m not interested in fitting a meta-regression (just a shortcut to simulate the difference).

# residual variance components
tau2 <- 0.3
omega2 <- 0.1
zeta2 <- 0.1

b0 <- 0.1
b1 <- 0.1

# random effects
b0_i <- rnorm(K, 0, sqrt(tau2))
b0_ij <- rnorm(sum(J), 0, sqrt(omega2))
b0_ijz <- rnorm(nrow(dat), 0, sqrt(zeta2))

# add to dataframe
dat$b0_i <- b0_i[dat$paper]
dat$b0_ij <- rep(b0_ij, Z)
dat$b0_ijz <- b0_ijz
dat$vi <- runif(nrow(dat), 0.05, 0.1) # random sampling variances

head(dat)

  paper exp effect        b0_i       b0_ij      b0_ijz         vi
1     1   1      1 -0.04382947 -0.23146442 -0.31300488 0.08194773
2     1   1      2 -0.04382947 -0.23146442  0.18361480 0.05893909
3     2   1      1  0.15294661  0.09919803  0.45521629 0.06651081
4     2   1      2  0.15294661  0.09919803  0.01169994 0.06953170
5     2   2      1  0.15294661  0.32233407  0.20481287 0.09750496
6     2   3      1  0.15294661  0.15360518  0.04071347 0.06072859

Now the crucial part, we need to create the block-variance-covariance matrix of the sampling errors. Crucially, sampling errors are correlated ONLY when we have multiple outcomes because experiments within the same paper have different subjects. The matrix should reflect this feature.

# create block-matrix
V <- vcalc(vi, 
           cluster = paper, # 1st level cluster
           subgroup = exp,  # independent experiments
           obs = effect, # correlated effects
           rho = 0.7, # correlation
           data = dat)

# splitting the matrix just for reference
Vb <- blsplit(V, dat$paper, fun = round, 3)

For example, the second paper has 3 experiments. The first experiment has two outcomes while the second and the third have only one outcome.

dat[dat$paper == 2, ]

  paper exp effect      b0_i      b0_ij     b0_ijz         vi
3     2   1      1 0.1529466 0.09919803 0.45521629 0.06651081
4     2   1      2 0.1529466 0.09919803 0.01169994 0.06953170
5     2   2      1 0.1529466 0.32233407 0.20481287 0.09750496
6     2   3      1 0.1529466 0.15360518 0.04071347 0.06072859

Thus sampling errors are correlated only for experiment 1, let’s see the matrix:

Vb[[2]]

      [,1]  [,2]  [,3]  [,4]
[1,] 0.067 0.048 0.000 0.000
[2,] 0.048 0.070 0.000 0.000
[3,] 0.000 0.000 0.098 0.000
[4,] 0.000 0.000 0.000 0.061

Now we can use V to simulate the sampling errors:

# sampling errors
e_ij <- MASS::mvrnorm(1, mu = rep(0, nrow(V)), Sigma = V)

Finally we can simulate the different effect size for each outcome:

# moderator (ACC, RT)
dat$x <- ifelse(dat$effect == 1, 1, 0)

# observed effect size
dat$yi <- with(dat, (b0 + b0_i + b0_ij + b0_ijz) + b1*x + e_ij)
dat$x <- factor(dat$x)
dat$exp <- factor(dat$exp)

head(dat)

  paper exp effect        b0_i       b0_ij      b0_ijz         vi x         yi
1     1   1      1 -0.04382947 -0.23146442 -0.31300488 0.08194773 1 -0.5126605
2     1   1      2 -0.04382947 -0.23146442  0.18361480 0.05893909 0 -0.1950280
3     2   1      1  0.15294661  0.09919803  0.45521629 0.06651081 1  0.9493747
4     2   1      2  0.15294661  0.09919803  0.01169994 0.06953170 0  0.1391810
5     2   2      1  0.15294661  0.32233407  0.20481287 0.09750496 1  1.2265437
6     2   3      1  0.15294661  0.15360518  0.04071347 0.06072859 1  0.6342383

Model without aggregation

To compare the model with and without aggregation, the trick is fitting a model without the lowest nesting level in the random effect structure. Normally, the simulated dataset above should be modeled (for the random part) as 1|paper/exp/effect. Here we drop the effect term:

fit0 <- rma.mv(yi, V, random = ~1|paper/exp, data = dat, sparse = TRUE)

summary(fit0)


Multivariate Meta-Analysis Model (k = 293; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-463.2780   926.5559   932.5559   943.5862   932.6393   

Variance Components:

            estim    sqrt  nlvls  fixed     factor 
sigma^2.1  0.3506  0.5921    100     no      paper 
sigma^2.2  0.1651  0.4063    193     no  paper/exp 

Test for Heterogeneity:
Q(df = 292) = 2486.0170, p-val < .0001

Model Results:

estimate      se    zval    pval    ci.lb   ci.ub    
  0.0010  0.0699  0.0147  0.9883  -0.1360  0.1381    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Using the V matrix we are fixing the correlation between multiple outcomes (i.e., the lowest nesting level) to 0.7.

Model with aggregation

To aggregate we use the Borenstein method implemented in the aggregate function:

# with aggregation
dat <- escalc(yi = yi, vi = vi, data = dat)
datl <- split(dat, dat$paper)
datl <- lapply(datl, function(x) aggregate(x, cluster = exp, rho = 0.7))
datagg <- do.call(rbind, datl)

head(datagg)


    paper exp effect        b0_i       b0_ij      b0_ijz     vi x      yi 
1       1   1    1.5 -0.04382947 -0.23146442 -0.06469504 0.0565 1 -0.2700 
2.1     2   1    1.5  0.15294661  0.09919803  0.23345811 0.0578 1  0.5742 
2.2     2   2    1.0  0.15294661  0.32233407  0.20481287 0.0975 1  1.2265 
2.3     2   3    1.0  0.15294661  0.15360518  0.04071347 0.0607 1  0.6342 
3       3   1    1.0 -0.00408931  0.23980867 -0.20986505 0.0874 1  0.0607 
4.1     4   1    1.5  0.25987891 -0.17661214  0.19485017 0.0701 1  0.4567

Then we fit the same model but without using V (sampling errors are never correlated now) instead we use v and the same random structure:

fit1 <- rma.mv(yi, vi, random = ~1|paper/exp, data = datagg, sparse = TRUE)

summary(fit1)


Multivariate Meta-Analysis Model (k = 193; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-196.5763   393.1527   399.1527   408.9251   399.2803   

Variance Components:

            estim    sqrt  nlvls  fixed     factor 
sigma^2.1  0.3506  0.5921    100     no      paper 
sigma^2.2  0.1651  0.4063    193     no  paper/exp 

Test for Heterogeneity:
Q(df = 192) = 1825.2162, p-val < .0001

Model Results:

estimate      se    zval    pval    ci.lb   ci.ub    
  0.0010  0.0699  0.0147  0.9883  -0.1360  0.1381    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The two models are exactly the same.