Wolfgang Viechtbauer

Maastricht University

Maastricht University

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.
This document provides the R code to reproduce the analyses conducted in
chapter 7 (which is focused on the computational aspects of a
meta-analysis, including the equal-, random-, and the mixed-effects
model) using the `metafor`

package. To read more about the
package, see the package
website and the package
documentation.

The package can be installed with:

Once the package is installed, we can load it with:

The data used for the analyses 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))
dat
```

```
## 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
```

The `yi`

values are standardized mean differences based on
studies examining the effectiveness of ‘challenge programs’ to treat
juvenile delinquency (presumably with positive values indicating that
delinquent behaviors were reduced in the group of juveniles
participating in such a challenge program compared to those who did
not). The `vi`

values are the corresponding sampling
variances. The dummy variable `random`

indicates whether
juveniles were randomly assigned to the two conditions with a study (1 =
yes, 0 = no). Finally, `intensity`

is a variable coded for
the purposes of the meta-analysis which indicates (based on assessment
of the coder) the intensity of the challenge program (coded from 1 =
‘very low’ to 7 = ‘very high’). These last two variables are potential
moderators of the treatment effectiveness.

We can fit an equal-effects model to these data with:

```
## Equal-Effects Model (k = 10)
##
## 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.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 for \(\tau^2\)) can be fitted with:

```
## 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.

A forest plot showing the results of the individual studies and the results from the random-effects model can be created as follows:

```
par(mar=c(4,4,2,2))
forest(res.re, xlim=c(-2,3), header=c("Study", "SMD [95% CI]"), top=2,
xlab="Standardized Mean Difference")
```

Next, we can carry out the ANOVA-type analysis described on pages
135-138. In particular, a fixed-effects meta-regression model with the
`random`

variable as moderator can be fitted with:

```
## Fixed-Effects with Moderators Model (k = 10)
##
## I^2 (residual heterogeneity / unaccounted variability): 0.00%
## H^2 (unaccounted variability / sampling variability): 0.96
## R^2 (amount of heterogeneity accounted for): 41.40%
##
## Test for Residual Heterogeneity:
## QE(df = 8) = 7.6901, p-val = 0.4643
##
## Test of Moderators (coefficient 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 equal-effects models (i.e., \(7.6901 + 7.0739 = 14.7640\)).

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

```
## Equal-Effects Model (k = 6)
##
## I^2 (total heterogeneity / total variability): 22.34%
## H^2 (total variability / sampling variability): 1.29
##
## 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
```

```
## Equal-Effects Model (k = 4)
##
## I^2 (total heterogeneity / total variability): 0.00%
## H^2 (total variability / sampling variability): 0.42
##
## 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:

```
## 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 meta-regression model with
moderators. Here, the authors are including both the `random`

and the `intensity`

moderators in the model. This can be done
with:

```
## Fixed-Effects with Moderators Model (k = 10)
##
## I^2 (residual heterogeneity / unaccounted variability): 8.89%
## H^2 (unaccounted variability / sampling variability): 1.10
## R^2 (amount of heterogeneity accounted for): 33.09%
##
## Test for Residual Heterogeneity:
## QE(df = 7) = 7.6832, p-val = 0.3614
##
## Test of Moderators (coefficients 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 meta-regression model is described on pages 140-142. These results can be obtained with:

```
## 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 (coefficients 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).

In essence, the model describes the (linear) relationship between the treatment intensity and the standardized mean differences and allows the intercept of the regression line to differ depending on whether random assignment was used within the studies or not. We can visualize the results as follows:

```
par(mar=c(5,4,2,2))
regplot(res.me, mod="intensity", xlim=c(1,7), pred=FALSE, ci=FALSE,
xlab="Intensity", ylab="Standardized Mean Difference",
pch=ifelse(dat$random == 0, 1, 19))
abline(a=coef(res.me)[1], b=coef(res.me)[3], lwd=2, lty="dotted")
abline(a=coef(res.me)[1] + coef(res.me)[2], b=coef(res.me)[3], lwd=2)
legend("topleft", inset=.02, lty=c("dotted", "solid"), lwd=2,
pch=c(21,19) ,legend=c("No Random Assignment", "Random Assignment"))
```

The points are drawn in size proportional to the weight the various studies received according to the model.

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 here). The value is given in the output above. We can also carry out this computation explicitly by doing a full versus reduced model comparison with (note that a warning will be issued with respect to the likelihood ratio test that is also part of this output since LRTs should be based on ML/REML estimation):

`## Warning: LRTs should be based on ML/REML estimation.`

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

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.

This documented is licensed under the following license: CC Attribution-Noncommercial-Share Alike 4.0 International.