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:

install.packages("metafor")

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

library(metafor)

A few additional notes:

  1. Results are only reproduced for chapters containing worked examples.
  2. 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).
  3. 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 metafor package 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.
  4. 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
res <- rma(yi, vi, data=dat, method="DL")
res
## 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
predict(res, transf=exp, digits=2) # summary effect size (risk ratio)
##  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
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
predict(res, transf=exp, digits=2) # summary effect size (risk ratio)
##  pred ci.lb ci.ub pi.lb pi.ub 
##  0.79  0.72  0.87  0.65  0.96
formatC(res$pval, format="f", digits=7) # p-value of the summary estimate
## [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
# note: since n1i=n2i in this example, the results are exactly the same
# 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
dat <- escalc("MC", m1i=105, m2i=100, sd1i=10, sd2i=10, ni=50, ri=0.5)
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
# 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

# risk ratio (log transformed)

dat <- escalc("RR", ai=5, n1i=100, ci=10, n2i=100)
summary(dat)
##        yi     vi    sei      zi   pval   ci.lb  ci.ub 
## 1 -0.6931 0.2800 0.5292 -1.3099 0.1902 -1.7303 0.3440
# odds ratio (log transformed)

dat <- escalc("OR", ai=5, n1i=100, ci=10, n2i=100)
summary(dat)
##        yi     vi    sei      zi   pval   ci.lb  ci.ub 
## 1 -0.7472 0.3216 0.5671 -1.3175 0.1877 -1.8588 0.3643
# risk difference

dat <- escalc("RD", ai=5, n1i=100, ci=10, n2i=100)
summary(dat)
##        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

# r-to-z transformed correlation coefficient

dat <- escalc("ZCOR", ri=0.50, ni=100)
summary(dat)
##       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 analysis

res <- rma(yi, vi, data=dat, method="EE", digits=2)
res
## 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 analysis

res <- rma(yi, vi, data=dat, method="DL", digits=2)
res
## 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 analysis

res <- rma(yi, vi, data=dat, method="EE", digits=2)
res
## 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
predict(res, transf=exp)
##  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 analysis

res <- rma(yi, vi, data=dat, method="DL", digits=2)
res
## 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 analysis

res <- rma(yi, vi, data=dat, method="EE", digits=2)
res
## 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
predict(res, transf=transf.ztor)
##  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 analysis

res <- rma(yi, vi, data=dat, method="DL", digits=2)
res
## 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
predict(res, transf=transf.ztor)
##  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

  1. The prediction interval is computed in a slightly different way in the book compared to how this is done (by default) in the metafor package. For more details, see here. However, by using argument pi.type="Riley", we can obtain the same results as those provided in the book.
  2. The confidence intervals for \(\tau^2\) and \(I^2\) are also computed in a different way. By default, the metafor package 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 argument type="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 analysis

res <- rma(yi, vi, data=dat, method="DL")
res
## 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
predict(res, pi.type="Riley")
##    pred     se  ci.lb  ci.ub   pi.lb  pi.ub 
##  0.3582 0.1052 0.1520 0.5645 -0.2525 0.9690
confint(res, type="HT")
##        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 analysis

res <- rma(yi, vi, data=dat, method="DL")
res
## 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
predict(res, transf=exp, pi.type="Riley")
##    pred  ci.lb  ci.ub  pi.lb  pi.ub 
##  0.5676 0.3554 0.9065 0.1499 2.1492
confint(res, type="HT")
##        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 analysis

res <- rma(yi, vi, data=dat, method="DL")
res
## 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
predict(res, transf=transf.ztor, pi.type="Riley")
##    pred  ci.lb  ci.ub   pi.lb  pi.ub 
##  0.4875 0.2714 0.6568 -0.3271 0.8865
confint(res, type="HT")
##        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 for all 10 studies

res <- rma(g, var, data=dat, method="EE", slab=study)
res
## 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
# compare effects of A and B

dif <- rma(g, var, mods = ~ type, method="FE", data=dat)
dif
## 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
# note: this allows for a different tau^2 value within the two subgroups
# 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 anyway

RE 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
# note: fitted model without intercept so the coefficients are the estimated effects for A and B
# 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 B

20) 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

res <- rma(yi, vi, data=dat, method="EE")
res
## 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

res <- rma(yi, vi, data=dat, method="DL")
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     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
predict(res, transf=exp, digits=2)
##  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.

dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)
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
# 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
# treat each subgroup as a separate study

res <- rma(es, var, data=dat, method="EE")
res
## 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
# use study as the unit of analysis

res <- rma(es, var, data=agg, method="EE")
res
## 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
# note: correlations from Table 24.3
# 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 using the aggregated data

rma(es, var, data=agg, method="EE")
## 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 using the differences data

rma(es, var, data=dif, method="EE")
## 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
# note: the coefficient for 'outcomeReading' is the estimate of the difference

30) Publication Bias

# copy data into 'dat' and examine data

dat <- dat.hackshaw1998
dat
##    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
predict(res, transf=exp, digits=3)
##   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")

# Rosenthal's Fail-safe N

fsn(yi, vi, data=dat)
## Fail-safe N Calculation Using the Rosenthal Approach
## 
## Observed Significance Level: <.0001
## Target Significance Level:   0.05
## 
## Fail-safe N: 393
# Orwin's Fail-safe N

fsn(yi, vi, data=dat, type="Orwin", target=log(1.05))
## Fail-safe N Calculation Using the Orwin Approach
## 
## Average Effect Size: 0.1858
## Target Effect Size:  0.0488
## 
## Fail-safe N: 104
# trim and fill analysis

taf <- trimfill(res)
taf
## 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
predict(taf, transf=exp, digits=3)
##   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")

# cumulative meta-analysis

rcm <- cumul(res, order=vi)
rcm
##                    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

# copy data into 'dat' and examine data

dat <- dat.vanhowe1999
dat
##           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)
# proper meta-analysis using a random-effects model

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.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
predict(res, transf=exp, digits=2)
##  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
binom.test(table(sign(dat$yi)))
## 
##  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
# compute one-sided p-values

dat$pval <- pnorm(dat$yi / sqrt(dat$vi))
dat
##           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)
}
# apply Fisher's method for pooling the p-values

fisher(dat$pval)
## 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
# apply Stouffer's method for pooling the p-values

stouffer(dat$pval)
## 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
# Mantel-Haenszel method

rma.mh(measure="OR", ai=events1, n1i=n1, ci=events2, n2i=n2, data=dat)
## 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
# Peto's method

rma.peto(ai=events1, n1i=n1, ci=events2, n2i=n2, data=dat)
## 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
# note: minor difference probably due to rounding in the book

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
predict(res, transf=exp, digits=3)
##   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)

# cumulative meta-analysis

rcm <- cumul(res)
# 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.