11.2 Fitting a three-level model

The metafor package is particularly well suited for fitting various three-level models in meta-analyses.

11.2.1 Data preparation

For this chapter, we continue to use the curry dataset from metaforest, which we have assigned to the object df. It has already been prepared for multilevel model fitting: Let’s have a peek at the dataset.

# Load the package, if you haven't yet
library(metaforest)
# Assign the curry dataset to df, if you haven't yet
df <- curry
# Examine the first 6 rows of the data
head(df)

We see that the first two columns of the dataset are ID variables:

study_id and effect_id. In these two columns, we stored the identifiers corresponding to the individual effect sizes, and we will use the studies as clusters on different levels of our multilevel model.

Multilevel models can accommodate different sources of heterogeneity, too. Maybe you have effect sizes originating in different countries, or maybe you want to cluster by (team of) authors. Just replace the study_id with your cluster id, in those cases. It is even possible to have more levels added (e.g., effect sizes within studies within countries), although that is beyond the scope of this tutorial.

11.2.2 Model fitting

We are now prepared to fit our model. I will save the model to the object m_multi. For model fitting, we use the rma.mv function in metafor. This is a more technical interface than the standard rma function. 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
yi	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
data	The data.frame 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 the random effects in our model. As we actually have two formulas in three-level models, we have to bundle them together in a list(). The notation for rma.mv is very similar to other packages specialized for mixed-effects models such as lme4 (Bates et al. 2015).

Within the list(), formulae are separated with commas. Within the formula, we first define the random effects variance as ~ 1, followed by a vertical bar |. Behind the vertical bar, we assign this random effect to a grouping variable such as studies, measures, regions and so forth. We only have to define the formulae for level 2 and level 3. The sampling variance in level 1 is assumed to be known; this is a necessary assumption in three-level meta-analysis, because we don’t have the raw data. Without this assumption, the model would not be identified. We type in the formulae in the order of our multilevel model. As in our example, ID_2 is nested in ID_1, we first define level 2 as ~ 1 | ID_2 and then level 3 as ~ 1 | ID_1. Now, let’s put it all together:

m_multi <- rma.mv(d, vi, random = list(~ 1 | effect_id, ~ 1 | study_id), data = df) 
m_multi

## 
## Multivariate Meta-Analysis Model (k = 56; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed     factor 
## sigma^2.1  0.0000  0.0000     56     no  effect_id 
## sigma^2.2  0.0838  0.2895     27     no   study_id 
## 
## Test for Heterogeneity:
## Q(df = 55) = 156.9109, p-val < .0001
## 
## Model Results:
## 
## estimate      se    zval    pval   ci.lb   ci.ub 
##   0.2836  0.0650  4.3638  <.0001  0.1562  0.4109  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Let’s go through the output one by one. First, let’s have a look at the Variance Components. Here, we see the variance calculated for each level of our model, \(\sigma^2_1\) (level 2) and \(\sigma^2_2\) (level 3). Under nlvls, we see the number of levels (i.e., subgroups) on each level. The number 56 corresponds to the number of unique effect sizes, and 27 is the number of unique studies. The heading factor is the best clue as to what levels these are: effect_id is the within-studies variance componend, and study_id is the between-studies variance component.

Under Model Results, we see the estimate of our pooled effect, which is 0.2836. As our y column contained effect sizes expressed as Hedges’ g, the effect is therefore \(g = 0.2836.\), with a confidence interval of \(0.16-0.41\).

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

11.2.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. Cheung (2014) provides a formula to estimate the sampling variance, which is needed to calculate the variance proportions for each level. We have have programmed a function called variance_decomposition, which implements this formula. The code for the function can be seen below. Parts of this code were adapted from Assink and Wibbelink (Assink and Wibbelink 2016).

R doesn’t know this function yet, so we have to let R learn it copying and pasting the code underneath in its entirety into the console on the bottom left pane of RStudio, and then hit Enter ⏎.

variance_decomposition <- function(m){
    n <- m$k
    vector.inv.var <- 1/(diag(m$V))
    sum.inv.var <- sum(vector.inv.var)
    sum.sq.inv.var <- (sum.inv.var)^2
    vector.inv.var.sq <- 1/(diag(m$V)^2)
    sum.inv.var.sq <- sum(vector.inv.var.sq)
    num <- (n-1)*sum.inv.var
    den <- sum.sq.inv.var - sum.inv.var.sq
    est.samp.var <- num/den
    if(length(m$sigma2) > 2) stop("Cannot handle more than three levels.")
    total_var <- (sum(m$sigma2)+est.samp.var)/100
    Variance <- c(est.samp.var, m$sigma2)/total_var
    names(Variance) <- c("Level1", m$s.names)
    Variance
}

Let’s have a look:

variance_decomposition(m_multi)

##       Level1    effect_id     study_id 
## 2.445403e+01 1.443707e-08 7.554597e+01

From the outputs, we see that about 24.45% of the overall variance can be attributed to level 1, 0% to level 2 (within-study variance), and as much as 75.55% to level 3 (between-study variance).

11.2.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 at 0. 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 2, level3). So, if we want to hold level 2 constant while leaving the rest as is, we use sigma2 = c(0,NA), and vice versa for the third level. We can then use ANOVAs to compare the fit of these models.

m_within_null <- rma.mv(yi = d, V = vi, random = list(~ 1 | effect_id, ~ 1 | study_id), sigma2 = c(0, NA), data = df) 
m_between_null <- rma.mv(yi = d, V = vi, random = list(~ 1 | effect_id, ~ 1 | study_id), sigma2 = c(NA, 0), data = df)
m_both_null <- rma.mv(yi = d, V = vi, random = list(~ 1 | effect_id, ~ 1 | study_id), sigma2 = c(0, 0), data = df) 

aov_within <- anova(m_multi, m_within_null) 
aov_between <- anova(m_multi, m_between_null) 
aov_bothnull <- anova(m_multi, m_both_null) 

# Join these results in a table
aov_table <- rbind(
c(df=aov_between$p.f, aov_between$fit.stats.f[c(3:4, 1)], LRT = NA, p = NA),
c(df=aov_within$p.r, aov_within$fit.stats.r[c(3:4, 1)], LRT = aov_within$LRT, p = aov_within$pval),
c(df=aov_between$p.r, aov_between$fit.stats.r[c(3:4, 1)], LRT = aov_between$LRT, p = aov_between$pval),
c(df=aov_bothnull$p.r, aov_bothnull$fit.stats.r[c(3:4, 1)], LRT = aov_bothnull$LRT, p = aov_bothnull$pval)
)
rownames(aov_table) <- c("Three-level model", "Within-studies variance constrained", "Between-studies variance constrained", "Both variance components constrained")

	df	AIC	BIC	ll	LRT	p
Three-level model	3	22.14	28.16	-8.07
Within-studies variance constrained	2	20.14	24.15	-8.07	0.00	1.000
Between-studies variance constrained	2	32.74	36.76	-14.37	12.61	0.000
Both variance components constrained	1	78.24	80.25	-38.12	60.10	0.000

This table allows us to easily compare the models, using fit indices (AIC and BIC; lower is better), and likelihood ratio tests and p-values. The results unanimously indicate that the model with within-studies variance constrained to 0 fits the best, and that this constraint does not lead to significantly worse fit (p = 1.00). Still, the decision which model to use should be based at least as much on theory, as on statistical criteria.

References

Assink, Mark, and Carlijn JM Wibbelink. 2016. “Fitting Three-Level Meta-Analytic Models in R: A Step-by-Step Tutorial.” The Quantitative Methods for Psychology 12 (3): 154–74.

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.