Picture by Daniela Holzer

How to set up the 3-level structure for meta-analyses

Understanding the rma.mv() function

Intro

This blog post is inspired by a recent mistake of mine while conducting a multilevel meta-analysis with the R package metafor.

Luckily, I was made aware of this during the preprint stage of a current paper and not post publication: Open Science ftw. In this post I will explain what went wrong and what I learnt about the proper way for conducting multilevel meta-analyses with metafor.

In particular, I explore how the way we define our level 2 and level 3 influences the outcome of our multilevel meta-analysis, especially on its heterogeneity estimates. If you are interested more generally in how to conduct multilevel meta-analyses in R, I recently wrote a very short primer on this here. You can find a case study of multi-leveling up two papers without such analyses here.


1. How to set up level 2 IDs

Nested structures arise whenever the true effects within the same level could be correlated, i.e. within studies, geographic regions like districts, labs… This was the case in my meta-analysis, as effect sizes from the same published study are probably more similar than effect sizes from different published studies. Or put differently: the effect sizes are nested within the studies.

You can see in the fake data frame generic_data below that, for instance, 3 effect sizes/rows come from Study A, 2 effect sizes from Study B, etc. In order for metafor to conduct a multilevel meta-analysis properly, we have to specify this nested structure, and we could do so in multiple ways.

generic_data <- tribble(
  ~study_name , ~effect_size,  ~variance,
  "Study A",    0.30,      .033,
  "Study A",    0.37,      .030,
  "Study A",    0.39,      .020,
  "Study B",    0.08,      .050,
  "Study B",   -0.03,      .020,
  "Study C",    0.20,      .027,
  "Study D",    0.85,      .119,
  "Study E",    0.50,      .033,
  "Study E",    0.56,      .032,
  "Study E",    0.53,      .034,
  "Study E",    0.50,      .050,
  "Study E",    0.51,      .068
) %>% 
  mutate(es_id_good = row_number()) %>% # each effect size gets unique ID
  group_by(study_name) %>% 
  mutate(es_id_bad = row_number(), # each effect size gets nested ID
         study_id = group_indices()) 

generic_data
# A tibble: 12 x 6
# Groups:   study_name [5]
   study_name effect_size variance es_id_good es_id_bad study_id
   <chr>            <dbl>    <dbl>      <int>     <int>    <int>
 1 Study A           0.3     0.033          1         1        1
 2 Study A           0.37    0.03           2         2        1
 3 Study A           0.39    0.02           3         3        1
 4 Study B           0.08    0.05           4         1        2
 5 Study B          -0.03    0.02           5         2        2
 6 Study C           0.2     0.027          6         1        3
 7 Study D           0.85    0.119          7         1        4
 8 Study E           0.5     0.033          8         1        5
 9 Study E           0.56    0.032          9         2        5
10 Study E           0.53    0.034         10         3        5
11 Study E           0.5     0.05          11         4        5
12 Study E           0.51    0.068         12         5        5

And here is where my mistake happened:

For the package metafor to work properly, each effect size needs to have a UNIQUE ID (es_id_good), and not a nested ID (es_id_bad).

To make things as clear as possible, I call the wrong way to code effect size IDs “bad” and the right way “good” throughout this article. You can see that es_id_good is unique for every row, but es_id_bad has repeated values.

This repetition of effect size IDs only makes sense for multivariate meta-analyses! Effect sizes with the same value of the ID receive the same random effect, while effects sizes with unique IDs are assumed to be independent. For example, we would give a random effect the same IDs, e.g. ranging from 1-5, if we wanted to specify IQ domains:

  • 1 = full scale
  • 2 = verbal
  • 3 = numeric
  • 4 = spatial
  • 5 = other

But if you want to include effect sizes that are nested within studies, you need to use unique effect size IDs!

If we simply calculate the 3-level model with the bad effect size IDs, we would get the following output:

bad <- rma.mv(yi = effect_size, 
              V = variance,                      
              random = list(~ 1 | es_id_bad,     # level 2
                            ~ 1 | study_id),     # level 3
              data = generic_data,
              method = "REML")
bad

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

Variance Components:

            estim    sqrt  nlvls  fixed     factor 
sigma^2.1  0.0000  0.0000      5     no  es_id_bad 
sigma^2.2  0.0446  0.2113      5     no   study_id 

Test for Heterogeneity:
Q(df = 11) = 15.5978, p-val = 0.1567

Model Results:

estimate      se    zval    pval   ci.lb   ci.ub 
  0.3303  0.1157  2.8559  0.0043  0.1036  0.5570  ** 

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

You can see under nlvls that we have only 5 cluster for level 2 and 5 cluster for level 3. But we included k = 12 effect sizes! We didn´t nest the effect sizes in the cluster properly!

When we use the correct es_id the output looks different:

good <- rma.mv(yi = effect_size, 
              V = variance,                      # sampling variances
              random = list(~ 1 | es_id_good, # indicate level 2
                            ~ 1 | study_id),     # indicate level 3
              data = generic_data,
              method = "REML")
good

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

Variance Components:

            estim    sqrt  nlvls  fixed      factor 
sigma^2.1  0.0000  0.0000     12     no  es_id_good 
sigma^2.2  0.0446  0.2113      5     no    study_id 

Test for Heterogeneity:
Q(df = 11) = 15.5978, p-val = 0.1567

Model Results:

estimate      se    zval    pval   ci.lb   ci.ub 
  0.3303  0.1157  2.8559  0.0043  0.1036  0.5570  ** 

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

Now you can see that under nlvls we have 12 effect sizes included for level 2 and 5 cluster for level 3. This is what we want to see for our multilevel meta-analyses.

But for some reason, the variability estimates stayed the same. This is simply because of the way I set up the fake data set. Normally this would not happen. Here, the effect sizes that are nested within these studies are very similar, but the studies are dissimilar!

In such a scenario we would expect that we find:

  • low sigma12/within-study variability on level 2

  • high sigma22/between-study variability on level 3

In our example we get the same within-study variability of zero no matter how we define the effect size ID! This is simply the case because there is no within cluster variability to start with!

Example when it makes it difference

Lets have a look what happens if the effect sizes nested within studies are different, but studies are similar.

Expectation:

  • on level 2: high within-study variability/sigma12

  • on level 3: low between-study variability/sigma22

But this time, the variability estimate for sigma12 should be different and the way we code level 2 should make a difference.

Here I create a new data frame with high within-study variability and low between-study variability:

data_b <- tribble(
  ~study_name , ~effect_size,  ~variance,
  "Study A",          0.30,      .033,
  "Study A",          0.60,      .030,
  "Study A",          0.00,      .020,
  "Study B",          0.62,      .050,
  "Study B",         -0.03,      .020,
  "Study C",          0.30,      .027,
  "Study D",          0.35,      .119,
  "Study E",          0.30,      .033,
  "Study E",          0.56,      .032,
  "Study E",          0.13,      .034,
  "Study E",          0.60,      .050,
  "Study E",          0.01,      .068
) %>% 
  mutate(es_id_good = row_number()) %>% 
  group_by(study_name) %>% 
  mutate(es_id_bad = row_number(),
         study_id = group_indices())

If we use the correct definition for effect size IDs we get the following output:

b_good <- rma.mv(yi = effect_size, 
                 V = variance,                  
                 random = list(~ 1 | es_id_good, #  level 2
                               ~ 1 | study_id),  #  level 3
                 data = data_b,
                 method = "REML")
b_good

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

Variance Components:

            estim    sqrt  nlvls  fixed      factor 
sigma^2.1  0.0313  0.1768     12     no  es_id_good 
sigma^2.2  0.0000  0.0000      5     no    study_id 

Test for Heterogeneity:
Q(df = 11) = 20.7415, p-val = 0.0362

Model Results:

estimate      se    zval    pval   ci.lb   ci.ub 
  0.2952  0.0754  3.9172  <.0001  0.1475  0.4429  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
b_good$sigma2[1]
[1] 0.0312623

As expected, our level 2 variability is high, while on level 3 it is basically non existent:
+ sigma12 = 0.0312623
+ sigma22 = 1.052799110^{-12}

b_bad <- rma.mv(yi = effect_size, 
            V = variance,                  # sampling variances
            random = list(~ 1 | es_id_bad,     # indicate level 2
                          ~ 1 | study_id), # indicate level 3
            tdist = TRUE,                  # Knapp-Hartung adjustment
            data = data_b,
            method = "REML",
            slab = data_b$study_name)

b_bad

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

Variance Components:

            estim    sqrt  nlvls  fixed     factor 
sigma^2.1  0.0162  0.1272      5     no  es_id_bad 
sigma^2.2  0.0000  0.0000      5     no   study_id 

Test for Heterogeneity:
Q(df = 11) = 20.7415, p-val = 0.0362

Model Results:

estimate      se    tval    pval   ci.lb   ci.ub 
  0.2662  0.0835  3.1865  0.0087  0.0823  0.4500  ** 

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

The same pattern emerges if we specify level 2 incorrectly, our level 2 variability is high, while on level 3 it is basically non existent.

sigma12 = 0.0161855
sigma22 = 8.398503910^{-12}

BUT, the sigmas are different depending whether you specified the es_id correctly or incorrectly!


2. How to set up level 3 IDs

Now that we correctly specified level 2, here are some insights on specifying level 3.

It is enough to use a unique study_name—or anything unique for that matter—to define level-3. It could just be a unique character string or a unique number.

study_name on level 3

lvl3_study_name <- rma.mv(yi = effect_size, 
                          V = variance,                     
                          random = list(~ 1 | es_id_good, #  level 2
                                        ~ 1 | study_name),   #  level 3
                          data = generic_data,
                          method = "REML")
lvl3_study_name

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

Variance Components:

            estim    sqrt  nlvls  fixed      factor 
sigma^2.1  0.0000  0.0000     12     no  es_id_good 
sigma^2.2  0.0446  0.2113      5     no  study_name 

Test for Heterogeneity:
Q(df = 11) = 15.5978, p-val = 0.1567

Model Results:

estimate      se    zval    pval   ci.lb   ci.ub 
  0.3303  0.1157  2.8559  0.0043  0.1036  0.5570  ** 

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

study_id on level 3

lvl3_study_id <- rma.mv(yi = effect_size, 
                        V = variance,                      
                        random = list(~ 1 | es_id_good, # level 2
                                      ~ 1 | study_id),  # level 3
                        data = generic_data,
                        method = "REML")
lvl3_study_id

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

Variance Components:

            estim    sqrt  nlvls  fixed      factor 
sigma^2.1  0.0000  0.0000     12     no  es_id_good 
sigma^2.2  0.0446  0.2113      5     no    study_id 

Test for Heterogeneity:
Q(df = 11) = 15.5978, p-val = 0.1567

Model Results:

estimate      se    zval    pval   ci.lb   ci.ub 
  0.3303  0.1157  2.8559  0.0043  0.1036  0.5570  ** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
identical(lvl3_study_id$sigma2, lvl3_study_name$sigma2)
[1] TRUE

As you can see, the heterogeneity estimates are equivalent, so it is not necessary to use the study id here. BUT, the study name needs to be unique, so if you include papers from the same author who has multiple papers from the same year, you need to either call it something like: Study A 2017a, Study A 2017b… or use the study ID. OR what I will do from now on out is to just use the DOI.

Using DOIs as the identifier for level 3 makes it as clear as possible, which study is included in each analysis.


3. How to define the random-effects structure of the model

Additionally, I dove into different ways to define the random argument. The random argument is used to define a formula that specifies the random-effects structure of the model. This can either be a single one-sided formula or list of one-sided formulas.


Regular univariate meta-analysis

Lets start with an easy example. Normally, if we want to calculate a “standard/2-level” random effects meta-analysis, we would use rma.uni():

rma.uni(yi = effect_size, 
        vi = variance, 
        data = generic_data)

Random-Effects Model (k = 12; tau^2 estimator: REML)

tau^2 (estimated amount of total heterogeneity): 0.0155 (SE = 0.0211)
tau (square root of estimated tau^2 value):      0.1245
I^2 (total heterogeneity / total variability):   31.09%
H^2 (total variability / sampling variability):  1.45

Test for Heterogeneity:
Q(df = 11) = 15.5978, p-val = 0.1567

Model Results:

estimate      se    zval    pval   ci.lb   ci.ub 
  0.3553  0.0653  5.4393  <.0001  0.2273  0.4834  *** 

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

But we can do the same with rma.mv(). The random = ~ 1 | es_id_good argument adds random effects corresponding to the effect size ID to the model.

uni <- rma.mv(yi = effect_size, 
        V = variance, 
        random = ~ 1 | es_id_good, # each es_id gets random effect 
        data=generic_data)

uni 

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

Variance Components:

            estim    sqrt  nlvls  fixed      factor 
sigma^2    0.0155  0.1245     12     no  es_id_good 

Test for Heterogeneity:
Q(df = 11) = 15.5978, p-val = 0.1567

Model Results:

estimate      se    zval    pval   ci.lb   ci.ub 
  0.3553  0.0653  5.4393  <.0001  0.2273  0.4834  *** 

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

The only differences:
+ how the heterogeneity is called: tau2 for the model calculated with rma.uni() and sigma2 with rma.mv()
+ no Higgins I2 for the rma.mv() is provided, this needs to be calculated separately.


Single one-sided formula/nesting notation

Now we use the syntax random = ~ 1 | study_id/es_id_good to add random effects not only at the effect size level, but also at the study level.

single <- rma.mv(yi = effect_size, 
                 V = variance,                       
                 random = ~ 1 | study_id/es_id_good, # level3/level 2
                 data = generic_data,
                 method = "REML")

single 

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

Variance Components:

            estim    sqrt  nlvls  fixed               factor 
sigma^2.1  0.0446  0.2113      5     no             study_id 
sigma^2.2  0.0000  0.0000     12     no  study_id/es_id_good 

Test for Heterogeneity:
Q(df = 11) = 15.5978, p-val = 0.1567

Model Results:

estimate      se    zval    pval   ci.lb   ci.ub 
  0.3303  0.1157  2.8559  0.0043  0.1036  0.5570  ** 

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


Update based on a tweet by Wolfgang Viechtbauer: You can also use the nested ID (es_id_bad) in this type of nested notation!

https://twitter.com/wviechtb/status/1251175269156741121

nested_notation <- rma.mv(yi = effect_size, 
                 V = variance,                       
                 random = ~ 1 | study_id/es_id_bad, # level3/level 2
                 data = generic_data,
                 method = "REML")

nested_notation 

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

Variance Components:

            estim    sqrt  nlvls  fixed              factor 
sigma^2.1  0.0446  0.2113      5     no            study_id 
sigma^2.2  0.0000  0.0000     12     no  study_id/es_id_bad 

Test for Heterogeneity:
Q(df = 11) = 15.5978, p-val = 0.1567

Model Results:

estimate      se    zval    pval   ci.lb   ci.ub 
  0.3303  0.1157  2.8559  0.0043  0.1036  0.5570  ** 

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

List of one-sided formulas

You could also do it like this, which is equivalent to the previous method, but makes it IMO clearer, which level is level 2 and which is level 3.

list <- rma.mv(yi = effect_size, 
              V = variance,                      
              random = list(~ 1 | es_id_good,    # level 2
                            ~ 1 | study_id),     # level 3
              data = generic_data,
              method = "REML")
list

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

Variance Components:

            estim    sqrt  nlvls  fixed      factor 
sigma^2.1  0.0000  0.0000     12     no  es_id_good 
sigma^2.2  0.0446  0.2113      5     no    study_id 

Test for Heterogeneity:
Q(df = 11) = 15.5978, p-val = 0.1567

Model Results:

estimate      se    zval    pval   ci.lb   ci.ub 
  0.3303  0.1157  2.8559  0.0043  0.1036  0.5570  ** 

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

You can see that this yields the exact same results, but be aware that the order of sigma12 and sigma22 is now switched!


inner | outer notation

You could also use inner | outer notation, which is described in more detail here:

“In addition or alternatively to specifying one or multiple ~ 1 | id terms, the random argument can also contain a formula of the form ~ inner | outer. Effects or outcomes with different values/levels of the outer grouping variable/factor are assumed to be independent, while effects or outcomes with the same value/level of the outer grouping variable/factor share correlated random effects corresponding to the levels of the inner grouping variable/factor (note that the inner grouping variable must be either a factor or a character variable).”

generic_data <- generic_data %>% 
  mutate(es_id_good = as.factor(es_id_good)) # needs to be a factor
inner_outer <- rma.mv(yi = effect_size, 
       V = variance,                          # sampling variances
       random = ~ es_id_good | study_id, # indicate level 3 |ßlevel 2
       data = generic_data,
       method = "REML")
inner_outer

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

Variance Components:

outer factor: study_id   (nlvls = 5)
inner factor: es_id_good (nlvls = 12)

            estim    sqrt  fixed 
tau^2      0.0446  0.2113     no 
rho        1.0000             no 

Test for Heterogeneity:
Q(df = 11) = 15.5978, p-val = 0.1567

Model Results:

estimate      se    zval    pval   ci.lb   ci.ub 
  0.3303  0.1157  2.8559  0.0043  0.1036  0.5570  ** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
inner_outer$sigma2s
[1] 1

As you can see, we get the same results, yet the variability is not separated for level-2 and level-3, which would be nice to have, so I wouldn´t recommend using this notation!


Conclusion

  1. The way you specify nested effect sizes at level 2 makes a big difference. Due to this subtle mistake, I had to check ~15.000 lines of code to make sure level 2 was specified correctly for each analysis. All analyses, tables, and figures had to be changed, so do not mess this up. Trust me!
  2. You can specify level 3 in various ways, I would recommend using either DOIs or the unique study names.
  3. There are multiple notations for multilevel models, I prefer lists of one-sided formulas, as they make things as clear as possible. You could also use the nested notation with single one-sided formulas.