## 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:

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 sigma

_{1}^{2}/within-study variability on level 2high sigma

_{2}^{2}/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/sigma

_{1}^{2}on level 3: low between-study variability/sigma

_{2}^{2}

But this time, the variability estimate for sigma_{1}^{2} 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:

+ sigma_{1}^{2} = 0.0312623

+ sigma_{2}^{2} = 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.

sigma_{1}^{2} = 0.0161855

sigma_{2}^{2} = 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.

## 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: tau^{2} for the model calculated with `rma.uni()`

and sigma^{2} with `rma.mv()`

+ no Higgins I^{2} 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
```

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 sigma_{1}^{2} and sigma_{2}^{2} 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

- 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!

- You can specify level 3 in various ways, I would recommend using either DOIs or the unique study names.

- 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**.