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. 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, as the sampling error in level 1 is assumed to be known. 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:

full.model<-rma.mv(y,
v,
random = list(~ 1 | ID_2,
~ 1 | ID_1),
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.0561  0.2368     80     no    ID_2
## sigma^2.2  0.1575  0.3968     14     no    ID_1
##
## 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’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. For 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. 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.distribution.3lm, 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.distribution.3lm <- function(data, m){

# Calculate estimated sampling variance and proportions across levels
data <- data
m <- m
n <- length(data$v) vector.inv.var <- 1/(data$v)
sum.inv.var <- sum(vector.inv.var)
sum.sq.inv.var <- (sum.inv.var)^2
vector.inv.var.sq <- 1/(data$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 level1<-((est.samp.var)/(m$sigma2[1]+m$sigma2[2]+est.samp.var)*100) level2<-((m$sigma2[1])/(m$sigma2[1]+m$sigma2[2]+est.samp.var)*100)
level3<-((m$sigma2[2])/(m$sigma2[1]+m$sigma2[2]+est.samp.var)*100) Level<-c("level 1", "level 2", "level 3") Variance<-c(level1, level2, level3) df<-data.frame(Level, Variance) df1<-df colnames(df1) <- c("Level", "% of total variance") # Generate plot df$distribution<-"Distribution"
df$Level<-factor(df$Level, levels(df\$Level)[c(3,2,1)])
g <- ggplot(df, aes(fill=Level, y=Variance, x=distribution)) +
geom_bar(stat="identity", position="fill", width = 0.3) + coord_flip(ylim = c(1,0)) + scale_y_continuous(labels = scales::percent)+
theme(axis.title.x=element_blank(),
axis.title.y=element_blank(),
axis.ticks.y=element_blank(),
axis.text.y = element_blank(),
axis.line.x = element_line(colour = "black",
size = 0.5, linetype = "solid"),
legend.position = "bottom",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
legend.background = element_rect(linetype="solid",
colour ="black"),
legend.title = element_blank(),
legend.key.size = unit(0.75,"cm")) + guides(fill = guide_legend(reverse = TRUE))
return(list(g, df1))
}

Please be aware that the function needs ggplot2 installed and loaded to work. There are only two parameters: data, the data we used for our three-level model (in which v signifies the effect size), and m, the model we just fitted using rma.mv. Let’s have a look:

library(ggplot2)
variance.distribution.3lm(data = mlm.data,
m = full.model)
## [[1]]

##
## [[2]]
##     Level % of total variance
## 1 level 1            8.937849
## 2 level 2           23.916209
## 3 level 3           67.145942

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.

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 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 verse for the third level.

Model without level 2

model.l2.removed<-rma.mv(y,
v,
random = list(~ 1 | ID_2,
~ 1 | ID_1),
tdist = TRUE,
data = mlm.data,
method = "REML",
sigma2 = c(0, NA))
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.0000  0.0000     80    yes    ID_2
## sigma^2.2  0.1620  0.4025     14     no    ID_1
##
## 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’s do the same when holding level 3 constant.

Model without level 3

model.l3.removed<-rma.mv(y,
v,
random = list(~ 1 | ID_2,
~ 1 | ID_1),
tdist = TRUE,
data = mlm.data,
method = "REML",
sigma2 = c(NA, 0))
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.3119  0.5585     80     no    ID_2
## sigma^2.2  0.0000  0.0000     14    yes    ID_1
##
## 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. doi:10.18637/jss.v067.i01.

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.