The metafor Package

The book Practical meta-analysis by Lipsey and Wilson (2001) is an excellent introduction to meta-analysis and covers in less than 250 pages the most important aspects of conducting a meta-analysis. Here, I will reproduce the analyses conducted in chapter 7, which is focused on the computational aspects of a meta-analysis, including the fixed-, random-, and the mixed-effects model. The data are provided in the book on page 130 in Table 7.1. We can create the same dataset with:

dat <- data.frame(
id = c(100, 308, 1596, 2479, 9021, 9028, 161, 172, 537, 7049),
yi = c(-0.33, 0.32, 0.39, 0.31, 0.17, 0.64, -0.33, 0.15, -0.02, 0.00),
vi = c(0.084, 0.035, 0.017, 0.034, 0.072, 0.117, 0.102, 0.093, 0.012, 0.067),
random = c(0, 0, 0, 0, 0, 0, 1, 1, 1, 1),
intensity = c(7, 3, 7, 5, 7, 7, 4, 4, 5, 6))

The contents of dat are then:

     id    yi    vi random intensity
1   100 -0.33 0.084      0         7
2   308  0.32 0.035      0         3
3  1596  0.39 0.017      0         7
4  2479  0.31 0.034      0         5
5  9021  0.17 0.072      0         7
6  9028  0.64 0.117      0         7
7   161 -0.33 0.102      1         4
8   172  0.15 0.093      1         4
9   537 -0.02 0.012      1         5
10 7049  0.00 0.067      1         6

After loading the metafor package with library(metafor), we can fit a fixed-effects model to these data with:

res <- rma(yi, vi, data=dat, method="FE")
res

Fixed-Effects Model (k = 10)
 
Test for Heterogeneity: 
Q(df = 9) = 14.7640, p-val = 0.0976
 
Model Results:
 
estimate       se     zval     pval    ci.lb    ci.ub          
  0.1549   0.0609   2.5450   0.0109   0.0356   0.2742        * 
 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

These results match what is reported on page 132-133 (see Exhibit 7.3).

A random-effects model (using the DerSimonian-Laird estimator) can be fitted with:

res <- rma(yi, vi, data=dat, method="DL")
res

Random-Effects Model (k = 10; tau^2 estimator: DL)
 
tau^2 (estimated amount of total heterogeneity): 0.0260 (SE = 0.0328)
tau (square root of estimated tau^2 value):      0.1611
I^2 (total heterogeneity / total variability):   39.04%
H^2 (total variability / sampling variability):  1.64
 
Test for Heterogeneity: 
Q(df = 9) = 14.7640, p-val = 0.0976
 
Model Results:
 
estimate       se     zval     pval    ci.lb    ci.ub          
  0.1534   0.0858   1.7893   0.0736  -0.0146   0.3215        . 
 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Again, these results match what is shown in Exhibit 7.3 and on pages 134-135.

Next, we can carry out the ANOVA-type analysis described on pages 135-138. In particular, a fixed-effects model with the random variable as moderator can be fitted with:

res <- rma(yi, vi, mods = ~ random, data=dat, method="FE")
res

Fixed-Effects with Moderators Model (k = 10)
 
Test for Residual Heterogeneity: 
QE(df = 8) = 7.6901, p-val = 0.4643
 
Test of Moderators (coefficient(s) 2): 
QM(df = 1) = 7.0739, p-val = 0.0078
 
Model Results:
 
         estimate      se     zval    pval    ci.lb    ci.ub     
intrcpt    0.2984  0.0813   3.6687  0.0002   0.1390   0.4578  ***
random    -0.3261  0.1226  -2.6597  0.0078  -0.5664  -0.0858   **
 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Note that the $Q_E$ and $Q_M$ statistics add up to the $Q$ statistic given earlier for the fixed/random-effects models (i.e., $7.6901 + 7.0739 = 14.764$).

The tests for heterogeneity within the two subgroups formed by the moderator can be obtained with:

rma(yi, vi, data=dat, method="FE", subset=random==0)

Fixed-Effects Model (k = 6)
 
Test for Heterogeneity: 
Q(df = 5) = 6.4382, p-val = 0.2659
 
Model Results:
 
estimate       se     zval     pval    ci.lb    ci.ub          
  0.2984   0.0813   3.6687   0.0002   0.1390   0.4578      *** 
 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

and

rma(yi, vi, data=dat, method="FE", subset=random==1)

Fixed-Effects Model (k = 4)
 
Test for Heterogeneity: 
Q(df = 3) = 1.2519, p-val = 0.7406
 
Model Results:
 
estimate       se     zval     pval    ci.lb    ci.ub          
 -0.0277   0.0917  -0.3017   0.7628  -0.2075   0.1521          
 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Here, we note that the two separate $Q$ statistics add up to the $Q_E$ statistic given earlier (i.e., $6.4382 + 1.2519 = 7.6901$).

Finally, the predicted/estimated effect for random = 0 and random = 1 can be obtained with:

predict(res, newmods=c(0,1))

     pred     se   ci.lb  ci.ub
1  0.2984 0.0813  0.1390 0.4578
2 -0.0277 0.0917 -0.2075 0.1521

All of the values shown in Exhibit 7.4 have now been computed.

The "weighted regression analysis" described on pages 138-140 is again nothing else than a fixed-effects model with moderators. Here, the authors are including both the random and the intensity moderators in the model. This can be done with:

res <- rma(yi, vi, mods = ~ random + intensity, data=dat, method="FE")
res

Fixed-Effects with Moderators Model (k = 10)
 
Test for Residual Heterogeneity: 
QE(df = 7) = 7.6832, p-val = 0.3614
 
Test of Moderators (coefficient(s) 2,3): 
QM(df = 2) = 7.0807, p-val = 0.0290
 
Model Results:
 
           estimate      se     zval    pval    ci.lb    ci.ub   
intrcpt      0.3223  0.2998   1.0752  0.2823  -0.2652   0.9099   
random      -0.3298  0.1304  -2.5286  0.0115  -0.5854  -0.0742  *
intensity   -0.0041  0.0493  -0.0829  0.9339  -0.1007   0.0925   
 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Exhibit 7.6 (page 141) shows the same results.

Finally, the mixed-effects model is described on pages 140-142. These results can be obtained with:

res <- rma(yi, vi, mods = ~ random + intensity, data=dat, method="DL")
res

Mixed-Effects Model (k = 10; tau^2 estimator: DL)
 
tau^2 (estimated amount of residual heterogeneity):     0.0049 (SE = 0.0294)
tau (square root of estimated tau^2 value):             0.0698
I^2 (residual heterogeneity / unaccounted variability): 8.89%
H^2 (unaccounted variability / sampling variability):   1.10
R^2 (amount of heterogeneity accounted for):            81.20%
 
Test for Residual Heterogeneity: 
QE(df = 7) = 7.6832, p-val = 0.3614
 
Test of Moderators (coefficient(s) 2,3): 
QM(df = 2) = 5.5711, p-val = 0.0617
 
Model Results:
 
           estimate      se     zval    pval    ci.lb    ci.ub   
intrcpt      0.3311  0.3198   1.0351  0.3006  -0.2958   0.9579   
random      -0.3269  0.1439  -2.2712  0.0231  -0.6090  -0.0448  *
intensity   -0.0068  0.0528  -0.1292  0.8972  -0.1103   0.0967   
 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

These are the same results as shown in Exhibit 7.7 (page 141).

The amount of heterogeneity accounted for by the moderator(s) included in the model is typically computed by comparing the amount of heterogeneity from the random-effects model with the amount of residual heterogeneity from the mixed-effects model (for more details, see the FAQs). The value is given in the output above. We can also carry out this computation explicitly with:

res.ME <- rma(yi, vi, mods = ~ random + intensity, data=dat, method="DL")
res.RE <- rma(yi, vi, data=dat, method="DL")
anova(res.RE, res.ME)

        df    AIC    BIC  logLik    LRT   pval      QE  tau^2   R^2
Full     4 4.9040 6.1143  1.5480                7.6832 0.0049      
Reduced  2 6.1347 6.7399 -1.0674 5.2308 0.0731 14.7640 0.0260 81.2%

The value under R^2 is the amount of variance (heterogeneity) accounted for. In this particular example, we estimate that 81.2% of the total amount of heterogeneity is accounted for by the two moderators. Note that Lipsey and Wilson calculate $R^2$ in a different way than is usually done in meta-analyses, so the value for $R^2$ shown in Exhibit 7.7 (page 141) is different from the value computed above.

Lipsey, M. W., & Wilson, D. B. (2001). Practical meta-analysis. Thousand Oaks, CA: Sage.