R Code Corresponding to the Book Introduction to Meta-Analysis by Borenstein et al. (2009)
Wolfgang Viechtbauer
Maastricht University
Maastricht University
01 October, 2023
General Notes / Setup
As the title promises, the book Introduction to
Meta-Analysis by Borenstein et al. (2009) is a comprehensive
introduction to standard meta-analytic methodology, with focus on the
statistical methods. This document provides the R code to reproduce all
examples and illustrations from the book 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:
A few additional notes:
- Results are only reproduced for chapters containing worked examples.
- Occasionally, there are some minor discrepancies between the results shown in the book and those obtained below. Where such discrepancies arise, they are noted (and the reasons for them if they are known).
- I did not attempt to reproduce the figures exactly as they appear in
the book. There are some fundamental differences in the aesthetics of
the figures shown in the book and the functions in the
metaforpackage for producing the corresponding figures. I made some adjustments to the defaults here and there to make the figures below look similar to the ones shown in the book and made sure to include all relevant elements. - The results are generally given without discussion or context. The code below is not a substitute for reading the book, but is meant to be used together with it. In other words, readers of the book interested in replicating the results with R can see here how this is possible.
1) How a Meta-Analysis Works
# analysis of the Canon et al. (2006) data using (log) risk ratios
dat <- dat.cannon2006[,1:6] # keep only variables we really need
dat <- escalc(measure="RR", ai=ep1t, n1i=nt, ci=ep1c, n2i=nc, data=dat, slab=trial)
dat## trial pop nt nc ep1t ep1c yi vi
## 1 PROVE IT post-ACS 2099 2063 147 172 -0.1744 0.0117
## 2 A-TO-Z post-ACS 2265 2232 205 235 -0.1513 0.0082
## 3 TNT stable CAD 4995 5006 334 418 -0.2221 0.0050
## 4 IDEAL stable CAD 4439 4449 411 463 -0.1169 0.0041
## Random-Effects Model (k = 4; tau^2 estimator: DL)
##
## tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.0053)
## tau (square root of estimated tau^2 value): 0
## I^2 (total heterogeneity / total variability): 0.00%
## H^2 (total variability / sampling variability): 1.00
##
## Test for Heterogeneity:
## Q(df = 3) = 1.2425, p-val = 0.7428
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.1634 0.0393 -4.1635 <.0001 -0.2404 -0.0865 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## pred ci.lb ci.ub pi.lb pi.ub
## 0.85 0.79 0.92 0.79 0.92
dat$weights <- paste0(round(weights(res)), "%") # weights in % (rounded)
dat$pvals <- round(summary(dat)$pval, digits=3) # p-values of the individual trials# Figure 1.1
par(mar=c(4,4,2,2))
forest(res, xlim=c(-1,2), atransf=exp, at=log(c(2/3, 1, 3/2)),
header=TRUE, top=2, mlab="Summary", efac=c(0,1,3),
ilab=data.frame(dat$weights, dat$pvals), ilab.xpos=c(0.8,1.2), ilab.pos=2)
text(0.8, -1, "100%", pos=2)
text(1.2, -1, formatC(res$pval, format="f", digits=5), pos=2)
text(0.8, 6, "Weight", pos=2, font=2)
text(1.2, 6, "P-Value", pos=2, font=2)
2) Why Perform a Meta-Analysis
# analysis of the Lau et al. (1992) data using (log) risk ratios
dat <- dat.lau1992
dat <- escalc(measure="RR", ai=ai, n1i=n1i, ci=ci, n2i=n2i, data=dat, slab=trial)
dat## trial year ai n1i ci n2i yi vi
## 1 Fletcher 1959 1 12 4 11 -1.4733 1.0758
## 2 Dewar 1963 4 21 7 21 -0.5596 0.2976
## 3 European 1 1969 20 83 15 84 0.2997 0.0927
## 4 European 2 1971 69 373 94 357 -0.3530 0.0196
## 5 Heikinheimo 1971 22 219 17 207 0.2015 0.0949
## 6 Italian 1971 19 164 18 157 0.0104 0.0957
## 7 Australian 1 1973 26 264 32 253 -0.2502 0.0620
## 8 Frankfurt 2 1973 13 102 29 104 -0.7829 0.0920
## 9 NHLBI SMIT 1974 7 53 3 54 0.8660 0.4388
## 10 Frank 1975 6 55 6 53 -0.0370 0.2963
## 11 Valere 1975 11 49 9 42 0.0465 0.1578
## 12 Klein 1976 4 14 1 9 0.9445 1.0675
## 13 UK Collab 1976 38 302 40 293 -0.0815 0.0446
## 14 Austrian 1977 37 352 65 376 -0.4975 0.0369
## 15 Australian 2 1977 25 123 31 107 -0.3545 0.0548
## 16 Lasierra 1977 1 13 3 11 -1.2657 1.1655
## 17 N Ger Collab 1977 63 249 51 234 0.1492 0.0272
## 18 Witchitz 1977 5 32 5 26 -0.2076 0.3303
## 19 European 3 1979 18 156 30 159 -0.4918 0.0762
## 20 ISAM 1986 54 859 63 882 -0.1277 0.0321
## 21 GISSI-1 1986 628 5860 758 5852 -0.1895 0.0026
## 22 Olson 1986 1 28 2 24 -0.8473 1.4226
## 23 Baroffio 1986 0 29 6 30 -2.5322 2.0883
## 24 Schreiber 1986 1 19 3 19 -1.0986 1.2281
## 25 Cribier 1986 1 21 1 23 0.0910 1.9089
## 26 Sainsous 1986 3 49 6 49 -0.6931 0.4592
## 27 Durand 1987 3 35 4 29 -0.4757 0.5203
## 28 White 1987 2 107 12 112 -1.7461 0.5651
## 29 Bassand 1987 4 52 7 55 -0.5035 0.3554
## 30 Vlay 1988 1 13 2 12 -0.7732 1.3397
## 31 Kennedy 1988 12 191 17 177 -0.4244 0.1313
## 32 ISIS-2 1988 791 8592 1029 8595 -0.2627 0.0020
## 33 Wisenberg 1988 2 41 5 25 -1.4110 0.6356
## Random-Effects Model (k = 33; tau^2 estimator: DL)
##
## tau^2 (estimated amount of total heterogeneity): 0.0077 (SE = 0.0127)
## tau (square root of estimated tau^2 value): 0.0876
## I^2 (total heterogeneity / total variability): 16.87%
## H^2 (total variability / sampling variability): 1.20
##
## Test for Heterogeneity:
## Q(df = 32) = 38.4942, p-val = 0.1991
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.2312 0.0468 -4.9345 <.0001 -0.3230 -0.1394 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## pred ci.lb ci.ub pi.lb pi.ub
## 0.79 0.72 0.87 0.65 0.96
## [1] "0.0000008"
# Figure 2.1
par(mar=c(4,4,2,2))
forest(res, xlim=c(-10,9), atransf=exp, at=log(c(0.01, 0.1, 1, 10, 100)),
header=TRUE, top=2, ilab=year, ilab.xpos=-6, digits=3L, cex=0.8,
efac=c(0,1), fonts=c("sans", "inconsolata"))
# note: using a monospaced font for the annotations on the right, so they are
# fully aligned; if the "Inconsolata" font is not available, can also use "mono"
text(-6, 35, "Year", font=2, cex=0.8)
4) Effect Sizes Based on Means
# mean difference assuming sigma^2_1 = sigma^2_1
dat <- escalc("MD", m1i=103, m2i=100, sd1i=5.5, sd2i=4.5, n1i=50, n2i=50, vtype="HO")
summary(dat) # note: summary() so we can also see the standard error (sei)## yi vi sei zi pval ci.lb ci.ub
## 1 3.0000 1.0100 1.0050 2.9851 0.0028 1.0303 4.9697
# mean difference not assuming sigma^2_1 = sigma^2_1
dat <- escalc("MD", m1i=103, m2i=100, sd1i=5.5, sd2i=4.5, n1i=50, n2i=50)
summary(dat)## yi vi sei zi pval ci.lb ci.ub
## 1 3.0000 1.0100 1.0050 2.9851 0.0028 1.0303 4.9697
# mean change
# note: if the pre- and post-test means and SDs are not available, but we know
# the mean and SD of the change scores, then set m2i=0, sd2i=0, and ri=0
dat <- escalc("MC", m1i=5, m2i=0, sd1i=10, sd2i=0, ni=50, ri=0)
summary(dat)## yi vi sei zi pval ci.lb ci.ub
## 1 5.0000 2.0000 1.4142 3.5355 0.0004 2.2282 7.7718
## yi vi sei zi pval ci.lb ci.ub
## 1 5.0000 2.0000 1.4142 3.5355 0.0004 2.2282 7.7718
# standardized mean difference (Hedges' g)
dat <- escalc("SMD", m1i=103, m2i=100, sd1i=5.5, sd2i=4.5, n1i=50, n2i=50, vtype="LS2")
summary(dat)## yi vi sei zi pval ci.lb ci.ub
## 1 0.5924 0.0411 0.2028 2.9208 0.0035 0.1949 0.9900
# note: by default, the sampling variance is computed in a slightly different
# way in the book compared to the metafor package; by using vtype="LS2", the
# same equation as given in the book is used; but the difference is usually
# negligible (the results above differ slightly from those given in the book
# due to intermittent rounding in the book)# standardized mean change (using raw score standardization)
dat <- escalc("SMCR", m1i=103, m2i=100, sd1i=5.5/sqrt(2*(1-0.7)), ni=50, ri=0.7, vtype="LS2")
summary(dat)## yi vi sei zi pval ci.lb ci.ub
## 1 0.4160 0.0134 0.1156 3.5986 0.0003 0.1894 0.6426
# note: 5.5 is the SD of the change scores, which we can convert to the SD of
# the raw scores by dividing it by sqrt(2*(1-r)) (this assumes that the SD was
# the same at the pre- and post-test)
# note: by default, the sampling variance is computed in a slightly different
# way in the book compared to the metafor package; by using vtype="LS2", the
# same equation as given in the book is used; but the difference is usually
# negligible (the results above differ slightly from those given in the book
# because the equation given in the book actually contains a mistake)# response ratio (log transformed ratio of means)
dat <- escalc("ROM", m1i=61.515, m2i=51.015, sd1i=19.475, sd2i=19.475, n1i=10, n2i=10)
summary(dat)## yi vi sei zi pval ci.lb ci.ub
## 1 0.1872 0.0246 0.1568 1.1934 0.2327 -0.1202 0.4945
5) Effect Sizes Based on Binary Data
## yi vi sei zi pval ci.lb ci.ub
## 1 -0.6931 0.2800 0.5292 -1.3099 0.1902 -1.7303 0.3440
## yi vi sei zi pval ci.lb ci.ub
## 1 -0.7472 0.3216 0.5671 -1.3175 0.1877 -1.8588 0.3643
## yi vi sei zi pval ci.lb ci.ub
## 1 -0.0500 0.0014 0.0371 -1.3484 0.1775 -0.1227 0.0227
6) Effect Sizes Based on Correlations
## yi vi sei zi pval ci.lb ci.ub
## 1 0.5493 0.0103 0.1015 5.4100 <.0001 0.3503 0.7483
14) Worked Examples (Part 1)
Example for Continuous Data
# Table 14.1: Dataset 1
dat <- read.table(header=TRUE, text = "
study mean1 sd1 n1 mean2 sd2 n2
Carroll 94 22 60 92 20 60
Grant 98 21 65 92 22 65
Peck 98 28 40 88 26 40
Donat 94 19 200 82 17 200
Stewart 98 21 50 88 22 45
Young 96 21 85 92 22 85")
dat <- escalc("SMD", m1i=mean1, sd1i=sd1, n1i=n1, m2i=mean2, sd2i=sd2, n2i=n2,
slab=study, data=dat, vtype="LS2")
dat## study mean1 sd1 n1 mean2 sd2 n2 yi vi
## 1 Carroll 94 22 60 92 20 60 0.0945 0.0329
## 2 Grant 98 21 65 92 22 65 0.2774 0.0307
## 3 Peck 98 28 40 88 26 40 0.3665 0.0499
## 4 Donat 94 19 200 82 17 200 0.6644 0.0105
## 5 Stewart 98 21 50 88 22 45 0.4618 0.0427
## 6 Young 96 21 85 92 22 85 0.1852 0.0234
## Equal-Effects Model (k = 6)
##
## I^2 (total heterogeneity / total variability): 58.34%
## H^2 (total variability / sampling variability): 2.40
##
## Test for Heterogeneity:
## Q(df = 5) = 12.00, p-val = 0.03
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.41 0.06 6.47 <.01 0.29 0.54 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Note: The book uses the term ‘fixed-effect model’, but I will use the term ‘equal-effects model’ throughout this document. My reasons for using the latter term are explained here.
# Figure 14.1
par(mar=c(4,4,2,2))
dat$weights <- paste0(round(weights(res)), "%") # weights in % (rounded)
forest(res, xlim=c(-2,3), header=TRUE, top=2, mlab="Summary",
ilab=dat$weights, ilab.xpos=1.6, ilab.pos=2, efac=c(0,1,1.5))
text(1.6, -1, "100%", pos=2)
text(1.6, 8, "Weight", pos=2, font=2)
## Random-Effects Model (k = 6; tau^2 estimator: DL)
##
## tau^2 (estimated amount of total heterogeneity): 0.04 (SE = 0.04)
## tau (square root of estimated tau^2 value): 0.19
## I^2 (total heterogeneity / total variability): 58.34%
## H^2 (total variability / sampling variability): 2.40
##
## Test for Heterogeneity:
## Q(df = 5) = 12.00, p-val = 0.03
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.36 0.11 3.40 <.01 0.15 0.56 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Figure 14.2
par(mar=c(4,4,2,2))
dat$weights <- paste0(round(weights(res)), "%") # weights in % (rounded)
forest(res, xlim=c(-2,3), header=TRUE, top=2, mlab="Summary",
ilab=dat$weights, ilab.xpos=1.6, ilab.pos=2, efac=c(0,1,1.5))
text(1.6, -1, "100%", pos=2)
text(1.6, 8, "Weight", pos=2, font=2)
Example for Binary Data
# Table 14.4: Dataset 2
dat <- read.table(header=TRUE, text = "
study events1 nonevents1 n1 events2 nonevents2 n2
Saint 12 53 65 16 49 65
Kelly 8 32 40 10 30 40
Pilbeam 14 66 80 19 61 80
Lane 25 375 400 80 320 400
Wright 8 32 40 11 29 40
Day 16 49 65 18 47 65")
dat <- escalc("OR", ai=events1, n1i=n1, ci=events2, n2i=n2, slab=study, data=dat)
dat## study events1 nonevents1 n1 events2 nonevents2 n2 yi vi
## 1 Saint 12 53 65 16 49 65 -0.3662 0.1851
## 2 Kelly 8 32 40 10 30 40 -0.2877 0.2896
## 3 Pilbeam 14 66 80 19 61 80 -0.3842 0.1556
## 4 Lane 25 375 400 80 320 400 -1.3218 0.0583
## 5 Wright 8 32 40 11 29 40 -0.4169 0.2816
## 6 Day 16 49 65 18 47 65 -0.1595 0.1597
## Equal-Effects Model (k = 6)
##
## I^2 (total heterogeneity / total variability): 52.61%
## H^2 (total variability / sampling variability): 2.11
##
## Test for Heterogeneity:
## Q(df = 5) = 10.55, p-val = 0.06
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.72 0.15 -4.71 <.01 -1.03 -0.42 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## pred ci.lb ci.ub
## 0.48 0.36 0.66
# Figure 14.3
par(mar=c(4,4,2,2))
dat$weights <- paste0(round(weights(res)), "%") # weights in % (rounded)
forest(res, xlim=c(-3.5,5), header=TRUE, top=2, mlab="Summary",
atransf=exp, at=log(c(0.125, 0.25, 0.5, 1, 2, 4)), digits=c(2L,3L),
ilab=dat$weights, ilab.xpos=2.5, ilab.pos=2, efac=c(0,1,1.5))
text(2.5, -1, "100%", pos=2)
text(2.5, 8, "Weight", pos=2, font=2)
## Random-Effects Model (k = 6; tau^2 estimator: DL)
##
## tau^2 (estimated amount of total heterogeneity): 0.17 (SE = 0.21)
## tau (square root of estimated tau^2 value): 0.42
## I^2 (total heterogeneity / total variability): 52.61%
## H^2 (total variability / sampling variability): 2.11
##
## Test for Heterogeneity:
## Q(df = 5) = 10.55, p-val = 0.06
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.57 0.24 -2.37 0.02 -1.03 -0.10 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Figure 14.4
par(mar=c(4,4,2,2))
dat$weights <- paste0(round(weights(res)), "%") # weights in % (rounded)
forest(res, xlim=c(-3.5,5), header=TRUE, top=2, mlab="Summary",
atransf=exp, at=log(c(0.125, 0.25, 0.5, 1, 2, 4)), digits=c(2L,3L),
ilab=dat$weights, ilab.xpos=2.5, ilab.pos=2, efac=c(0,1,1.5))
text(2.5, -1, "100%", pos=2)
text(2.5, 8, "Weight", pos=2, font=2)
Example for Correlational Data
# Table 14.7: Dataset 3
dat <- read.table(header=TRUE, text = "
study correl n
Fonda 0.50 40
Newman 0.60 90
Grant 0.40 25
Granger 0.20 400
Milland 0.70 60
Finch 0.45 50")
dat <- escalc("ZCOR", ri=correl, ni=n, slab=study, data=dat)
dat## study correl n yi vi
## 1 Fonda 0.50 40 0.5493 0.0270
## 2 Newman 0.60 90 0.6931 0.0115
## 3 Grant 0.40 25 0.4236 0.0455
## 4 Granger 0.20 400 0.2027 0.0025
## 5 Milland 0.70 60 0.8673 0.0175
## 6 Finch 0.45 50 0.4847 0.0213
## Equal-Effects Model (k = 6)
##
## I^2 (total heterogeneity / total variability): 86.17%
## H^2 (total variability / sampling variability): 7.23
##
## Test for Heterogeneity:
## Q(df = 5) = 36.14, p-val < .01
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.38 0.04 9.54 <.01 0.30 0.45 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## pred ci.lb ci.ub
## 0.36 0.29 0.42
# Figure 14.5
par(mar=c(4,4,2,2))
dat$weights <- paste0(round(weights(res)), "%") # weights in % (rounded)
forest(res, xlim=c(-2,4), header=TRUE, top=2, mlab="Summary",
atransf=transf.ztor, at=transf.rtoz(c(-0.5, 0, 0.5, 0.85)),
ilab=dat$weights, ilab.xpos=2.2, ilab.pos=2, efac=c(0,1,1.5))
text(2.2, -1, "100%", pos=2)
text(2.2, 8, "Weight", pos=2, font=2)
## Random-Effects Model (k = 6; tau^2 estimator: DL)
##
## tau^2 (estimated amount of total heterogeneity): 0.08 (SE = 0.07)
## tau (square root of estimated tau^2 value): 0.29
## I^2 (total heterogeneity / total variability): 86.17%
## H^2 (total variability / sampling variability): 7.23
##
## Test for Heterogeneity:
## Q(df = 5) = 36.14, p-val < .01
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.53 0.13 4.10 <.01 0.28 0.79 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## pred ci.lb ci.ub pi.lb pi.ub
## 0.49 0.27 0.66 -0.08 0.82
# Figure 14.6
par(mar=c(4,4,2,2))
dat$weights <- paste0(round(weights(res)), "%") # weights in % (rounded)
forest(res, xlim=c(-2,4), header=TRUE, top=2, mlab="Summary",
atransf=transf.ztor, at=transf.rtoz(c(-0.5, 0, 0.5, 0.85)),
ilab=dat$weights, ilab.xpos=2.2, ilab.pos=2, efac=c(0,1,1.5))
text(2.2, -1, "100%", pos=2)
text(2.2, 8, "Weight", pos=2, font=2)
18) Worked Examples (Part 2)
Notes
- The prediction interval is computed in a slightly different way in
the book compared to how this is done (by default) in the
metaforpackage. For more details, see here. However, by using argumentpi.type="Riley", we can obtain the same results as those provided in the book. - The confidence intervals for \(\tau^2\) and \(I^2\) are also computed in a different way.
By default, the
metaforpackage uses the Q-profile method, which is an exact method (under the assumptions of the model), while the methods described in the book are based on large-sample approximations. However, we can still obtain the latter by using argumenttype="HT".
Example for Continuous Data
# Table 14.1: Dataset 1
dat <- read.table(header=TRUE, text = "
study mean1 sd1 n1 mean2 sd2 n2
Carroll 94 22 60 92 20 60
Grant 98 21 65 92 22 65
Peck 98 28 40 88 26 40
Donat 94 19 200 82 17 200
Stewart 98 21 50 88 22 45
Young 96 21 85 92 22 85")
dat <- escalc("SMD", m1i=mean1, sd1i=sd1, n1i=n1, m2i=mean2, sd2i=sd2, n2i=n2,
slab=study, data=dat, vtype="LS2")
dat## study mean1 sd1 n1 mean2 sd2 n2 yi vi
## 1 Carroll 94 22 60 92 20 60 0.0945 0.0329
## 2 Grant 98 21 65 92 22 65 0.2774 0.0307
## 3 Peck 98 28 40 88 26 40 0.3665 0.0499
## 4 Donat 94 19 200 82 17 200 0.6644 0.0105
## 5 Stewart 98 21 50 88 22 45 0.4618 0.0427
## 6 Young 96 21 85 92 22 85 0.1852 0.0234
## Random-Effects Model (k = 6; tau^2 estimator: DL)
##
## tau^2 (estimated amount of total heterogeneity): 0.0373 (SE = 0.0420)
## tau (square root of estimated tau^2 value): 0.1932
## I^2 (total heterogeneity / total variability): 58.34%
## H^2 (total variability / sampling variability): 2.40
##
## Test for Heterogeneity:
## Q(df = 5) = 12.0033, p-val = 0.0347
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.3582 0.1052 3.4038 0.0007 0.1520 0.5645 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## pred se ci.lb ci.ub pi.lb pi.ub
## 0.3582 0.1052 0.1520 0.5645 -0.2525 0.9690
## estimate ci.lb ci.ub
## tau^2 0.0373 0.0000 0.1312
## tau 0.1932 0.0000 0.3622
## I^2(%) 58.3447 0.0000 83.1242
## H^2 2.4007 1.0000 5.9256
# Figure 18.1
par(mar=c(4,4,2,2))
dat$weights <- paste0(round(weights(res)), "%") # weights in % (rounded)
forest(res, xlim=c(-2,3), header=TRUE, top=2, mlab="Summary",
ilab=dat$weights, ilab.xpos=1.6, ilab.pos=2, efac=c(0,1,1.5),
addpred=TRUE, pi.type="Riley")
text(1.6, -1, "100%", pos=2)
text(1.6, 8, "Weight", pos=2, font=2)
Example for Binary Data
# Table 14.4: Dataset 2
dat <- read.table(header=TRUE, text = "
study events1 nonevents1 n1 events2 nonevents2 n2
Saint 12 53 65 16 49 65
Kelly 8 32 40 10 30 40
Pilbeam 14 66 80 19 61 80
Lane 25 375 400 80 320 400
Wright 8 32 40 11 29 40
Day 16 49 65 18 47 65")
dat <- escalc("OR", ai=events1, n1i=n1, ci=events2, n2i=n2, slab=study, data=dat)
dat## study events1 nonevents1 n1 events2 nonevents2 n2 yi vi
## 1 Saint 12 53 65 16 49 65 -0.3662 0.1851
## 2 Kelly 8 32 40 10 30 40 -0.2877 0.2896
## 3 Pilbeam 14 66 80 19 61 80 -0.3842 0.1556
## 4 Lane 25 375 400 80 320 400 -1.3218 0.0583
## 5 Wright 8 32 40 11 29 40 -0.4169 0.2816
## 6 Day 16 49 65 18 47 65 -0.1595 0.1597
## Random-Effects Model (k = 6; tau^2 estimator: DL)
##
## tau^2 (estimated amount of total heterogeneity): 0.1729 (SE = 0.2148)
## tau (square root of estimated tau^2 value): 0.4158
## I^2 (total heterogeneity / total variability): 52.61%
## H^2 (total variability / sampling variability): 2.11
##
## Test for Heterogeneity:
## Q(df = 5) = 10.5512, p-val = 0.0610
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.5663 0.2388 -2.3711 0.0177 -1.0344 -0.0982 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## pred ci.lb ci.ub pi.lb pi.ub
## 0.5676 0.3554 0.9065 0.1499 2.1492
## estimate ci.lb ci.ub
## tau^2 0.1729 0.0000 0.6676
## tau 0.4158 0.0000 0.8171
## I^2(%) 52.6118 0.0000 81.0849
## H^2 2.1102 1.0000 5.2868
# Figure 18.2
par(mar=c(4,4,2,2))
dat$weights <- paste0(round(weights(res)), "%") # weights in % (rounded)
forest(res, xlim=c(-3.5,5), header=TRUE, top=2, mlab="Summary",
atransf=exp, at=log(c(0.125, 0.25, 0.5, 1, 2, 4)), digits=c(2L,3L),
ilab=dat$weights, ilab.xpos=2.5, ilab.pos=2, efac=c(0,1,1.5),
addpred=TRUE, pi.type="Riley")
text(2.5, -1, "100%", pos=2)
text(2.5, 8, "Weight", pos=2, font=2)
Example for Correlational Data
# Table 14.7: Dataset 3
dat <- read.table(header=TRUE, text = "
study correl n
Fonda 0.50 40
Newman 0.60 90
Grant 0.40 25
Granger 0.20 400
Milland 0.70 60
Finch 0.45 50")
dat <- escalc("ZCOR", ri=correl, ni=n, slab=study, data=dat)
dat## study correl n yi vi
## 1 Fonda 0.50 40 0.5493 0.0270
## 2 Newman 0.60 90 0.6931 0.0115
## 3 Grant 0.40 25 0.4236 0.0455
## 4 Granger 0.20 400 0.2027 0.0025
## 5 Milland 0.70 60 0.8673 0.0175
## 6 Finch 0.45 50 0.4847 0.0213
## Random-Effects Model (k = 6; tau^2 estimator: DL)
##
## tau^2 (estimated amount of total heterogeneity): 0.0819 (SE = 0.0727)
## tau (square root of estimated tau^2 value): 0.2861
## I^2 (total heterogeneity / total variability): 86.17%
## H^2 (total variability / sampling variability): 7.23
##
## Test for Heterogeneity:
## Q(df = 5) = 36.1437, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.5328 0.1298 4.1045 <.0001 0.2784 0.7872 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## pred ci.lb ci.ub pi.lb pi.ub
## 0.4875 0.2714 0.6568 -0.3271 0.8865
## estimate ci.lb ci.ub
## tau^2 0.0819 0.0338 0.1791
## tau 0.2861 0.1839 0.4232
## I^2(%) 86.1663 72.0135 93.1620
## H^2 7.2287 3.5732 14.6242
# Figure 18.3
par(mar=c(4,4,2,2))
dat$weights <- paste0(round(weights(res)), "%") # weights in % (rounded)
forest(res, xlim=c(-2,4), header=TRUE, top=2, mlab="Summary",
atransf=transf.ztor, at=transf.rtoz(c(-0.5, 0, 0.5, 0.9)),
ilab=dat$weights, ilab.xpos=2.2, ilab.pos=2, efac=c(0,1,1.5),
addpred=TRUE, pi.type="Riley")
text(2.2, -1, "100%", pos=2)
text(2.2, 8, "Weight", pos=2, font=2)
19) Subgroup Analyses
dat <- read.table(header=TRUE, text = "
study type g var
Thornhill A 0.110 0.010
Kendall A 0.224 0.030
Vandamm A 0.338 0.020
Leonard A 0.451 0.015
Professor A 0.480 0.010
Jeffries B 0.440 0.015
Fremont B 0.492 0.020
Doyle B 0.651 0.015
Stella B 0.710 0.025
Thorwald B 0.740 0.012")
dat## study type g var
## 1 Thornhill A 0.110 0.010
## 2 Kendall A 0.224 0.030
## 3 Vandamm A 0.338 0.020
## 4 Leonard A 0.451 0.015
## 5 Professor A 0.480 0.010
## 6 Jeffries B 0.440 0.015
## 7 Fremont B 0.492 0.020
## 8 Doyle B 0.651 0.015
## 9 Stella B 0.710 0.025
## 10 Thorwald B 0.740 0.012
EE models within subgroups
# equal-effects model for subgroup A
resA <- rma(g, var, data=dat, method="EE", subset=type=="A")
resA## Equal-Effects Model (k = 5)
##
## I^2 (total heterogeneity / total variability): 52.56%
## H^2 (total variability / sampling variability): 2.11
##
## Test for Heterogeneity:
## Q(df = 4) = 8.4316, p-val = 0.0770
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.3241 0.0535 6.0633 <.0001 0.2193 0.4289 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# equal-effects model for subgroup B
resB <- rma(g, var, data=dat, method="EE", subset=type=="B")
resB## Equal-Effects Model (k = 5)
##
## I^2 (total heterogeneity / total variability): 11.95%
## H^2 (total variability / sampling variability): 1.14
##
## Test for Heterogeneity:
## Q(df = 4) = 4.5429, p-val = 0.3375
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.6111 0.0571 10.7013 <.0001 0.4992 0.7230 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Equal-Effects Model (k = 10)
##
## I^2 (total heterogeneity / total variability): 65.96%
## H^2 (total variability / sampling variability): 2.94
##
## Test for Heterogeneity:
## Q(df = 9) = 26.4371, p-val = 0.0017
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.4581 0.0390 11.7396 <.0001 0.3816 0.5346 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Figure 19.1
par(mar=c(4,4,2,2))
forest(res, xlim=c(-1,2), cex=1, refline=0.5,
header=c("Type / Study", "SMD [95% CI]"),
xlab="Standardized Mean Difference", mlab="Combined",
at=seq(-0.2, 1.2, by=0.2), efac=c(0,1,1.5),
rows=c(3:7, 12:16), ylim=c(-1.2,21))
text(-1, 17.5, "Type A", pos=4, font=2)
text(-1, 8.5, "Type B", pos=4, font=2)
addpoly(resA, row=10.5, mlab="Combined")
addpoly(resB, row= 1.5, mlab="Combined")
# Table 19.2
sav <- list(resA, resB, res)
tab <- rbind(sapply(sav, function(x) coef(x)),
sapply(sav, function(x) vcov(x)),
sapply(sav, function(x) x$se),
sapply(sav, function(x) x$ci.lb),
sapply(sav, function(x) x$ci.ub),
sapply(sav, function(x) x$zval),
sapply(sav, function(x) x$pval),
sapply(sav, function(x) x$QE),
sapply(sav, function(x) x$k-1),
sapply(sav, function(x) x$QEp),
sapply(sav, function(x) x$I2))
colnames(tab) <- c("A", "B", "Combined")
rownames(tab) <- c("Y", "V", "SE", "LL", "UB", "Z", "pval", "Q", "df", "pval", "I2")
round(tab, digits=4)## A B Combined
## Y 0.3241 0.6111 0.4581
## V 0.0029 0.0033 0.0015
## SE 0.0535 0.0571 0.0390
## LL 0.2193 0.4992 0.3816
## UB 0.4289 0.7230 0.5346
## Z 6.0633 10.7013 11.7396
## pval 0.0000 0.0000 0.0000
## Q 8.4316 4.5429 26.4371
## df 4.0000 4.0000 9.0000
## pval 0.0770 0.3375 0.0017
## I2 52.5594 11.9506 65.9569
## Fixed-Effects with Moderators Model (k = 10)
##
## I^2 (residual heterogeneity / unaccounted variability): 38.34%
## H^2 (unaccounted variability / sampling variability): 1.62
## R^2 (amount of heterogeneity accounted for): 44.79%
##
## Test for Residual Heterogeneity:
## QE(df = 8) = 12.9745, p-val = 0.1127
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 13.4626, p-val = 0.0002
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.3241 0.0535 6.0633 <.0001 0.2193 0.4289 ***
## typeB 0.2870 0.0782 3.6691 0.0002 0.1337 0.4403 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# note: the estimate for 'typeB' is the estimated difference; the values
# for zval and pval are the test of the difference between A and B# Table 19.3
tab <- data.frame(Q = c(resA$QE, resB$QE, dif$QE, dif$QM, res$QE),
df = c(resA$k-1, resB$k-1, dif$k-dif$p, dif$m, res$k-1),
p = c(resA$QEp, resB$QEp, dif$QEp, dif$QMp, res$QEp))
rownames(tab) <- c("A", "B", "Within", "Between", "Total")
round(tab, digits=4)## Q df p
## A 8.4316 4 0.0770
## B 4.5429 4 0.3375
## Within 12.9745 8 0.1127
## Between 13.4626 1 0.0002
## Total 26.4371 9 0.0017
RE models with separate \(\tau^2\) estimates
# random-effect model within subgroup A
resA <- rma(g, var, data=dat, method="DL", subset=type=="A")
resA## Random-Effects Model (k = 5; tau^2 estimator: DL)
##
## tau^2 (estimated amount of total heterogeneity): 0.0164 (SE = 0.0225)
## tau (square root of estimated tau^2 value): 0.1282
## I^2 (total heterogeneity / total variability): 52.56%
## H^2 (total variability / sampling variability): 2.11
##
## Test for Heterogeneity:
## Q(df = 4) = 8.4316, p-val = 0.0770
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.3245 0.0799 4.0595 <.0001 0.1678 0.4812 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# random-effect model within subgroup B
resB <- rma(g, var, data=dat, method="DL", subset=type=="B")
resB## Random-Effects Model (k = 5; tau^2 estimator: DL)
##
## tau^2 (estimated amount of total heterogeneity): 0.0022 (SE = 0.0133)
## tau (square root of estimated tau^2 value): 0.0474
## I^2 (total heterogeneity / total variability): 11.95%
## H^2 (total variability / sampling variability): 1.14
##
## Test for Heterogeneity:
## Q(df = 4) = 4.5429, p-val = 0.3375
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.6100 0.0611 9.9833 <.0001 0.4903 0.7298 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# random-effect model for all 10 studies allowing for different tau^2 values within subgroups
# note: rma.mv() uses ML/REML estimation; here, we pass the tau^2 values as
# estimated by the DL method from the two subgroup models to the function
res <- rma.mv(g, var, random = ~ type | study, struct="DIAG",
data=dat, slab=study, tau2=c(resA$tau2, resB$tau2))
res## Multivariate Meta-Analysis Model (k = 10; method: REML)
##
## Variance Components:
##
## outer factor: study (nlvls = 10)
## inner factor: type (nlvls = 2)
##
## estim sqrt k.lvl fixed level
## tau^2.1 0.0164 0.1282 5 yes A
## tau^2.2 0.0022 0.0474 5 yes B
##
## Test for Heterogeneity:
## Q(df = 9) = 26.4371, p-val = 0.0017
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.5047 0.0485 10.3967 <.0001 0.4096 0.5999 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Figure 19.5
par(mar=c(4,4,2,2))
forest(res, xlim=c(-1,2), cex=1, refline=0.5,
header=c("Type / Study", "SMD [95% CI]"),
xlab="Standardized Mean Difference", mlab="Combined",
at=seq(-0.2, 1.2, by=0.2), efac=c(0,1,1.5),
rows=c(3:7, 12:16), ylim=c(-1.2,21))
text(-1, 17.5, "Type A", pos=4, font=2)
text(-1, 8.5, "Type B", pos=4, font=2)
addpoly(resA, row=10.5, mlab="Combined")
addpoly(resB, row= 1.5, mlab="Combined")
# Table 19.6
sav <- list(resA, resB, res)
tab <- rbind(sapply(sav, function(x) coef(x)),
sapply(sav, function(x) vcov(x)),
sapply(sav, function(x) x$se),
sapply(sav, function(x) x$ci.lb),
sapply(sav, function(x) x$ci.ub),
sapply(sav, function(x) x$zval),
sapply(sav, function(x) x$pval))
colnames(tab) <- c("A", "B", "Combined")
rownames(tab) <- c("Y", "V", "SE", "LL", "UB", "Z", "pval")
round(tab, digits=4)## A B Combined
## Y 0.3245 0.6100 0.5047
## V 0.0064 0.0037 0.0024
## SE 0.0799 0.0611 0.0485
## LL 0.1678 0.4903 0.4096
## UB 0.4812 0.7298 0.5999
## Z 4.0595 9.9833 10.3967
## pval 0.0000 0.0000 0.0000
# note: the Q-statistics in the table shown in the book are not traditional Q-test
# statistics, but incorporate the random-effects model weights (not shown here)# compare effects of A and B
dif <- rma.mv(g, var, mods = ~ type, random = ~ type | study, struct="DIAG",
data=dat, slab=study, tau2=c(resA$tau2, resB$tau2))
dif## Multivariate Meta-Analysis Model (k = 10; method: REML)
##
## Variance Components:
##
## outer factor: study (nlvls = 10)
## inner factor: type (nlvls = 2)
##
## estim sqrt k.lvl fixed level
## tau^2.1 0.0164 0.1282 5 yes A
## tau^2.2 0.0022 0.0474 5 yes B
##
## Test for Residual Heterogeneity:
## QE(df = 8) = 12.9745, p-val = 0.1127
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 8.0547, p-val = 0.0045
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.3245 0.0799 4.0595 <.0001 0.1678 0.4812 ***
## typeB 0.2856 0.1006 2.8381 0.0045 0.0884 0.4828 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# note: the estimate for 'typeB' is the estimated difference; the values
# for zval and pval are the test of the difference between A and B
# note: for the decomposition of the Q-statistics to work as shown in Table
# 19.7, one must compute these statistics incorporating the random-effects model
# weights; this gets confusing really quickly, because that's not how the
# Q-statistics for heterogeneity are typically computed, so I'll skip this part;
# however, note that the QM-test in the output above is the Q_Between test shown
# in the table and this is what we are interested in anywayRE models with pooled \(\tau^2\) estimate
# mixed-effect meta-regression model assuming a common tau^2 value for the 2 subgroups
dif <- rma(g, var, mods = ~ type - 1, data=dat, method="DL", slab=study)
dif## Mixed-Effects Model (k = 10; tau^2 estimator: DL)
##
## tau^2 (estimated amount of residual heterogeneity): 0.0097 (SE = 0.0128)
## tau (square root of estimated tau^2 value): 0.0986
## I^2 (residual heterogeneity / unaccounted variability): 38.34%
## H^2 (unaccounted variability / sampling variability): 1.62
##
## Test for Residual Heterogeneity:
## QE(df = 8) = 12.9745, p-val = 0.1127
##
## Test of Moderators (coefficients 1:2):
## QM(df = 2) = 91.2268, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## typeA 0.3247 0.0706 4.6001 <.0001 0.1864 0.4631 ***
## typeB 0.6083 0.0727 8.3705 <.0001 0.4659 0.7507 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# random-effect model for all 10 studies setting tau^2 to the value from the model above
res <- rma(g, var, data=dat, tau2=dif$tau2, slab=study)
res## Random-Effects Model (k = 10; user-specified tau^2 value)
##
## tau^2 (specified amount of total heterogeneity): 0.0097 (SE = 0.0119)
## tau (square root of specified tau^2 value): 0.0986
## I^2 (total heterogeneity / total variability): 38.67%
## H^2 (total variability / sampling variability): 1.63
##
## Test for Heterogeneity:
## Q(df = 9) = 26.4371, p-val = 0.0017
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.4624 0.0506 9.1321 <.0001 0.3632 0.5617 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Figure 19.8
par(mar=c(4,4,2,2))
forest(res, xlim=c(-1,2), cex=1, refline=0.5,
header=c("Type / Study", "SMD [95% CI]"),
xlab="Standardized Mean Difference", mlab="Combined",
at=seq(-0.2, 1.2, by=0.2), efac=c(0,1,1.5),
rows=c(3:7, 12:16), ylim=c(-1.2,21))
text(-1, 17.5, "Type A", pos=4, font=2)
text(-1, 8.5, "Type B", pos=4, font=2)
addpoly(coef(dif)[1], vcov(dif)[1,1], row=10.5, mlab="Combined")
addpoly(coef(dif)[2], vcov(dif)[2,2], row= 1.5, mlab="Combined")
# Table 19.11
tab <- rbind(c(coef(dif)[1], coef(dif)[2], coef(res)),
c(vcov(dif)[1,1], vcov(dif)[2,2], vcov(res)),
c(dif$se[1], dif$se[2], res$se),
c(dif$ci.lb[1], dif$ci.lb[2], res$ci.lb),
c(dif$ci.ub[1], dif$ci.ub[2], res$ci.ub),
c(dif$zval[1], dif$zval[2], res$zval),
c(dif$pval[1], dif$pval[2], res$pval))
colnames(tab) <- c("A", "B", "Combined")
rownames(tab) <- c("Y", "V", "SE", "LL", "UB", "Z", "pval")
round(tab, digits=4)## A B Combined
## Y 0.3247 0.6083 0.4624
## V 0.0050 0.0053 0.0026
## SE 0.0706 0.0727 0.0506
## LL 0.1864 0.4659 0.3632
## UB 0.4631 0.7507 0.5617
## Z 4.6001 8.3705 9.1321
## pval 0.0000 0.0000 0.0000
# note: the Q-statistics in the table shown in the book are not traditional Q-test
# statistics, but incorporate the random-effects model weights (not shown here)# compare effects of A and B
dif <- rma(g, var, mods = ~ type, data=dat, method="DL", slab=study)
dif## Mixed-Effects Model (k = 10; tau^2 estimator: DL)
##
## tau^2 (estimated amount of residual heterogeneity): 0.0097 (SE = 0.0128)
## tau (square root of estimated tau^2 value): 0.0986
## I^2 (residual heterogeneity / unaccounted variability): 38.34%
## H^2 (unaccounted variability / sampling variability): 1.62
## R^2 (amount of heterogeneity accounted for): 67.45%
##
## Test for Residual Heterogeneity:
## QE(df = 8) = 12.9745, p-val = 0.1127
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 7.8324, p-val = 0.0051
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.3247 0.0706 4.6001 <.0001 0.1864 0.4631 ***
## typeB 0.2835 0.1013 2.7987 0.0051 0.0850 0.4821 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# note: the estimate for 'typeB' is the estimated difference; the values
# for zval and pval are the test of the difference between A and B20) Meta-Regression
# compute log risk ratios for the BCG vaccine data
dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg,
slab=paste(author, year, sep=", "))
dat## trial author year tpos tneg cpos cneg ablat alloc yi vi
## 1 1 Aronson 1948 4 119 11 128 44 random -0.8893 0.3256
## 2 2 Ferguson & Simes 1949 6 300 29 274 55 random -1.5854 0.1946
## 3 3 Rosenthal et al 1960 3 228 11 209 42 random -1.3481 0.4154
## 4 4 Hart & Sutherland 1977 62 13536 248 12619 52 random -1.4416 0.0200
## 5 5 Frimodt-Moller et al 1973 33 5036 47 5761 13 alternate -0.2175 0.0512
## 6 6 Stein & Aronson 1953 180 1361 372 1079 44 alternate -0.7861 0.0069
## 7 7 Vandiviere et al 1973 8 2537 10 619 19 random -1.6209 0.2230
## 8 8 TPT Madras 1980 505 87886 499 87892 13 random 0.0120 0.0040
## 9 9 Coetzee & Berjak 1968 29 7470 45 7232 27 random -0.4694 0.0564
## 10 10 Rosenthal et al 1961 17 1699 65 1600 42 systematic -1.3713 0.0730
## 11 11 Comstock et al 1974 186 50448 141 27197 18 systematic -0.3394 0.0124
## 12 12 Comstock & Webster 1969 5 2493 3 2338 33 systematic 0.4459 0.5325
## 13 13 Comstock et al 1976 27 16886 29 17825 33 systematic -0.0173 0.0714
## Equal-Effects Model (k = 13)
##
## I^2 (total heterogeneity / total variability): 92.12%
## H^2 (total variability / sampling variability): 12.69
##
## Test for Heterogeneity:
## Q(df = 12) = 152.2330, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.4303 0.0405 -10.6247 <.0001 -0.5097 -0.3509 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Figure 20.1
par(mar=c(4,4,2,2))
forest(res, xlim=c(-8,7), atransf=exp, at=log(c(0.1, 1, 10)), digits=3L, cex=0.8,
header=TRUE, top=2, mlab="Summary", efac=c(0,1,1.5), order="obs")
# meta-regression model with absolute latitude
reg <- rma(yi, vi, mods = ~ ablat, data=dat, method="FE", digits=5)
reg## Fixed-Effects with Moderators Model (k = 13)
##
## I^2 (residual heterogeneity / unaccounted variability): 64.21%
## H^2 (unaccounted variability / sampling variability): 2.79
## R^2 (amount of heterogeneity accounted for): 77.98%
##
## Test for Residual Heterogeneity:
## QE(df = 11) = 30.73309, p-val = 0.00121
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 121.49992, p-val < .00001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.34356 0.08105 4.23899 0.00002 0.18471 0.50242 ***
## ablat -0.02924 0.00265 -11.02270 <.00001 -0.03444 -0.02404 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Table 20.3
tab <- data.frame(Q = c(reg$QM, reg$QE, res$QE),
df = c(reg$m, reg$k-reg$m, res$k-1),
pval = c(reg$QMp, reg$QEp, res$QEp))
rownames(tab) <- c("Model (Q_model)", "Residual (Q_resid)", "Total (Q)")
round(tab, digits=5)## Q df pval
## Model (Q_model) 121.49992 1 0.00000
## Residual (Q_resid) 30.73309 12 0.00121
## Total (Q) 152.23301 12 0.00000
# Figure 20.2
regplot(reg, xlim=c(0,70), ylim=c(-3,1), xlab="Latitude", ylab="ln(RR)", las=1,
pch=1, plim=c(NA,8), ci=FALSE)
title("Regression of log risk ratio on latitude (fixed-effect)")
## Random-Effects Model (k = 13; tau^2 estimator: DL)
##
## tau^2 (estimated amount of total heterogeneity): 0.3088 (SE = 0.2299)
## tau (square root of estimated tau^2 value): 0.5557
## I^2 (total heterogeneity / total variability): 92.12%
## H^2 (total variability / sampling variability): 12.69
##
## Test for Heterogeneity:
## Q(df = 12) = 152.2330, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.7141 0.1787 -3.9952 <.0001 -1.0644 -0.3638 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## pred ci.lb ci.ub pi.lb pi.ub
## 0.49 0.34 0.70 0.16 1.54
# Figure 20.5
par(mar=c(4,4,2,2))
forest(res, xlim=c(-8,7), atransf=exp, at=log(c(0.1, 1, 10)), digits=3L, cex=0.8,
header=TRUE, top=2, mlab="Summary", efac=c(0,1,1.5), order="obs")
# mixed-effects meta-regression model with absolute latitude
reg <- rma(yi, vi, mods = ~ ablat, data=dat, method="DL", digits=5)
reg## Mixed-Effects Model (k = 13; tau^2 estimator: DL)
##
## tau^2 (estimated amount of residual heterogeneity): 0.06330 (SE = 0.05483)
## tau (square root of estimated tau^2 value): 0.25160
## I^2 (residual heterogeneity / unaccounted variability): 64.21%
## H^2 (unaccounted variability / sampling variability): 2.79
## R^2 (amount of heterogeneity accounted for): 79.50%
##
## Test for Residual Heterogeneity:
## QE(df = 11) = 30.73309, p-val = 0.00121
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 18.84523, p-val = 0.00001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.25954 0.23231 1.11724 0.26389 -0.19577 0.71486
## ablat -0.02923 0.00673 -4.34111 0.00001 -0.04243 -0.01603 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Figure 20.6
regplot(reg, xlim=c(0,70), ylim=c(-3,1), xlab="Latitude", ylab="ln(RR)", las=1,
pch=1, plim=c(NA,8), ci=FALSE)
title("Regression of log risk ratio on latitude (random-effect)")
21) Subgrouping / Meta-Regression
The following results are not actually shown in the chapter. They are
provided here to show how these types of analyses can be conducted using
the metafor package.
## trial author year tpos tneg cpos cneg ablat alloc yi vi
## 1 1 Aronson 1948 4 119 11 128 44 random -0.8893 0.3256
## 2 2 Ferguson & Simes 1949 6 300 29 274 55 random -1.5854 0.1946
## 3 3 Rosenthal et al 1960 3 228 11 209 42 random -1.3481 0.4154
## 4 4 Hart & Sutherland 1977 62 13536 248 12619 52 random -1.4416 0.0200
## 5 5 Frimodt-Moller et al 1973 33 5036 47 5761 13 alternate -0.2175 0.0512
## 6 6 Stein & Aronson 1953 180 1361 372 1079 44 alternate -0.7861 0.0069
## 7 7 Vandiviere et al 1973 8 2537 10 619 19 random -1.6209 0.2230
## 8 8 TPT Madras 1980 505 87886 499 87892 13 random 0.0120 0.0040
## 9 9 Coetzee & Berjak 1968 29 7470 45 7232 27 random -0.4694 0.0564
## 10 10 Rosenthal et al 1961 17 1699 65 1600 42 systematic -1.3713 0.0730
## 11 11 Comstock et al 1974 186 50448 141 27197 18 systematic -0.3394 0.0124
## 12 12 Comstock & Webster 1969 5 2493 3 2338 33 systematic 0.4459 0.5325
## 13 13 Comstock et al 1976 27 16886 29 17825 33 systematic -0.0173 0.0714
# random-effects model with Knapp-Hartung method
res <- rma(yi, vi, data=dat, method="DL", test="knha")
res## Random-Effects Model (k = 13; tau^2 estimator: DL)
##
## tau^2 (estimated amount of total heterogeneity): 0.3088 (SE = 0.2299)
## tau (square root of estimated tau^2 value): 0.5557
## I^2 (total heterogeneity / total variability): 92.12%
## H^2 (total variability / sampling variability): 12.69
##
## Test for Heterogeneity:
## Q(df = 12) = 152.2330, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## -0.7141 0.1807 -3.9520 12 0.0019 -1.1078 -0.3204 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# mixed-effects meta-regression model with Knapp-Hartung method
res <- rma(yi, vi, mods = ~ ablat, data=dat, method="DL", test="knha")
res## Mixed-Effects Model (k = 13; tau^2 estimator: DL)
##
## tau^2 (estimated amount of residual heterogeneity): 0.0633 (SE = 0.0548)
## tau (square root of estimated tau^2 value): 0.2516
## I^2 (residual heterogeneity / unaccounted variability): 64.21%
## H^2 (unaccounted variability / sampling variability): 2.79
## R^2 (amount of heterogeneity accounted for): 79.50%
##
## Test for Residual Heterogeneity:
## QE(df = 11) = 30.7331, p-val = 0.0012
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 11) = 13.5589, p-val = 0.0036
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## intrcpt 0.2595 0.2739 0.9477 11 0.3636 -0.3433 0.8623
## ablat -0.0292 0.0079 -3.6822 11 0.0036 -0.0467 -0.0118 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# random-effects model with permutation test
res <- rma(yi, vi, data=dat, method="DL")
permutest(res, exact=TRUE)## Running 8192 iterations for an exact permutation test.
##
## Model Results:
##
## estimate se zval pval¹ ci.lb ci.ub
## -0.7141 0.1787 -3.9952 0.0017 -1.0644 -0.3638 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## 1) p-value based on permutation testing
# mixed-effects meta-regression model with permutation test
# note: an exact permutation test would require 389188800 iterations, which will
# take forever; instead, we can use a large number of random permutations; for
# reproducibility, we set the 'seed' of the random number generator beforehand
res <- rma(yi, vi, mods = ~ ablat, data=dat, method="DL")
set.seed(1234)
permutest(res)## Running 1000 iterations for an approximate permutation test.
##
## Test of Moderators (coefficient 2):¹
## QM(df = 1) = 18.8452, p-val = 0.0070
##
## Model Results:
##
## estimate se zval pval¹ ci.lb ci.ub
## intrcpt 0.2595 0.2323 1.1172 0.6300 -0.1958 0.7149
## ablat -0.0292 0.0067 -4.3411 0.0070 -0.0424 -0.0160 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## 1) p-values based on permutation testing
23) Subgroups within Studies
# Table 23.1: Independent subgroups – five fictional studies
dat <- read.table(header=TRUE, text = "
study subgroup es var
'Study 1' 'Stage 1' 0.3 0.05
'Study 1' 'Stage 2' 0.1 0.05
'Study 2' 'Stage 1' 0.2 0.02
'Study 2' 'Stage 2' 0.1 0.02
'Study 3' 'Stage 1' 0.4 0.05
'Study 3' 'Stage 2' 0.2 0.05
'Study 4' 'Stage 1' 0.2 0.01
'Study 4' 'Stage 2' 0.1 0.01
'Study 5' 'Stage 1' 0.4 0.06
'Study 5' 'Stage 2' 0.3 0.06
")
dat## study subgroup es var
## 1 Study 1 Stage 1 0.3 0.05
## 2 Study 1 Stage 2 0.1 0.05
## 3 Study 2 Stage 1 0.2 0.02
## 4 Study 2 Stage 2 0.1 0.02
## 5 Study 3 Stage 1 0.4 0.05
## 6 Study 3 Stage 2 0.2 0.05
## 7 Study 4 Stage 1 0.2 0.01
## 8 Study 4 Stage 2 0.1 0.01
## 9 Study 5 Stage 1 0.4 0.06
## 10 Study 5 Stage 2 0.3 0.06
## Equal-Effects Model (k = 10)
##
## I^2 (total heterogeneity / total variability): 0.00%
## H^2 (total variability / sampling variability): 0.38
##
## Test for Heterogeneity:
## Q(df = 9) = 3.4462, p-val = 0.9440
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.1855 0.0492 3.7710 0.0002 0.0891 0.2819 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# compute combined effect across subgroups within studies
dat <- escalc(yi=es, vi=var, data=dat, var.names=c("es","var"))
agg <- aggregate(dat, cluster=study, struct="ID", select=c("study","es","var"))
agg## study es var
## 1 Study 1 0.2000 0.0250
## 2 Study 2 0.1500 0.0100
## 3 Study 3 0.3000 0.0250
## 4 Study 4 0.1500 0.0050
## 5 Study 5 0.3500 0.0300
## Equal-Effects Model (k = 5)
##
## I^2 (total heterogeneity / total variability): 0.00%
## H^2 (total variability / sampling variability): 0.45
##
## Test for Heterogeneity:
## Q(df = 4) = 1.8129, p-val = 0.7701
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.1855 0.0492 3.7710 0.0002 0.0891 0.2819 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
24) Multiple Outcomes / Time-Points
# Table 24.1: Multiple outcomes – five fictional studies
dat <- read.table(header=TRUE, text = "
study outcome es var cor
'Study 1' Reading 0.3 0.05 0.5
'Study 1' Math 0.1 0.05 0.5
'Study 2' Reading 0.2 0.02 0.6
'Study 2' Math 0.1 0.02 0.6
'Study 3' Reading 0.4 0.05 0.6
'Study 3' Math 0.2 0.05 0.6
'Study 4' Reading 0.2 0.01 0.4
'Study 4' Math 0.1 0.01 0.4
'Study 5' Reading 0.4 0.06 0.8
'Study 5' Math 0.3 0.06 0.8
")
dat## study outcome es var cor
## 1 Study 1 Reading 0.3 0.05 0.5
## 2 Study 1 Math 0.1 0.05 0.5
## 3 Study 2 Reading 0.2 0.02 0.6
## 4 Study 2 Math 0.1 0.02 0.6
## 5 Study 3 Reading 0.4 0.05 0.6
## 6 Study 3 Math 0.2 0.05 0.6
## 7 Study 4 Reading 0.2 0.01 0.4
## 8 Study 4 Math 0.1 0.01 0.4
## 9 Study 5 Reading 0.4 0.06 0.8
## 10 Study 5 Math 0.3 0.06 0.8
# compute combined effect across outcomes and corresponding sampling variance
dat <- escalc(yi=es, vi=var, data=dat, var.names=c("es","var"))
agg <- aggregate(dat, cluster=study, struct="CS", select=c("study","es","var"),
rho=tapply(dat$cor, dat$study, head, 1))
agg## study es var
## 1 Study 1 0.2000 0.0375
## 2 Study 2 0.1500 0.0160
## 3 Study 3 0.3000 0.0400
## 4 Study 4 0.1500 0.0070
## 5 Study 5 0.3500 0.0540
## Equal-Effects Model (k = 5)
##
## I^2 (total heterogeneity / total variability): 0.00%
## H^2 (total variability / sampling variability): 0.27
##
## Test for Heterogeneity:
## Q(df = 4) = 1.0897, p-val = 0.8959
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.1819 0.0602 3.0193 0.0025 0.0638 0.3000 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# compute difference between outcomes and corresponding sampling variances
dif <- do.call(rbind, lapply(split(dat, dat$study), function(x) {
c(es = x$es[1] - x$es[2], var = x$var[1] + x$var[2] - 2*x$cor[1]*sqrt(x$var[1]*x$var[2]))}))
dif## es var
## Study 1 0.2 0.050
## Study 2 0.1 0.016
## Study 3 0.2 0.040
## Study 4 0.1 0.012
## Study 5 0.1 0.024
## Equal-Effects Model (k = 5)
##
## I^2 (total heterogeneity / total variability): 0.00%
## H^2 (total variability / sampling variability): 0.09
##
## Test for Heterogeneity:
## Q(df = 4) = 0.3629, p-val = 0.9854
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.1194 0.0656 1.8199 0.0688 -0.0092 0.2479 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# could also approach this as a multivariate meta-analysis (not shown in the book)
V <- lapply(split(dat, dat$study), function(x) {
Vi <- matrix(x$var, nrow=2, ncol=2)
Vi[1,2] <- Vi[2,1] <- x$cor[1]*sqrt(Vi[1,1]*Vi[2,2])
Vi})
rma.mv(es, V, mods = ~ outcome, data=dat)## Multivariate Meta-Analysis Model (k = 10; method: REML)
##
## Variance Components: none
##
## Test for Residual Heterogeneity:
## QE(df = 8) = 1.4526, p-val = 0.9935
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 3.3121, p-val = 0.0688
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.1222 0.0686 1.7818 0.0748 -0.0122 0.2566 .
## outcomeReading 0.1194 0.0656 1.8199 0.0688 -0.0092 0.2479 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
30) Publication Bias
## study author year country design cases or or.lb or.ub yi vi
## 1 1 Garfinkel 1981 USA cohort 153 1.18 0.90 1.54 0.1655 0.0188
## 2 2 Hirayama 1984 Japan cohort 200 1.45 1.02 2.08 0.3716 0.0330
## 3 3 Butler 1988 USA cohort 8 2.02 0.48 8.56 0.7031 0.5402
## 4 4 Cardenas 1997 USA cohort 150 1.20 0.80 1.60 0.1823 0.0313
## 5 5 Chan 1982 Hong Kong case-control 84 0.75 0.43 1.30 -0.2877 0.0797
## 6 6 Correa 1983 USA case-control 22 2.07 0.81 5.25 0.7275 0.2273
## 7 7 Trichopolous 1983 Greece case-control 62 2.13 1.19 3.83 0.7561 0.0889
## 8 8 Buffler 1984 USA case-control 41 0.80 0.34 1.90 -0.2231 0.1927
## 9 9 Kabat 1984 USA case-control 24 0.79 0.25 2.45 -0.2357 0.3390
## 10 10 Lam 1985 Hong Kong case-control 60 2.01 1.09 3.72 0.6981 0.0981
## 11 11 Garfinkel 1985 USA case-control 134 1.23 0.81 1.87 0.2070 0.0456
## 12 12 Wu 1985 USA case-control 29 1.20 0.50 3.30 0.1823 0.2317
## 13 13 Akiba 1986 Japan case-control 94 1.52 0.87 2.63 0.4187 0.0796
## 14 14 Lee 1986 UK case-control 32 1.03 0.41 2.55 0.0296 0.2174
## 15 15 Koo 1987 Hong Kong case-control 86 1.55 0.90 2.67 0.4383 0.0770
## 16 16 Pershagen 1987 Sweden case-control 70 1.03 0.61 1.74 0.0296 0.0715
## 17 17 Humble 1987 USA case-control 20 2.34 0.81 6.75 0.8502 0.2926
## 18 18 Lam 1987 Hong Kong case-control 199 1.65 1.16 2.35 0.5008 0.0324
## 19 19 Gao 1987 China case-control 246 1.19 0.82 1.73 0.1740 0.0363
## 20 20 Brownson 1987 USA case-control 19 1.52 0.39 5.96 0.4187 0.4839
## 21 21 Geng 1988 China case-control 54 2.16 1.08 4.29 0.7701 0.1238
## 22 22 Shimizu 1988 Japan case-control 90 1.08 0.64 1.82 0.0770 0.0711
## 23 23 Inoue 1988 Japan case-control 22 2.55 0.74 8.78 0.9361 0.3982
## 24 24 Kalandidi 1990 Greece case-control 90 1.62 0.90 2.91 0.4824 0.0896
## 25 25 Sobue 1990 Japan case-control 144 1.06 0.74 1.52 0.0583 0.0337
## 26 26 Wu-Williams 1990 China case-control 417 0.79 0.62 1.02 -0.2357 0.0161
## 27 27 Liu 1991 China case-control 54 0.74 0.32 1.69 -0.3011 0.1802
## 28 28 Jockel 1991 Germany case-control 23 2.27 0.75 6.82 0.8198 0.3171
## 29 29 Brownson 1992 USA case-control 431 0.97 0.78 1.21 -0.0305 0.0125
## 30 30 Stockwell 1992 USA case-control 210 1.60 0.80 3.00 0.4700 0.1137
## 31 31 Du 1993 China case-control 75 1.19 0.66 2.13 0.1740 0.0893
## 32 32 Liu 1993 China case-control 38 1.66 0.73 3.78 0.5068 0.1760
## 33 33 Fontham 1994 USA case-control 651 1.26 1.04 1.54 0.2311 0.0100
## 34 34 Kabat 1995 USA case-control 67 1.10 0.62 1.96 0.0953 0.0862
## 35 35 Zaridze 1995 Russia case-control 162 1.66 1.12 2.45 0.5068 0.0399
## 36 36 Sun 1996 China case-control 230 1.16 0.80 1.69 0.1484 0.0364
## 37 37 Wang 1996 China case-control 135 1.11 0.67 1.84 0.1044 0.0664
# note: the 'yi' values are log odds ratios, not log risk ratios as suggested in
# the book; also, there are slight differences between the data used in the book
# and those in 'dat.hackshaw1998'; the differences are negligible though# equal-effects model
res <- rma(yi, vi, data=dat, method="EE", slab=paste0(author, ", ", year))
res## Equal-Effects Model (k = 37)
##
## I^2 (total heterogeneity / total variability): 24.21%
## H^2 (total variability / sampling variability): 1.32
##
## Test for Heterogeneity:
## Q(df = 36) = 47.4979, p-val = 0.0952
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.1858 0.0373 4.9797 <.0001 0.1126 0.2589 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## pred ci.lb ci.ub
## 1.204 1.119 1.295
dat$weights <- formatC(weights(res), format="f", digits=2)
dat$pvals <- formatC(summary(dat)$pval, format="f", digits=3)# Figure 30.1
par(mar=c(4,4,2,2))
forest(res, xlim=c(-5,8), atransf=exp, at=log(c(0.1, 1, 10)),
header=TRUE, top=2, mlab="Summary",
ilab=data.frame(dat$weights, dat$pvals), ilab.xpos=c(3.5,5), ilab.pos=2,
efac=c(0,1,1.5), order="prec", cex=0.7)
text(3.5, 39, "Weight", pos=2, font=2, cex=0.7)
text(5.0, 39, "P-Value", pos=2, font=2, cex=0.7)
# Figure 30.2
funnel(res, xlim=c(-2,2), ylim=c(0,0.8))
title("Funnel Plot of Standard Error by Log Risk Ratio")
## Fail-safe N Calculation Using the Rosenthal Approach
##
## Observed Significance Level: <.0001
## Target Significance Level: 0.05
##
## Fail-safe N: 393
## Fail-safe N Calculation Using the Orwin Approach
##
## Average Effect Size: 0.1858
## Target Effect Size: 0.0488
##
## Fail-safe N: 104
## Estimated number of missing studies on the left side: 7 (SE = 4.0405)
##
## Equal-Effects Model (k = 44)
##
## I^2 (total heterogeneity / total variability): 30.44%
## H^2 (total variability / sampling variability): 1.44
##
## Test for Heterogeneity:
## Q(df = 43) = 61.8207, p-val = 0.0313
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.1556 0.0364 4.2706 <.0001 0.0842 0.2270 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## pred ci.lb ci.ub
## 1.168 1.088 1.255
# Figure 30.3
funnel(taf, xlim=c(-2,2), ylim=c(0,0.8), pch=21, pch.fill=19)
title("Funnel Plot of Standard Error by Log Risk Ratio")
## estimate se zval pvals ci.lb ci.ub Q Qp I2 H2
## Fontham, 1994 0.2311 0.1001 2.3078 0.0210 0.0348 0.4274 0.0000 1.0000 0.0000 1.0000
## Brownson, 1992 0.1149 0.0747 1.5392 0.1238 -0.0314 0.2612 3.0306 0.0817 67.0036 3.0306
## Wu-Williams, 1990 0.0249 0.0644 0.3863 0.6993 -0.1013 0.1510 8.6954 0.0129 76.9994 4.3477
## Garfinkel, 1981 0.0503 0.0583 0.8632 0.3880 -0.0639 0.1645 9.5586 0.0227 68.6146 3.1862
## Cardenas, 1997 0.0632 0.0553 1.1425 0.2533 -0.0452 0.1717 10.0616 0.0394 60.2448 2.5154
## Lam, 1987 0.1009 0.0529 1.9086 0.0563 -0.0027 0.2046 15.4550 0.0086 67.6481 3.0910
## Hirayama, 1984 0.1221 0.0508 2.4037 0.0162 0.0225 0.2216 17.4983 0.0076 65.7109 2.9164
## Sobue, 1990 0.1175 0.0489 2.4013 0.0163 0.0216 0.2135 17.6104 0.0139 60.2508 2.5158
## Gao, 1987 0.1210 0.0474 2.5530 0.0107 0.0281 0.2139 17.6927 0.0237 54.7837 2.2116
## Sun, 1996 0.1226 0.0460 2.6653 0.0077 0.0325 0.2128 17.7121 0.0387 49.1874 1.9680
## Zaridze, 1995 0.1420 0.0448 3.1671 0.0015 0.0541 0.2299 21.2273 0.0196 52.8908 2.1227
## Garfinkel, 1985 0.1447 0.0439 3.2988 0.0010 0.0587 0.2307 21.3162 0.0302 48.3959 1.9378
## Wang, 1996 0.1436 0.0433 3.3200 0.0009 0.0588 0.2284 21.3400 0.0456 43.7676 1.7783
## Shimizu, 1988 0.1419 0.0427 3.3234 0.0009 0.0582 0.2256 21.4009 0.0654 39.2548 1.6462
## Pershagen, 1987 0.1391 0.0422 3.2992 0.0010 0.0565 0.2217 21.5730 0.0878 35.1040 1.5409
## Koo, 1987 0.1459 0.0417 3.4991 0.0005 0.0642 0.2275 22.7096 0.0905 33.9487 1.5140
## Akiba, 1986 0.1517 0.0412 3.6784 0.0002 0.0709 0.2325 23.6245 0.0980 32.2736 1.4765
## Chan, 1982 0.1425 0.0408 3.4924 0.0005 0.0625 0.2225 25.9972 0.0745 34.6084 1.5292
## Kabat, 1995 0.1416 0.0404 3.5038 0.0005 0.0624 0.2208 26.0226 0.0992 30.8292 1.4457
## Trichopolous, 1983 0.1527 0.0400 3.8126 0.0001 0.0742 0.2312 30.1928 0.0494 37.0711 1.5891
## Du, 1993 0.1531 0.0397 3.8561 0.0001 0.0753 0.2309 30.1978 0.0667 33.7699 1.5099
## Kalandidi, 1990 0.1587 0.0393 4.0345 0.0001 0.0816 0.2359 31.3873 0.0674 33.0939 1.4946
## Lam, 1985 0.1671 0.0390 4.2809 0.0000 0.0906 0.2437 34.3079 0.0457 35.8748 1.5595
## Stockwell, 1992 0.1711 0.0388 4.4129 0.0000 0.0951 0.2472 35.1040 0.0508 34.4805 1.5263
## Geng, 1988 0.1783 0.0385 4.6261 0.0000 0.1028 0.2539 37.9668 0.0349 36.7869 1.5820
## Liu, 1993 0.1811 0.0384 4.7172 0.0000 0.1058 0.2563 38.5748 0.0406 35.1909 1.5430
## Liu, 1991 0.1772 0.0382 4.6342 0.0000 0.1022 0.2521 39.8544 0.0403 34.7625 1.5329
## Buffler, 1984 0.1742 0.0381 4.5726 0.0000 0.0995 0.2488 40.6798 0.0442 33.6280 1.5067
## Lee, 1986 0.1732 0.0380 4.5626 0.0000 0.0988 0.2476 40.7753 0.0563 31.3310 1.4563
## Correa, 1983 0.1767 0.0378 4.6693 0.0000 0.1025 0.2509 42.1187 0.0548 31.1469 1.4524
## Wu, 1985 0.1767 0.0377 4.6846 0.0000 0.1028 0.2507 42.1188 0.0699 28.7729 1.4040
## Humble, 1987 0.1800 0.0376 4.7826 0.0000 0.1062 0.2537 43.6614 0.0652 28.9991 1.4084
## Jockel, 1991 0.1828 0.0375 4.8690 0.0000 0.1092 0.2564 44.9464 0.0641 28.8040 1.4046
## Kabat, 1984 0.1811 0.0375 4.8329 0.0000 0.1077 0.2545 45.4610 0.0729 27.4102 1.3776
## Inoue, 1988 0.1837 0.0374 4.9123 0.0000 0.1104 0.2571 46.8874 0.0696 27.4859 1.3790
## Brownson, 1987 0.1844 0.0374 4.9375 0.0000 0.1112 0.2576 47.0012 0.0846 25.5338 1.3429
## Butler, 1988 0.1858 0.0373 4.9797 0.0000 0.1126 0.2589 47.4979 0.0952 24.2072 1.3194
# Figure 30.4
par(mar=c(4,4,2,2))
forest(rcm, xlim=c(-1.5,1.5), atransf=exp, at=log(c(0.5, 1, 2)),
header=TRUE, top=2, efac=c(0,1,1.5), cex=0.7)
33) Simpson’s Paradox
## study category non.pos non.neg cir.pos cir.neg
## 1 Bwayo 1 high-risk group 92 86 160 612
## 2 Bwayo 2 high-risk group 22 46 37 200
## 3 Kreiss high-risk group 59 18 254 168
## 4 Hira high-risk group 418 172 10 10
## 5 Cameron high-risk group 18 61 6 208
## 6 Pepin high-risk group 5 42 13 243
## 7 Greenblatt high-risk group 11 28 8 68
## 8 Diallo high-risk group 38 46 212 873
## 9 Simonsen high-risk group 17 70 21 232
## 10 Tyndall high-risk group 85 93 105 527
## 11 Nasio high-risk group 86 78 137 580
## 12 Mehendal high-risk group 837 3411 38 253
## 13 Bollinger high-risk group 50 241 1 14
## 14 Chiasson high-risk group 36 797 14 542
## 15 Sassan high-risk group 75 18 415 221
## 16 Hunter partner study 43 330 165 3765
## 17 Carael partner study 90 105 34 45
## 18 Chao partner study 442 4844 75 232
## 19 Moss partner study 24 16 12 17
## 20 Allen partner study 275 612 132 324
## 21 Sedlin partner study 32 26 33 40
## 22 Konde-Luc partner study 153 1516 6 127
## 23 Barongo population survey 55 1356 42 642
## 24 Grosskurth population survey 158 4604 61 1026
## 25 Van de Perre population survey 46 224 6 26
## 26 Seed population survey 171 422 51 192
## 27 Malamba population survey 111 114 21 47
## 28 Quigley population survey 101 272 48 121
## 29 Urassa 1 population survey 56 1301 42 600
## 30 Urassa 2 population survey 105 2040 32 426
## 31 Urassa 3 population survey 38 309 19 158
## 32 Urassa 4 population survey 112 716 54 692
## 33 Urassa 5 population survey 101 365 48 136
# naive pooling by summing up the counts and then computing the odds ratio and
# corresponding confidence interval based on the aggregated table
agg <- escalc(measure="OR", ai=sum(dat$non.pos), bi=sum(dat$non.neg),
ci=sum(dat$cir.pos), di=sum(dat$cir.neg))
summary(agg, transf=exp, digits=2)## yi ci.lb ci.ub
## 1 0.94 0.89 0.99
# calculate log odds ratios and corresponding sampling variances
dat <- escalc(measure="OR", ai=non.pos, bi=non.neg, ci=cir.pos, di=cir.neg, data=dat, slab=study)## Random-Effects Model (k = 33; tau^2 estimator: DL)
##
## tau^2 (estimated amount of total heterogeneity): 0.6397 (SE = 0.2103)
## tau (square root of estimated tau^2 value): 0.7998
## I^2 (total heterogeneity / total variability): 92.36%
## H^2 (total variability / sampling variability): 13.09
##
## Test for Heterogeneity:
## Q(df = 32) = 419.0087, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.5132 0.1499 3.4248 0.0006 0.2195 0.8069 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## pred ci.lb ci.ub pi.lb pi.ub
## 1.67 1.25 2.24 0.34 8.23
# Figure 33.1
par(mar=c(4,4,2,2))
forest(res, xlim=c(-5,7), atransf=exp, at=log(c(0.1, 1, 10, 100)),
header=TRUE, top=2, mlab="Summary", efac=c(0,1), cex=0.7,
fonts=c("sans", "inconsolata"))
# mixed-effects meta-regression model with category as predictor
sub <- rma(yi, vi, mods = ~ category - 1, data=dat, method="DL")
sub## Mixed-Effects Model (k = 33; tau^2 estimator: DL)
##
## tau^2 (estimated amount of residual heterogeneity): 0.3445 (SE = 0.1217)
## tau (square root of estimated tau^2 value): 0.5869
## I^2 (residual heterogeneity / unaccounted variability): 86.15%
## H^2 (unaccounted variability / sampling variability): 7.22
##
## Test for Residual Heterogeneity:
## QE(df = 30) = 216.5475, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## QM(df = 3) = 39.1767, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## categoryhigh-risk group 1.0973 0.1772 6.1928 <.0001 0.7500 1.4446 ***
## categorypartner study 0.2188 0.2477 0.8833 0.3771 -0.2668 0.7044
## categorypopulation survey -0.0410 0.1914 -0.2142 0.8304 -0.4161 0.3341
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# note: by removing the intercept, the three coefficients directly provide the
# estimated average log odds ratios for the three categories# estimated average odds ratios and corresponding 95% CIs
predict(sub, newmods=diag(3), transf=exp, digits=3)## pred ci.lb ci.ub pi.lb pi.ub
## 1 2.996 2.117 4.240 0.901 9.963
## 2 1.245 0.766 2.023 0.357 4.338
## 3 0.960 0.660 1.397 0.286 3.219
# Figure 33.2
par(mar=c(4,4,2,2))
est <- coef(sub)
var <- diag(vcov(sub))
cat <- c("High risk", "Partner", "Random")
forest(est, var, slab=cat, header=c("Studies", "OR [95% CI]"), top=2,
xlim=c(-2,4), ylim=c(-1.5,5), xlab="Odds Ratio (log scale)",
atransf=exp, at=log(c(0.5, 1, 2, 4, 8)), digits=3L, efac=0)
abline(h=0)
res <- rma(est, var, method="DL")
addpoly(res, row=-1, efac=3, mlab="Summary")
36) Methods Based on p-Values
# sign test
dat <- dat.lau1992
dat <- escalc(measure="RR", ai=ai, n1i=n1i, ci=ci, n2i=n2i, data=dat, slab=trial)
table(sign(dat$yi))##
## -1 1
## 25 8
##
## Exact binomial test
##
## data: table(sign(dat$yi))
## number of successes = 25, number of trials = 33, p-value = 0.004551
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
## 0.5774107 0.8890767
## sample estimates:
## probability of success
## 0.7575758
## trial year ai n1i ci n2i yi vi pval
## 1 Fletcher 1959 1 12 4 11 -1.4733 1.0758 7.773373e-02
## 2 Dewar 1963 4 21 7 21 -0.5596 0.2976 1.524947e-01
## 3 European 1 1969 20 83 15 84 0.2997 0.0927 8.374752e-01
## 4 European 2 1971 69 373 94 357 -0.3530 0.0196 5.892704e-03
## 5 Heikinheimo 1971 22 219 17 207 0.2015 0.0949 7.434715e-01
## 6 Italian 1971 19 164 18 157 0.0104 0.0957 5.134679e-01
## 7 Australian 1 1973 26 264 32 253 -0.2502 0.0620 1.574346e-01
## 8 Frankfurt 2 1973 13 102 29 104 -0.7829 0.0920 4.919634e-03
## 9 NHLBI SMIT 1974 7 53 3 54 0.8660 0.4388 9.044458e-01
## 10 Frank 1975 6 55 6 53 -0.0370 0.2963 4.728727e-01
## 11 Valere 1975 11 49 9 42 0.0465 0.1578 5.466124e-01
## 12 Klein 1976 4 14 1 9 0.9445 1.0675 8.196760e-01
## 13 UK Collab 1976 38 302 40 293 -0.0815 0.0446 3.496831e-01
## 14 Austrian 1977 37 352 65 376 -0.4975 0.0369 4.805083e-03
## 15 Australian 2 1977 25 123 31 107 -0.3545 0.0548 6.495541e-02
## 16 Lasierra 1977 1 13 3 11 -1.2657 1.1655 1.205253e-01
## 17 N Ger Collab 1977 63 249 51 234 0.1492 0.0272 8.171786e-01
## 18 Witchitz 1977 5 32 5 26 -0.2076 0.3303 3.589391e-01
## 19 European 3 1979 18 156 30 159 -0.4918 0.0762 3.740345e-02
## 20 ISAM 1986 54 859 63 882 -0.1277 0.0321 2.379301e-01
## 21 GISSI-1 1986 628 5860 758 5852 -0.1895 0.0026 9.268991e-05
## 22 Olson 1986 1 28 2 24 -0.8473 1.4226 2.387337e-01
## 23 Baroffio 1986 0 29 6 30 -2.5322 2.0883 3.986429e-02
## 24 Schreiber 1986 1 19 3 19 -1.0986 1.2281 1.607541e-01
## 25 Cribier 1986 1 21 1 23 0.0910 1.9089 5.262489e-01
## 26 Sainsous 1986 3 49 6 49 -0.6931 0.4592 1.531781e-01
## 27 Durand 1987 3 35 4 29 -0.4757 0.5203 2.547720e-01
## 28 White 1987 2 107 12 112 -1.7461 0.5651 1.009383e-02
## 29 Bassand 1987 4 52 7 55 -0.5035 0.3554 1.991752e-01
## 30 Vlay 1988 1 13 2 12 -0.7732 1.3397 2.520674e-01
## 31 Kennedy 1988 12 191 17 177 -0.4244 0.1313 1.207104e-01
## 32 ISIS-2 1988 791 8592 1029 8595 -0.2627 0.0020 2.189447e-09
## 33 Wisenberg 1988 2 41 5 25 -1.4110 0.6356 3.837903e-02
Although the computations underlying Fisher’s and Stouffer’s methods are not difficult, we can make use of the poolr package to automate things.
# install the 'poolr' package (if it is not already installed) and load it
if (!require(poolr)) {
install.packages("poolr")
library(poolr)
}## combined p-values with: Fisher's method
## number of p-values combined: 33
## test statistic: 178.536 ~ chi-square(df = 66)
## adjustment: none
## combined p-value: 2.644007e-12
## combined p-values with: Stouffer's method
## number of p-values combined: 33
## test statistic: 5.862 ~ N(0,1)
## adjustment: none
## combined p-value: 2.280421e-09
37) Methods for Dichotomous Data
# Table 37.2 (like Table 14.4)
dat <- read.table(header=TRUE, text = "
study events1 nonevents1 n1 events2 nonevents2 n2
Saint 12 53 65 16 49 65
Kelly 8 32 40 10 30 40
Pilbeam 14 66 80 19 61 80
Lane 25 375 400 80 320 400
Wright 8 32 40 11 29 40
Day 16 49 65 18 47 65")
dat## study events1 nonevents1 n1 events2 nonevents2 n2
## 1 Saint 12 53 65 16 49 65
## 2 Kelly 8 32 40 10 30 40
## 3 Pilbeam 14 66 80 19 61 80
## 4 Lane 25 375 400 80 320 400
## 5 Wright 8 32 40 11 29 40
## 6 Day 16 49 65 18 47 65
## Equal-Effects Model (k = 6)
##
## I^2 (total heterogeneity / total variability): 52.78%
## H^2 (total variability / sampling variability): 2.12
##
## Test for Heterogeneity:
## Q(df = 5) = 10.5887, p-val = 0.0602
##
## Model Results (log scale):
##
## estimate se zval pval ci.lb ci.ub
## -0.7539 0.1502 -5.0211 <.0001 -1.0482 -0.4596
##
## Model Results (OR scale):
##
## estimate ci.lb ci.ub
## 0.4705 0.3506 0.6315
##
## Cochran-Mantel-Haenszel Test: CMH = 25.2668, df = 1, p-val < 0.0001
## Tarone's Test for Heterogeneity: X^2 = 10.7482, df = 5, p-val = 0.0566
## Equal-Effects Model (k = 6)
##
## I^2 (total heterogeneity / total variability): 49.25%
## H^2 (total variability / sampling variability): 1.97
##
## Test for Heterogeneity:
## Q(df = 5) = 9.8526, p-val = 0.0795
##
## Model Results (log scale):
##
## estimate se zval pval ci.lb ci.ub
## -0.7322 0.1436 -5.0984 <.0001 -1.0137 -0.4507
##
## Model Results (OR scale):
##
## estimate ci.lb ci.ub
## 0.4808 0.3629 0.6372
38) Psychometric Meta-Analysis
# Table 38.1: Fictional data for psychometric meta-analysis
dat <- read.table(header=TRUE, text = "
study n r rel
'University 1' 130 0.24 0.75
'University 2' 90 0.11 0.75
'Private 1' 30 0.05 0.60
'Private 2' 25 0.17 0.60
'Volunteer 1' 50 0.38 0.90
'Volunteer 2' 65 0.50 0.90")# psychometric meta-analysis with the observed (attenuated) correlations
dat <- escalc(measure="COR", ri=r, ni=n, data=dat, vtype="AV")
res <- rma(yi, vi, weights=n, data=dat, method="HS")
res## Random-Effects Model (k = 6; tau^2 estimator: HS)
##
## tau^2 (estimated amount of total heterogeneity): 0.0069 (SE = 0.0107)
## tau (square root of estimated tau^2 value): 0.0828
## I^2 (total heterogeneity / total variability): 31.92%
## H^2 (total variability / sampling variability): 1.47
##
## Test for Heterogeneity:
## Q(df = 5) = 9.0059, p-val = 0.1088
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.2522 0.0615 4.0994 <.0001 0.1316 0.3727 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# psychometric meta-analysis with the corrected (unattenuated) correlations
dat$yi.c <- dat$yi / sqrt(dat$rel)
dat$vi.c <- dat$vi / (dat$rel)
res.u <- rma(yi.c, vi.c, weights=1/vi.c, data=dat, method="HS")
res.u## Random-Effects Model (k = 6; tau^2 estimator: HS)
##
## tau^2 (estimated amount of total heterogeneity): 0.0042 (SE = 0.0115)
## tau (square root of estimated tau^2 value): 0.0649
## I^2 (total heterogeneity / total variability): 18.18%
## H^2 (total variability / sampling variability): 1.22
##
## Test for Heterogeneity:
## Q(df = 5) = 7.4283, p-val = 0.1907
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.2949 0.0624 4.7296 <.0001 0.1727 0.4172 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# explained variance
S2 <- sum(dat$n * (dat$r - coef(res))^2) / sum(dat$n)
round((S2 - res.u$tau2) / S2, digits=4)## [1] 0.7955
42) Cumulative Meta-Analysis
# analysis of the Lau et al. (1992) data using (log) risk ratios
dat <- dat.lau1992
dat <- escalc(measure="RR", ai=ai, n1i=n1i, ci=ci, n2i=n2i, data=dat, slab=trial)
res <- rma(yi, vi, data=dat, method="DL")
res## Random-Effects Model (k = 33; tau^2 estimator: DL)
##
## tau^2 (estimated amount of total heterogeneity): 0.0077 (SE = 0.0127)
## tau (square root of estimated tau^2 value): 0.0876
## I^2 (total heterogeneity / total variability): 16.87%
## H^2 (total variability / sampling variability): 1.20
##
## Test for Heterogeneity:
## Q(df = 32) = 38.4942, p-val = 0.1991
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.2312 0.0468 -4.9345 <.0001 -0.3230 -0.1394 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## pred ci.lb ci.ub pi.lb pi.ub
## 0.794 0.724 0.870 0.653 0.964
# Figure 42.1
par(mar=c(4,4,2,2))
forest(res, xlim=c(-10,9), atransf=exp, at=log(c(0.01, 0.1, 1, 10, 100)),
header=TRUE, top=2, ilab=dat$year, ilab.xpos=-6, digits=3L, cex=0.8,
efac=c(0,1), fonts=c("sans", "inconsolata"))
text(-6, 35, "Year", font=2, cex=0.8)
# Figure 42.2
par(mar=c(4,4,2,2))
forest(rcm, xlim=c(-5,3), atransf=exp, at=log(c(0.1, 0.25, 0.5, 1, 2)),
header=TRUE, top=2, ilab=year, ilab.xpos=-3, digits=3L, cex=0.8,
efac=c(0,1))
text(-3, 35, "Year", font=2, cex=0.8)
License
This documented is licensed under the following license: CC Attribution-Noncommercial-Share Alike 4.0 International.