# The metafor Package

A Meta-Analysis Package for R

### Sidebar

analyses:konstantopoulos2011

## Konstantopoulos (2011)

### The Methods and Data

Konstantopoulos (2011) describes a three-level meta-analytic model that is applicable when multiple effect size estimates are nested within a higher level grouping variable (e.g., paper, lab, research group, species). One would then expect the underlying true effects to be more similar for the same level of the grouping variable than true effects arising from different levels (i.e., the true effects within the same level may be correlated). We can account for the correlation in the true effects by adding a random effect to the model at the level corresponding to the grouping variable.

In the paper, the use of such a multilevel model is illustrated based on a meta-analysis of studies examining the effects of modified school calendars on student achievement (Cooper et al., 2003). In particular, instead of following the more traditional school calendar with a long summer break (in addition to a short winter and spring break), some schools have switched to a modified school calendar comprising more frequent but shorter intermittent breaks (e.g., 9 weeks of school followed by 3 weeks off), while keeping the total number of days at school approximately the same.

The dataset (taken from Konstantopoulos, 2011) contains the results from 56 studies, each comparing the level of academic achievement in a group of students following a modified school calendar with that of a group of students following a more traditional school calendar. The difference between the two groups was quantified in terms of a standardized mean difference (with positive values indicating a higher mean level of achievement in the group following the modified school calendar).

The studies were conducted at various schools that were clustered within districts. The data therefore have a multilevel structure, with studies nested within districts. A multilevel meta-analysis of these data can be used to estimate and account for the amount of heterogeneity between districts and between schools within districts.

The data (provided in Table 1 of the paper) can be loaded with:

library(metafor)
dat <- dat.konstantopoulos2011

(I copy the dataset into dat, which is a bit shorter and therefore easier to type further below). Let's examine the first 15 rows of the dataset:

head(dat, 15)
   district school study year     yi    vi
1        11      1     1 1976 -0.180 0.118
2        11      2     2 1976 -0.220 0.118
3        11      3     3 1976  0.230 0.144
4        11      4     4 1976 -0.300 0.144
5        12      1     5 1989  0.130 0.014
6        12      2     6 1989 -0.260 0.014
7        12      3     7 1989  0.190 0.015
8        12      4     8 1989  0.320 0.024
9        18      1     9 1994  0.450 0.023
10       18      2    10 1994  0.380 0.043
11       18      3    11 1994  0.290 0.012
12       27      1    12 1976  0.160 0.020
13       27      2    13 1976  0.650 0.004
14       27      3    14 1976  0.360 0.004
15       27      4    15 1976  0.600 0.007 

So, 4 studies were conducted in district 11, 4 studies in district 12, 3 studies in district 18, and so on. Variables yi and vi are the standardized mean differences and corresponding sampling variances.

### Descriptives

On page 69, some descriptive information about the dataset is provided. Let's first reproduce this information. The studies were conducted in various years, ranging from 1976 to 2000:

table(dat$year) 1976 1989 1994 1995 1997 2000 19 4 5 2 15 11  The percentage of studies conducted in the various years can be obtained with: 100*round(prop.table(table(dat$year)),2)
1976 1989 1994 1995 1997 2000
34    7    9    4   27   20

We can obtain some descriptives about the effect size estimates with:

round(c(summary(dat$yi), SD=sd(dat$yi)), 2)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.      SD
-0.52   -0.05    0.04    0.12    0.28    1.19    0.33

To get a better sense of the data structure, we can examine how many studies were carried out in each of the 11 districts with:

table(dat$district)  11 12 18 27 56 58 71 86 91 108 644 4 4 3 4 4 11 3 8 6 5 4 So, as noted in the article, the data have an unbalanced structure, with the number of studies/schools per district ranging from 3 to 11 (range(table(dat$district))) with an average of 5.1 (round(mean(table(dat$district)), 1)). To obtain the descriptives about the effect size estimates per district (Table 3 in the paper), we can use: round(aggregate(yi ~ district, data=dat, FUN=function(x) c(mean=mean(x), sd=sd(x), min=min(x), max=max(x))), 3)  district yi.mean yi.sd yi.min yi.max 1 11 -0.118 0.237 -0.300 0.230 2 12 0.095 0.250 -0.260 0.320 3 18 0.373 0.080 0.290 0.450 4 27 0.442 0.227 0.160 0.650 5 56 0.062 0.103 -0.060 0.190 6 58 -0.077 0.126 -0.300 0.070 7 71 0.823 0.465 0.300 1.190 8 86 -0.029 0.042 -0.100 0.020 9 91 0.247 0.271 0.000 0.660 10 108 0.016 0.487 -0.520 0.700 11 644 0.180 0.295 -0.050 0.610 ### Two-Level Model First, a standard (two-level) random-effects model is fitted to the data. Here, we treat the 56 studies as independent (which we later will see is not justified). We can fit such a model with: res <- rma(yi, vi, data=dat) print(res, digits=3) Random-Effects Model (k = 56; tau^2 estimator: REML) tau^2 (estimated amount of total heterogeneity): 0.088 (SE = 0.020) tau (square root of estimated tau^2 value): 0.297 I^2 (total heterogeneity / total variability): 94.70% H^2 (total variability / sampling variability): 18.89 Test for Heterogeneity: Q(df = 55) = 578.864, p-val < .001 Model Results: estimate se zval pval ci.lb ci.ub 0.128 0.044 2.916 0.004 0.042 0.214 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 These results correspond to those given on the left-hand side of Table 4. The CI for$\tau^2$given in the paper (i.e.,$0.059–0.146$) is based on a Satterthwaite approximation. A CI for$\tau^2$based on the Q-profile method can be obtained with: confint(res, digits=3)  estimate ci.lb ci.ub tau^2 0.088 0.056 0.139 tau 0.297 0.238 0.373 I^2(%) 94.705 91.946 96.560 H^2 18.885 12.416 29.072 The results from the mixed-effects meta-regression model including the study year as a potential moderator/covariate can be reproduced with: res <- rma(yi, vi, mods = ~ I(year-mean(year)), data=dat) Mixed-Effects Model (k = 56; tau^2 estimator: REML) tau^2 (estimated amount of residual heterogeneity): 0.089 (SE = 0.020) tau (square root of estimated tau^2 value): 0.298 I^2 (residual heterogeneity / unaccounted variability): 94.71% H^2 (unaccounted variability / sampling variability): 18.89 R^2 (amount of heterogeneity accounted for): 0.00% Test for Residual Heterogeneity: QE(df = 54) = 550.260, p-val < .001 Test of Moderators (coefficient 2): QM(df = 1) = 1.383, p-val = 0.240 Model Results: estimate se zval pval ci.lb ci.ub intrcpt 0.126 0.044 2.859 0.004 0.040 0.212 ** I(year - mean(year)) 0.005 0.004 1.176 0.240 -0.003 0.014 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 The corresponding results can be found on the right-hand side of Table 4. Note that the year variable was entered in mean centered form in the model. Again, a CI for$\tau^2$can be obtained with confint(res, digits=3) (output not shown). We could also fit the same two-level random/mixed-effects models with the rma.mv() function. The syntax for the random-effects model would be: res <- rma.mv(yi, vi, random = ~ 1 | study, data=dat) print(res, digits=3) Multivariate Meta-Analysis Model (k = 56; method: REML) Variance Components: estim sqrt nlvls fixed factor sigma^2 0.088 0.297 56 no study Test for Heterogeneity: Q(df = 55) = 578.864, p-val < .001 Model Results: estimate se zval pval ci.lb ci.ub 0.128 0.044 2.916 0.004 0.042 0.214 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 The random = ~ 1 | study argument adds random effects corresponding to the study level to the model (which is a unique value for every row in the dataset, so we are essentially adding random effects for every estimate to the model). Note that these are the same results as we obtained earlier. As before, moderators/covariates can be added to the model via the mods argument. ### Three-Level Model Let's move on now to the three-level model. Here, we want to add random effects not only at the study level, but also at the district level. The appropriate syntax for this is: res.ml <- rma.mv(yi, vi, random = ~ 1 | district/study, data=dat) print(res.ml, digits=3) Multivariate Meta-Analysis Model (k = 56; method: REML) Variance Components: estim sqrt nlvls fixed factor sigma^2.1 0.065 0.255 11 no district sigma^2.2 0.033 0.181 56 no district/study Test for Heterogeneity: Q(df = 55) = 578.864, p-val < .001 Model Results: estimate se zval pval ci.lb ci.ub 0.185 0.085 2.185 0.029 0.019 0.350 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 These results correspond to those given on the left-hand side of Table 5 in the paper. Somewhat confusingly, the results given in the table appear to be based on ML (instead of REML) estimation. With the argument method="ML", we could reproduce the results given in the paper more closely. Note: We would obtain the same results when using random = ~ 1 | district/school. While the school variable always starts at 1 within each district, using this notation means that random effects should be added for each level of district and for each level of school within each level of district. The latter is the same as adding random effects for every level of study, so the results will be identical. ### Profile Likelihood Plots Whenever we start fitting more complicated models with the rma.mv() function, it is a good idea to check the profile likelihood plots of the variance components of the model. The profile() function can be used to obtain such plots. Here, the model includes two variance components, which are denoted as$\sigma^2_1$(for between-district heterogeneity) and$\sigma^2_2$(for between-school-within-district heterogeneity). Likelihood profiles for these two components can be obtained with: par(mfrow=c(2,1)) profile(res.ml, sigma2=1) profile(res.ml, sigma2=2) Both profile likelihood plots are peaked at the respective parameter estimates (as indicated by the vertical dotted lines) and the log likelihoods quickly decrease (i.e., become more negative) as the values of the components are moved away from the actual REML estimates. Hence, we can be fairly confident that both variance components are identifiable. ### Intraclass Correlation of the True Effects The three-level model used in the meta-analysis above allows for the underlying true effects within districts to be correlated. In fact, it can be shown that the model implies an intraclass correlation coefficient (ICC) of $$\rho = \frac{\sigma^2_1}{\sigma^2_1 + \sigma^2_2}$$ for true effects within the same level of the grouping variable (where$\sigma^2_1$is the variance component corresponding to the grouping variable and$\sigma^2_2$is the variance component corresponding to the level nested within the grouping variable). In the present case, the ICC is estimated to be: round(res.ml$sigma2[1] / sum(res.ml$sigma2), 3) [1] 0.665 Therefore, the underlying true effects within districts are estimated to correlate quite strongly (the simpler two-level model we fitted at the beginning ignores this dependence). Also, it is worth noting that the sum of the two variance components can be interpreted as the total amount of heterogeneity in the true effects: round(sum(res.ml$sigma2), 3)
[1] 0.098

### Multivariate Parameterization of the Model

We can fit the same model using a multivariate parameterization, which directly provides us with the ICC. For this, we use the syntax:

res.mv <- rma.mv(yi, vi, random = ~ study | district, data=dat)
print(res.mv, digits=3)
Multivariate Meta-Analysis Model (k = 56; method: REML)

Variance Components:

outer factor: district (nlvls = 11)
inner factor: study    (nlvls = 56)

estim   sqrt  fixed
tau^2      0.098  0.313     no
rho        0.665            no

Test for Heterogeneity:
Q(df = 55) = 578.864, p-val < .001

Model Results:

estimate     se   zval   pval  ci.lb  ci.ub
0.185  0.085  2.185  0.029  0.019  0.350  *

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

The random = ~ study | district argument adds correlated random effects for the different studies within districts to the model, where the variance-covariance matrix of the random effects takes on a compound symmetric structure (struct="CS" is the default). Note that the estimate of $\rho$ that is obtained is exactly the same as the ICC value we computed earlier based on the multilevel model. Also, the estimate of $\tau^2$ obtained from the multivariate parameterization is the same as the total amount of heterogeneity computed earlier based on the multilevel model. Note that random = ~ school | district would again yield the same results.

As long as $\rho$ is estimated to be positive, the multilevel and multivariate parametrizations are in essence identical. In fact, the log likelihoods of the two models should be identical, which we can confirm with:

logLik(res.ml)
logLik(res.mv)
'log Lik.' -7.958724 (df=3)
'log Lik.' -7.958724 (df=3)

### Profile Likelihood Plots

Again, profile likelihood plots should be used to ensure the identifiability of the variance (and correlation) components when fitting complex meta-analytic models. In the multivariate model, there is one variance component, denoted as $\tau^2$, and one correlation component, denoted as $\rho$. Profile likelihood plots for these two parameters can be obtained with:

par(mfrow=c(2,1))
profile(res.mv, tau2=1)
profile(res.mv, rho=1)

Again, both plots indicate that the estimates obtained in fact correspond to the peaks of the respective likelihood profiles, with decreasing log likelihoods as the values of parameters are moved away from the actual estimates.

Since the log likelihood drops of quite dramatically when $\rho$ is set equal to a value very close to 1, the left-hand side of the profile gets 'squished' together at the top and it is more difficult to see the curvature around the estimate. One can change the x-axis limits with the xlim argument. For example, with profile(res.mv, rho=1, xlim=c(0.3,0.9)) we can obtain a nicer profile plot for $\rho$.

### Uncorrelated Sampling Errors

It is important to note that the models used above assume that the sampling errors of the effect size estimates are independent. This is typically an appropriate assumption as long as there is no overlap in the data/individuals used to compute the various estimates. However, when multiple estimates are obtained from the same group of individuals, then this assumption is most certainly violated. Similarly, when multiple treatment groups are contrasted against a common control group, then the reuse of the control group data automatically induces a correlation between the estimates. See Berkey et al. (1998) and Gleser & Olkin (2009) for examples involving correlated sampling errors.

### A Common Mistake in the Three-Level Model

When fitting the three-level model, a common mistake is to forget to add random effects at both levels. In other words, one could fit the model:

rma.mv(yi, vi, random = ~ 1 | district, data=dat, digits=3)
Multivariate Meta-Analysis Model (k = 56; method: REML)

Variance Components:

estim   sqrt  nlvls  fixed    factor
sigma^2    0.083  0.288     11     no  district

Test for Heterogeneity:
Q(df = 55) = 578.864, p-val < .001

Model Results:

estimate     se   zval   pval  ci.lb  ci.ub
0.196  0.090  2.179  0.029  0.020  0.372  *

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

This model only accounts for heterogeneity between districts, so it assumes that the true effects for different studies within the same district are all identical (or put differently, that the ICC is equal to 1). In the present case, the ICC in the three-level model is actually quite high, so the results are not dramatically different, but this may not be true in other applications.

### Overparameterized Models

Finally, for illustration purposes, it is instructive to examine what can happen when we fit an overparameterized model. For example, suppose we combine the multilevel and multivariate structures above in a single model:

res.op <- rma.mv(yi, vi, random = list(~ study | district, ~ 1 | district, ~ 1 | study), data=dat)
print(res.op, digits=3)
Multivariate Meta-Analysis Model (k = 56; method: REML)

Variance Components:

estim   sqrt  nlvls  fixed    factor
sigma^2.1  0.041  0.203     11     no  district
sigma^2.2  0.026  0.162     56     no     study

outer factor: district (nlvls = 11)
inner factor: study    (nlvls = 56)

estim   sqrt  fixed
tau^2      0.030  0.174     no
rho        0.784            no

Test for Heterogeneity:
Q(df = 55) = 578.864, p-val < .001

Model Results:

estimate     se   zval   pval  ci.lb  ci.ub
0.185  0.085  2.185  0.029  0.019  0.350  *

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

In essence, since the model is overparameterized, the estimates of the variance and correlation components are largely arbitrary. However, it is not obvious from the results above that this is the case. By examining the profile likelihood plots of the variance and correlation components, we can determine that we are in fact in trouble:

par(mfrow=c(2,2))
profile(res.op)

(all variance components and correlations are automatically profiled when passing the model object to the profile() function.)

Note that large parts of the profiles are flat (and the profile for $\sigma^2_1$ shows some gaps, which indicates that the model did not converge for some of the profile points). These are clear indications that the model is overparameterized. In fact, despite the added complexity of the model, the fit of the model is exactly the same as for the multilevel and multivariate models fitted earlier:

logLik(res.op)
'log Lik.' -7.958724 (df=5)

In the present case, the estimate of the fixed effect and corresponding standard error (as shown under the Model Results) are also the same as the values obtained earlier. So, for making inferences about the fixed effects, things happen to work out okay in the overparameterized model. However, that may not be the case in other cases, so such models should be avoided.

### References

Cooper, H., Valentine, J. C., Charlton, K., & Melson, A. (2003). The effects of modified school calendars on student achievement and on school and community attitudes. Review of Educational Research, 73(1), 1–52.

Konstantopoulos, S. (2011). Fixed effects and variance components estimation in three-level meta-analysis. Research Synthesis Methods, 2(1), 61–76.