## 12.1 Fitting a three-level model

### 12.1.1 Data preparation

Before we start with fitting the model, we need the `metafor`

package to be loaded from our library.

`library(metafor)`

For this chapter, i will use the `mlm.data`

dataset. It has already been prepared for model fitting using `metafor`

. Let’s have a peek at the dataset:

`mlm.data`

We see that the dataset has 80 rows and 5 columns (i.e., 5 variables). The first four are the most important ones:

**ID_1**and**ID_2**. In these two columns, we stored the identifiers corresponding to the individual effect sizes and clusters on different levels of our multilevel model. We have left these columns quite generic to make clear that such multilevel data can accommodate various types of data structures. For example, one can think of a scenario where ID_2 corresponds with the different effect sizes we find within a study (e.g., different effect sizes for different measures of depression). In this example, the effect sizes in ID_2 would then be nested within the larger clusters defined in ID_1, the studies. Another possible example would be that ID_2 in fact symbolizes the unique identifier for single studies, which are then nested in larger clusters (e.g., cultural regions) shown in ID_1.**y**and**v**. In these columns, the effect size data is stored. The format is similar to the one we know from the`metagen`

function in`meta`

. The column`y`

contains the effect size (\(d\), \(g\), \(r\), \(z\), or other continuous effect size metrics). The`v`

column contains the sampling variance which can be obtained by**squaring the standard error**.

There are two additional columns which contain data needed for **subgroup analyses**. To conduct subgroup analyses with `metafor`

, however, we need to store the data a little differently than in Chapter 7. This time, we do not use a factor vector containing the subgroup levels as strings, but **recode all subgroup levels** as individuals columns in the data. In our case, we are interested if the publication status has a moderating effect on the pooled results. We want to compare studies published in peer-review journals to dissertations. To do this, we need to create **two variables/columns for each subgroup level** (peer-review, dissertation). In each column we provide a code if the effect size belongs to this subgroup, defined by `1 = yes`

and `0 = no`

. Of course, the subgroup codings now have to be **mutually exclusive**, as we cannot calculate subgroup effects if a study is both part of the “peer-review” group and the “dissertation” group.

### 12.1.2 Model fitting

We are now prepared to fit our model. I will save the model to the object `full.model`

. For model fitting, we use the `rma.mv`

function in `metafor`

. The function has plenty of parameters one can specify (type `?metafor::rma.mv`

in the console to see them all). For now, we will only focus on the really essential ones.

Code | Description |
---|---|

y | The variable in our dataset containing the effect size data |

v | The variable in our dataset containing the sampling variance data corresponding to each y |

random | The model specifications for our multilevel model. We describe this in more detail below |

tdist | This is an old friend in disguise, the Knapp-Hartung adjustment for our confidence intervals (if set to TRUE) |

data | The dataframe containing our data |

method | The tau-squared estimator. Our options here are limited, and it is advised to use the restricted maximum likelihood (‘REML’) method, which is the default |

**Setting up the formula**

Under the `random`

parameter, we feed `rma.mv`

with the formula for our model. The notation for `rma.mv`

is quite similar to other packages specialized for mixed-effects models such as `lme4`

(Bates et al. 2015). We start defining the random effects by using `~ 1`

, followed by a **vertical bar** `|`

. Behind the vertical bar, we assign a random effect to a **grouping variable** such as studies, measures, regions or the like. This grouping variable often called a **random intercept**, because it tells our model to assume different values (intercepts) for each grouping. We have described before that in the multilevel model, we have two grouping variables: the one on **level 2**, and the one on **level 3**. We assume that these grouping variables are **nested**, in the sense that several effect sizes on level 2 (represented by `ID_2`

) together make up a larger cluster on level 3 (represented by `ID_1`

). There is a special way through which we can tell `rma.mv`

that we assume that our random effects have a nested structure. This is achieved by using a slash (`/`

) to seperate the higher- and lower-order grouping variables. Left to `/`

, we define the higher-order variable or cluster (here, `ID_1`

). To the right, we define the lower-order variable which is nested in the larger cluster (in our example, `ID_2`

). All in all, our formula therefore looks like this: `~ 1 | ID_1/ID_2`

. The full code for our call to `rma.mv`

looks like this:

```
full.model <- rma.mv(y,
v,
random = ~ 1 | ID_1/ID_2,
tdist = TRUE,
data = mlm.data,
method = "REML")
```

To see the results of our model, we can use the `summary()`

function.

`summary(full.model)`

```
##
## Multivariate Meta-Analysis Model (k = 80; method: REML)
##
## logLik Deviance AIC BIC AICc
## -44.6962 89.3924 95.3924 102.5007 95.7124
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.1575 0.3968 14 no ID_1
## sigma^2.2 0.0561 0.2368 80 no ID_1/ID_2
##
## Test for Heterogeneity:
## Q(df = 79) = 624.8063, p-val < .0001
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## 0.4345 0.1173 3.7058 0.0004 0.2011 0.6679 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

Let’s go through the output one by one. First, let us have a look at the `Variance Components`

. Here, we see the random-effects variance calculated for each level of our model, \(\sigma^2_1\) (level 3, `ID_1`

) and \(\sigma^2_2\) (level 2, `ID_2`

nested in `ID_1`

). Under `nlvls`

, we see the number of levels (i.e., subgroups) on each level. For `ID_1/ID_2`

, this is 80 because each effect size has its own ID.

Under `Model Results`

, we see the `estimate`

of our pooled effect, which is 0.4345. As our `y`

column contained effect sizes expressed as Hedges’ g, the effect is therefore \(g = 0.4545\), with a confidence interval of \(0.20-0.67\).

Under `Test for Heterogeneity`

, we see that the variance of effect sizes within our data overall was highly significant (\(p < 0.0001\)). This however, is not very informative as we are interested how variance can be attributed to the different levels in our model.

### 12.1.3 Distribution of total variance

As mentioned before, what we’re actually interested in is the **distribution of variance over the three levels of our model**. Viechtbauer (2019) provides a formula to estimate the sampling variance, which is needed to calculate the variance proportions for each level. We have have prepared a function called `mlm.variance.distribution`

, which implements this formula. The function is part of the `dmetar`

package. If you have the package installed already, you have to load it into your library first.

`library(dmetar)`

If you don’t want to use the `dmetar`

package, you can find the source code for this function here. In this case, *R* doesn’t know this function yet, so we have to let *R* learn it by **copying and pasting** the code **in its entirety** into the **console** in the bottom left pane of RStudio, and then hit **Enter ⏎**. The function then requires the `ggplot2`

package to work.

We only have to specify one parameter: `x`

, the model we just fitted using `rma.mv`

. Let’s have a look:

`mlm.variance.distribution(x = full.model)`

```
## % of total variance I2
## Level 1 8.937849 ---
## Level 2 23.916209 23.92
## Level 3 67.145942 67.15
## Total I2: 91.06%
```

From the outputs, we see that about 8.94% of the overall variance can be attributed to level 1, 23.92% to level 2, and as much as 67.15% to level 3. The variance proportions at level 2 and 3 are equivalent to the between-study heterogeneity \(I^2\) we covered in Chapter 6. However, in the context of three-level meta-analysis models, there are three types of \(I^2\) values: \(I^2_{Level2}\), which describes the amount of within-cluster heterogeneity, \(I^2_{Level3}\), representing the amount of between-cluster heterogeneity, and \(I^2_{total}\), which is the total amount of heterogeneity not attributable to sampling error. \(I^2_{total}\) is simply the sum of the \(I^2\) values at level 2 and level 3.

### 12.1.4 Comparing the fit

Of course, when fitting a three-level model, we are interested if this model actually **captures the variability in our data better than a two-level model**. `metafor`

allows us to easily do this by comparing models **in which one of the levels is removed**.

To do this, we can use the `rma.mv`

function again, but this time, we hold the variance component \(\sigma^2\) of one level constant. This can be done by specifying the `sigma2`

parameter. The parameter can be specified by providing a vector telling the function which level to hold constant, with the generic version being `c(level 3, level 2)`

. So, if we want to hold level 2 constant while leaving the rest as is, we use `sigma2 = c(NA, 0)`

, and vice verse for the third level.

**Model without level 2**

```
model.l2.removed<-rma.mv(y,
v,
random = ~ 1 | ID_1/ID_2,
tdist = TRUE,
data = mlm.data,
method = "REML",
sigma2 = c(NA, 0))
summary(model.l2.removed)
```

```
##
## Multivariate Meta-Analysis Model (k = 80; method: REML)
##
## logLik Deviance AIC BIC AICc
## -59.8748 119.7495 123.7495 128.4884 123.9074
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.1620 0.4025 14 no ID_1
## sigma^2.2 0.0000 0.0000 80 yes ID_1/ID_2
##
## Test for Heterogeneity:
## Q(df = 79) = 624.8063, p-val < .0001
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## 0.4217 0.1130 3.7337 0.0004 0.1969 0.6465 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

We see that the `Model Results`

have changed, but does the second level in fact capture a significant amount of variability in the data? To check, we can use the `anova`

function to compare the `full.model`

to the `model.l2.removed`

.

`anova(full.model, model.l2.removed)`

```
##
## df AIC BIC AICc logLik LRT pval QE
## Full 3 95.3924 102.5007 95.7124 -44.6962 624.8063
## Reduced 2 123.7495 128.4884 123.9074 -59.8748 30.3571 <.0001 624.8063
```

We see that the `Full`

model compared to the `Reduced`

model does indeed have a better fit, with the **Akaike Information Criterion** (`AIC`

) and the **Bayesian Information Criterion** (`BIC`

) being lower for this model. The difference is significant (\(p<0.0001\)), suggesting that it was a good idea to include this level into our analysis. Now, let us do the same when holding level 3 constant.

**Model without level 3**

```
model.l3.removed<-rma.mv(y,
v,
random = ~ 1 | ID_1/ID_2,
tdist = TRUE,
data = mlm.data,
method = "REML",
sigma2 = c(0, NA))
summary(model.l3.removed)
```

```
##
## Multivariate Meta-Analysis Model (k = 80; method: REML)
##
## logLik Deviance AIC BIC AICc
## -75.4151 150.8302 154.8302 159.5691 154.9881
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0000 0.0000 14 yes ID_1
## sigma^2.2 0.3119 0.5585 80 no ID_1/ID_2
##
## Test for Heterogeneity:
## Q(df = 79) = 624.8063, p-val < .0001
##
## Model Results:
##
## estimate se tval pval ci.lb ci.ub
## 0.6049 0.0691 8.7539 <.0001 0.4674 0.7424 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

`anova(full.model, model.l2.removed)`

```
##
## df AIC BIC AICc logLik LRT pval QE
## Full 3 95.3924 102.5007 95.7124 -44.6962 624.8063
## Reduced 2 123.7495 128.4884 123.9074 -59.8748 30.3571 <.0001 624.8063
```

We see the same pattern here, suggesting that overall, our three-level model is very useful.

### References

Bates, Douglas, Martin Mächler, Ben Bolker, and Steve Walker. 2015. “Fitting Linear Mixed-Effects Models Using lme4.” *Journal of Statistical Software* 67 (1): 1–48. https://doi.org/10.18637/jss.v067.i01.