The Handbook of Research Synthesis and Meta-Analysis by
Cooper et al. (2019), now in its third edition, has been one of the
quintessential texts on meta-analysis and the entire research synthesis
process as a whole. In this document, I provide the R code to reproduce
the worked examples and analyses from various chapters. Emphasis will be
on using the metafor
package, but several other packages
will also be used. To read more about the metafor
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:
To recreate Figure 1.1 (showing the number of citations to articles including the terms ‘research synthesis’, ‘systematic review’, or ‘meta-analysis’ in their titles), I redid the search in the Web of Science Core Collection, which yielded the following data:
dat <- read.table(header=TRUE, text = "
year hits
2020 15973
2019 32120
2018 26293
2017 23885
2016 21318
2015 18487
2014 14781
2013 12357
2012 9802
2011 7528
2010 6120
2009 5121
2008 4006
2007 3553
2006 2771
2005 2336
2004 1911
2003 1526
2002 1309
2001 1005
2000 891
1999 832
1998 729
1997 580
1996 466
1995 70
1994 25
1993 11
1992 5
1991 9
1990 20
1989 127
1988 104
")
dat <- dat[-1,] # remove current year (not complete)
dat <- dat[dat$year >= 1995,] # keep only 1995 or later
dat <- dat[nrow(dat):1,] # reverse order
We can then create a bar chart based on these data:
# Figure 1.1
par(mar=c(4,4,2,2))
barplot(dat$hits, names.arg=dat$year, las=2, space=0.4, col="#6c9ece", border=NA)
abline(h=seq(0, 30000, by=5000), col="gray")
barplot(dat$hits, space=0.4, col="#6c9ece", border=NA, add=TRUE, axes=FALSE)
Part of the code for this chapter (by Jack Vevea, Nicole Zelinksy, and Robert Orwin) is adapted from the chapter itself (see section 10.5). Below, we will make use of several additional packages that need to be installed (if they are not already installed). So let’s do this first.
# install the 'irr' package (if it is not already installed) and load it
if (!require(irr)) {
install.packages("irr")
library(irr)
}
# install the 'psych' package (if it is not already installed) and load it
if (!require(psych)) {
install.packages("psych")
library(psych)
}
# install the 'vcd' package (if it is not already installed) and load it
if (!require(vcd)) {
install.packages("vcd")
library(vcd)
}
# install the 'lme4' package (if it is not already installed) and load it
if (!require(lme4)) {
install.packages("lme4")
library(lme4)
}
# Table 10.1
dat <- read.table(header=TRUE, text = "
study c1 c2 c3
1 3 2 3
2 3 1 1
3 2 2 2
4 3 2 3
5 1 1 1
6 3 1 3
7 2 2 1
8 1 1 1
9 2 2 1
10 2 1 3
11 2 2 2
12 3 3 3
13 3 1 2
14 2 1 1
15 1 1 1
16 1 1 2
17 3 3 1
18 2 2 2
19 2 2 2
20 3 1 1
21 2 1 2
22 1 1 3
23 3 2 2
24 3 3 3
25 2 2 3")
# put ratings for coders 1, 2, and 3 into separate vectors
c1 <- dat$c1
c2 <- dat$c2
c3 <- dat$c3
# combine ratings for coders 1 and 2 into a matrix
c1c2 <- cbind(c1, c2)
# combine ratings from all three coders into a matrix
all3 <- cbind(c1, c2, c3)
## c1
## c2 1 2 3 Sum
## 1 5 3 4 12
## 2 0 7 3 10
## 3 0 0 3 3
## Sum 5 10 10 25
# note: first variable is for rows, second is for columns, so to reproduce
# panel A of Table 10.2, we have to use table(c2, c1)
## [1] 0.6
## [1] 0.36
## Percentage agreement (Tolerance=0)
##
## Subjects = 25
## Raters = 2
## %-agree = 60
# note: agree(c1c2) would have been sufficient, but due to the large number of
# additional packages being used, I will make it clear by using the :: operator
# which package a function belongs to (unless this is clear from the context)
## Percentage agreement (Tolerance=0)
##
## Subjects = 25
## Raters = 3
## %-agree = 36
## Cohen's Kappa for 2 Raters (Weights: unweighted)
##
## Subjects = 25
## Raters = 2
## Kappa = 0.425
##
## z = 3.51
## p-value = 0.000455
## Fleiss' Kappa for m Raters
##
## Subjects = 25
## Raters = 3
## Kappa = 0.315
##
## z = 3.85
## p-value = 0.00012
## Cohen's Kappa for 2 Raters (Weights: 0,1,2)
##
## Subjects = 25
## Raters = 2
## Kappa = 0.386
##
## z = 3.22
## p-value = 0.00126
We can also use the psych
package to compute Cohen’s
kappa, which also provides corresponding confidence intervals.
# unweighted and weighted Cohen's kappa for coders 1 and 2
W <- outer(1:3, 1:3, FUN = function(x,y) abs(x-y)) # create weight matrix
W
## [,1] [,2] [,3]
## [1,] 0 1 2
## [2,] 1 0 1
## [3,] 2 1 0
## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels = levels): upper or
## lower confidence interval exceed abs(1) and set to +/- 1.
## Call: cohen.kappa1(x = x, w = w, n.obs = n.obs, alpha = alpha, levels = levels)
##
## Cohen Kappa and Weighted Kappa correlation coefficients and confidence boundaries
## lower estimate upper
## unweighted kappa 0.185 0.425 0.666
## weighted kappa -1.000 0.386 1.000
##
## Number of subjects = 25
Note that the CI for weighted kappa is not correct! Using the
vcd
package, we can also compute Cohen’s kappa and obtain
the correct CI for weighted kappa.
## value ASE z Pr(>|z|) lower upper
## Unweighted 0.425 0.123 3.46 0.000532 0.185 0.666
## Weighted 0.386 0.128 3.02 0.002558 0.135 0.637
# note: the (default) weighting scheme used for computing weighted kappa by
# the function is the one described in the chapter
# Krippendorff's alpha for coders 1 and 2 when treating the data as ratings on
# a nominal, on an ordinal, or on a ratio scale
irr::kripp.alpha(t(c1c2))
## Krippendorff's alpha
##
## Subjects = 25
## Raters = 2
## alpha = 0.403
## Krippendorff's alpha
##
## Subjects = 25
## Raters = 2
## alpha = 0.278
## Krippendorff's alpha
##
## Subjects = 25
## Raters = 2
## alpha = 0.311
## [1] 0.4520228
## Mean of bivariate correlations R
##
## Subjects = 25
## Raters = 3
## R = 0.331
##
## z = 1.55
## p-value = 0.121
## Call: psych::ICC(x = c1c2)
##
## Intraclass correlation coefficients
## type ICC F df1 df2 p lower bound upper bound
## Single_raters_absolute ICC1 0.28 1.8 24 25 0.082 -0.118 0.60
## Single_random_raters ICC2 0.35 2.6 24 24 0.010 -0.025 0.65
## Single_fixed_raters ICC3 0.45 2.6 24 24 0.010 0.075 0.71
## Average_raters_absolute ICC1k 0.43 1.8 24 25 0.082 -0.268 0.75
## Average_random_raters ICC2k 0.52 2.6 24 24 0.010 -0.050 0.79
## Average_fixed_raters ICC3k 0.62 2.6 24 24 0.010 0.140 0.83
##
## Number of subjects = 25 Number of Judges = 2
## See the help file for a discussion of the other 4 McGraw and Wong estimates,
# note: this function computes 6 different types of ICCs; the first three are
# discussed in the chapter and correspond to the three different designs
# described on page 187
Using the lmer()
function from the lme4
package, we can also do these calculations manually.
# restructure data into 'long' format
dat <- data.frame(study = 1:25,
rater = rep(1:2, each=25),
rating = c(c1,c2))
# absolute agreement based on one-way random-effects model
res <- lmer(rating ~ (1 | study), data = dat)
vcs <- data.frame(VarCorr(res))
vcs$vcov[1] / (vcs$vcov[1] + vcs$vcov[2])
## [1] 0.2777018
# absolute agreement based on two-way random-effects model
res <- lmer(rating ~ (1 | study) + (1 | rater), data = dat)
vcs <- data.frame(VarCorr(res))
vcs$vcov[1] / (vcs$vcov[1] + vcs$vcov[2] + vcs$vcov[3])
## [1] 0.3545258
# absolute agreement based on two-way mixed-effects model
res <- lmer(rating ~ rater + (1 | study), data = dat)
vcs <- data.frame(VarCorr(res))
vcs$vcov[1] / (vcs$vcov[1] + vcs$vcov[2])
## [1] 0.4503106
# example data from page 199
dat <- data.frame(
study = 1:25,
rater = rep(1:3, each=25),
rating = c(3,3,2,3,NA,3,2,1,2,2,NA,3,3,2,1,1,3,2,2,3,2,1,3,NA,2,
2,1,NA,2,1,1,2,1,2,1,2,3,1,1,NA,1,3,2,2,1,1,1,2,3,2,
3,1,2,3,1,3,1,1,NA,3,2,3,2,1,1,2,1,2,2,1,2,3,2,3,3))
dat[c(1:4, 71:75),]
## study rater rating
## 1 1 1 3
## 2 2 1 3
## 3 3 1 2
## 4 4 1 3
## 71 21 3 2
## 72 22 3 3
## 73 23 3 2
## 74 24 3 3
## 75 25 3 3
# absolute agreement for all three raters (based on one-way random-effects model)
res <- lmer(rating ~ (1 | study), data = dat)
vcs <- data.frame(VarCorr(res))
vcs$vcov[1] / (vcs$vcov[1] + vcs$vcov[2])
## [1] 0.1794355
# data for Figure 11.1
dat <- read.table(header=TRUE, text = "
study md n var se pval
A 0.400 60 0.067 0.258 0.121
B 0.200 600 0.007 0.082 0.014
C 0.300 100 0.040 0.201 0.134
D 0.400 200 0.020 0.141 0.005
E 0.300 400 0.010 0.100 0.003
F -0.200 200 0.020 0.141 0.157")
dat
## study md n var se pval
## 1 A 0.4 60 0.067 0.258 0.121
## 2 B 0.2 600 0.007 0.082 0.014
## 3 C 0.3 100 0.040 0.201 0.134
## 4 D 0.4 200 0.020 0.141 0.005
## 5 E 0.3 400 0.010 0.100 0.003
## 6 F -0.2 200 0.020 0.141 0.157
# Figure 11.1
res <- rma(md, var, data=dat, method="EE", slab=study)
tmp <- dat[-1]
tmp$se <- fmtx(tmp$se, digits=3)
tmp$var <- fmtx(tmp$var, digits=3)
size <- sqrt(weights(res))
size <- 2.5 * size / max(size)
par(mar=c(5,4,2,2))
forest(res, xlim=c(-6.5,1.4), psize=size, header=TRUE, mlab="Combined",
efac=c(0,1,2), annotate=FALSE, xlab="Standardized Mean Difference",
ilab=tmp, ilab.xpos=c(-5.0, -4.1, -3.2, -2.3, -1.4))
text(-5.0, 8, "Mean\nDifference", font=2)
text(-4.1, 8, "Sample\nSize", font=2)
text(-3.2, 8, "Variance", font=2)
text(-2.3, 8, "Standard\nError", font=2)
text(-1.4, 8, "p-Value", font=2)
# 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
## 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)
summary(dat)
## yi vi sei zi pval ci.lb ci.ub
## 1 0.5924 0.0418 0.2043 2.8993 0.0037 0.1919 0.9929
# 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 exact same equation as
# given in the book is used; but the difference is usually negligible (the results below
# differ slightly from those given in the book due to intermittent rounding in the book)
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
# 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)
summary(dat)
## yi vi sei zi pval ci.lb ci.ub
## 1 0.4160 0.0137 0.1172 3.5502 0.0004 0.1863 0.6457
# 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 exact same equation as
# given in the book is used; but the difference is usually negligible (the results below
# differ slightly from those given in the book due to intermittent rounding in the book
# and also because the equation given in the book contains a mistake)
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
# standardized mean difference based on ANCOVA results
# note: not implemented in metafor, so we have to do the computations manually
Sw <- 5.5 / sqrt(1 - 0.7^2)
d <- (103 - 100) / Sw
Vd <- (50 + 50) * (1 - 0.7^2) / (50 * 50) + d^2 / (2*(50 + 50 - 2 - 1))
J <- metafor:::.cmicalc(50 + 50 - 2 - 1)
g <- J * d
Vg <- J^2 * Vd
round(g, digits=4)
## [1] 0.3865
## [1] 0.0209
## [1] 0.1444
## yi vi sei zi pval ci.lb ci.ub
## 1 0.5493 0.0103 0.1015 5.4100 <.0001 0.3503 0.7483
## [1] 0.5
## yi vi sei zi pval ci.lb ci.ub
## 1 -0.0500 0.0014 0.0371 -1.3484 0.1775 -0.1227 0.0227
## 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
# odds ratio (log transformed) for a case-control study
dat <- escalc("OR", ai=25, bi=20, ci=75, di=80)
summary(dat)
## yi vi sei zi pval ci.lb ci.ub
## 1 0.2877 0.1158 0.3403 0.8453 0.3980 -0.3794 0.9547
# Table 12.1: Data for the Gender Differences in Conformity Example
dat <- read.table(header=TRUE, text = "
study group stdingrp nitems pmaleauth n d v
1 1 1 2 141 25 -0.330 0.029
2 1 2 2 119 25 0.070 0.034
3 2 1 2 191 50 -0.300 0.022
4 3 1 38 254 100 0.350 0.016
5 3 2 30 64 100 0.700 0.066
6 3 3 45 20 100 0.850 0.218
7 3 4 45 90 100 0.400 0.045
8 3 5 45 60 100 0.480 0.069
9 3 6 5 80 100 0.370 0.051
10 3 7 5 125 100 -0.060 0.032")
# note: including the 'percent male authors' variable from Table 12.3
dat
## study group stdingrp nitems pmaleauth n d v
## 1 1 1 1 2 141 25 -0.33 0.029
## 2 2 1 2 2 119 25 0.07 0.034
## 3 3 2 1 2 191 50 -0.30 0.022
## 4 4 3 1 38 254 100 0.35 0.016
## 5 5 3 2 30 64 100 0.70 0.066
## 6 6 3 3 45 20 100 0.85 0.218
## 7 7 3 4 45 90 100 0.40 0.045
## 8 8 3 5 45 60 100 0.48 0.069
## 9 9 3 6 5 80 100 0.37 0.051
## 10 10 3 7 5 125 100 -0.06 0.032
## Equal-Effects Model (k = 10)
##
## I^2 (total heterogeneity / total variability): 71.68%
## H^2 (total variability / sampling variability): 3.53
##
## Test for Heterogeneity:
## Q(df = 9) = 31.776, p-val < .001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.124 0.060 2.074 0.038 0.007 0.241 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Note: The book uses the term ‘fixed-effects model’, but I will use the term ‘equal-effects model’ throughout this document. My reasons for using the latter term are explained here.
## Random-Effects Model (k = 10; tau^2 estimator: DL)
##
## tau^2 (estimated amount of total heterogeneity): 0.094 (SE = 0.066)
## tau (square root of estimated tau^2 value): 0.307
## I^2 (total heterogeneity / total variability): 71.68%
## H^2 (total variability / sampling variability): 3.53
##
## Test for Heterogeneity:
## Q(df = 9) = 31.776, p-val < .001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.189 0.118 1.594 0.111 -0.043 0.421
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# note: unfortunately, the estimate of tau^2 was not computed correctly in the
# book (c=242.1138, not 269.798) and hence all of the results given in the left
# column on page 251 are incorrect
# fixed-effects ANOVA-type analysis
res <- rma(d, v, mods = ~ factor(group) - 1, data=dat, method="FE")
print(res, digits=3)
## Fixed-Effects with Moderators Model (k = 10)
##
## I^2 (residual heterogeneity / unaccounted variability): 36.85%
## H^2 (unaccounted variability / sampling variability): 1.58
##
## Test for Residual Heterogeneity:
## QE(df = 7) = 11.084, p-val = 0.135
##
## Test of Moderators (coefficients 1:3):
## QM(df = 3) = 24.992, p-val < .001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## factor(group)1 -0.146 0.125 -1.166 0.244 -0.391 0.099
## factor(group)2 -0.300 0.148 -2.023 0.043 -0.591 -0.009 *
## factor(group)3 0.339 0.077 4.421 <.001 0.189 0.490 ***
##
## ---
## 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 effect for the three groups
## Fixed-Effects Model (k = 3)
##
## I^2 (total heterogeneity / total variability): 90.33%
## H^2 (total variability / sampling variability): 10.35
##
## Test for Heterogeneity:
## Q(df = 2) = 20.692, p-val < .001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.124 0.060 2.074 0.038 0.007 0.241 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# partitioning of the Q-statistics
res <- rma(d, v, mods = ~ factor(group), data=dat, method="FE")
res
## Fixed-Effects with Moderators Model (k = 10)
##
## I^2 (residual heterogeneity / unaccounted variability): 36.85%
## H^2 (unaccounted variability / sampling variability): 1.58
## R^2 (amount of heterogeneity accounted for): 55.15%
##
## Test for Residual Heterogeneity:
## QE(df = 7) = 11.0843, p-val = 0.1350
##
## Test of Moderators (coefficients 2:3):
## QM(df = 2) = 20.6916, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt -0.1459 0.1251 -1.1660 0.2436 -0.3911 0.0993
## factor(group)2 -0.1541 0.1940 -0.7943 0.4270 -0.5344 0.2262
## factor(group)3 0.4851 0.1468 3.3053 0.0009 0.1975 0.7728 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# not removing the intercept, so the QM-statistic is equal to Q-between
round(res$QM, digits=3) # Q-between
## [1] 20.692
## [1] 11.084
# Q-within for each group
res1 <- rma(d, v, data=dat, method="FE", subset=group==1)
res2 <- rma(d, v, data=dat, method="FE", subset=group==2)
res3 <- rma(d, v, data=dat, method="FE", subset=group==3)
round(res1$QE, digits=3)
## [1] 2.54
## [1] 0
## [1] 8.545
## [1] 11.084
# contrast between group 1 and 3
res <- rma(d, v, mods = ~ factor(group) - 1, data=dat, method="FE")
anova(res, X=c(-1,0,1), digits=3)
## Hypothesis:
## 1: -factor(group)1 + factor(group)3 = 0
##
## Results:
## estimate se zval pval
## 1: 0.485 0.147 3.305 <.001
##
## Test of Hypothesis:
## QM(df = 1) = 10.925, p-val < .001
## pred se ci.lb ci.ub
## 0.485 0.147 0.197 0.773
# mixed-effects model analysis
# distribution-free (method of moments) estimate of tau^2
res <- rma(d, v, mods = ~ factor(group), data=dat, method="DL")
round(res$tau2, digits=3)
## [1] 0.025
# maximum likelihood estimate of tau^2
res <- rma(d, v, mods = ~ factor(group), data=dat, method="ML")
round(res$tau2, digits=3)
## [1] 0
# restricted maximum likelihood estimate of tau^2
res <- rma(d, v, mods = ~ factor(group), data=dat, method="REML")
round(res$tau2, digits=3)
## [1] 0.027
## Mixed-Effects Model (k = 10; tau^2 estimator: DL)
##
## tau^2 (estimated amount of residual heterogeneity): 0.025 (SE = 0.037)
## tau (square root of estimated tau^2 value): 0.157
## I^2 (residual heterogeneity / unaccounted variability): 36.85%
## H^2 (unaccounted variability / sampling variability): 1.58
##
## Test for Residual Heterogeneity:
## QE(df = 7) = 11.084, p-val = 0.135
##
## Test of Moderators (coefficients 1:3):
## QM(df = 3) = 15.062, p-val = 0.002
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## factor(group)1 -0.139 0.168 -0.829 0.407 -0.467 0.190
## factor(group)2 -0.300 0.216 -1.387 0.165 -0.724 0.124
## factor(group)3 0.361 0.102 3.528 <.001 0.161 0.562 ***
##
## ---
## 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 effect for the three groups
# weighted grand mean effect size
rma(coef(res), diag(vcov(res)), method="FE", digits=3)
## Fixed-Effects Model (k = 3)
##
## I^2 (total heterogeneity / total variability): 82.69%
## H^2 (total variability / sampling variability): 5.78
##
## Test for Heterogeneity:
## Q(df = 2) = 11.557, p-val = 0.003
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.152 0.081 1.872 0.061 -0.007 0.310 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Hypothesis:
## 1: -factor(group)1 + factor(group)3 = 0
##
## Results:
## estimate se zval pval
## 1: 0.500 0.196 2.547 0.011
##
## Test of Hypothesis:
## QM(df = 1) = 6.485, p-val = 0.011
## pred se ci.lb ci.ub pi.lb pi.ub
## 0.500 0.196 0.115 0.885 0.007 0.993
# meta-regression model
res <- rma(d, v, mods = ~ I(log(nitems)), data=dat, method="FE")
print(res, digits=3)
## Fixed-Effects with Moderators Model (k = 10)
##
## I^2 (residual heterogeneity / unaccounted variability): 7.66%
## H^2 (unaccounted variability / sampling variability): 1.08
## R^2 (amount of heterogeneity accounted for): 69.33%
##
## Test for Residual Heterogeneity:
## QE(df = 8) = 8.664, p-val = 0.371
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 23.112, p-val < .001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt -0.323 0.111 -2.922 0.003 -0.540 -0.106 **
## I(log(nitems)) 0.210 0.044 4.808 <.001 0.125 0.296 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# note: when doing transformations on predictors (such as taking the log) in
# the model formula, then we need to wrap this inside the I() function
# mixed-effects meta-regression model
# distribution-free (method of moments) estimate of tau^2
res <- rma(d, v, mods = ~ I(log(nitems)), data=dat, method="DL")
round(res$tau2, digits=3)
## [1] 0.003
# maximum likelihood estimate of tau^2
res <- rma(d, v, mods = ~ I(log(nitems)), data=dat, method="ML")
round(res$tau2, digits=3)
## [1] 0
# restricted maximum likelihood estimate of tau^2
res <- rma(d, v, mods = ~ I(log(nitems)), data=dat, method="REML")
round(res$tau2, digits=3)
## [1] 0.004
# continuing with the distribution-free (method of moments) estimate of tau^2
res <- rma(d, v, mods = ~ I(log(nitems)), data=dat, method="DL")
print(res, digits=3)
## Mixed-Effects Model (k = 10; tau^2 estimator: DL)
##
## tau^2 (estimated amount of residual heterogeneity): 0.003 (SE = 0.021)
## tau (square root of estimated tau^2 value): 0.057
## I^2 (residual heterogeneity / unaccounted variability): 7.66%
## H^2 (unaccounted variability / sampling variability): 1.08
## R^2 (amount of heterogeneity accounted for): 96.53%
##
## Test for Residual Heterogeneity:
## QE(df = 8) = 8.664, p-val = 0.371
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 20.821, p-val < .001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt -0.321 0.117 -2.744 0.006 -0.550 -0.092 **
## I(log(nitems)) 0.212 0.046 4.563 <.001 0.121 0.303 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Mixed-Effects Model (k = 10; tau^2 estimator: DL)
##
## tau^2 (estimated amount of residual heterogeneity): 0.003 (SE = 0.021)
## tau (square root of estimated tau^2 value): 0.057
## I^2 (residual heterogeneity / unaccounted variability): 7.66%
## H^2 (unaccounted variability / sampling variability): 1.08
## R^2 (amount of heterogeneity accounted for): 96.53%
##
## Test for Residual Heterogeneity:
## QE(df = 8) = 8.664, p-val = 0.371
##
## Number of estimates: 10
## Number of clusters: 10
## Estimates per cluster: 1
##
## Test of Moderators (coefficient 2):¹
## F(df1 = 1, df2 = 8) = 25.315, p-val = 0.001
##
## Model Results:
##
## estimate se¹ tval¹ df¹ pval¹ ci.lb¹ ci.ub¹
## intrcpt -0.321 0.124 -2.586 8 0.032 -0.607 -0.035 *
## I(log(nitems)) 0.212 0.042 5.031 8 0.001 0.115 0.309 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## 1) results based on cluster-robust inference (var-cov estimator: CR1,
## approx t/F-tests and confidence intervals, df: residual method)
# Table 13.1
dat <- dat.kalaian1996[c(1:8,17:29,9:16,30:31),c(2,4:5,7,6,8)]
dat$study <- rep(1:26, times=rle(dat$study)$lengths)
names(dat) <- c("study", "nt", "nc", "d", "outcome", "v")
dat
## study nt nc d outcome v
## 1 1 28 22 0.22 verbal 0.0817
## 2 2 39 40 0.09 verbal 0.0507
## 3 3 22 17 0.14 verbal 0.1045
## 4 4 48 43 0.14 verbal 0.0442
## 5 5 25 74 -0.01 verbal 0.0535
## 6 6 37 35 0.14 verbal 0.0557
## 7 7 24 70 0.18 verbal 0.0561
## 8 8 16 19 0.01 verbal 0.1151
## 17 9 43 37 0.01 verbal 0.0503
## 18 10 19 13 0.67 verbal 0.1366
## 19 11 16 11 -0.38 verbal 0.1561
## 20 12 20 12 -0.24 verbal 0.1342
## 21 13 39 28 0.29 verbal 0.0620
## 22 14 38 25 0.26 math 0.0669
## 23 15 18 13 -0.41 math 0.1352
## 24 16 19 13 0.08 math 0.1297
## 25 17 37 22 0.30 math 0.0732
## 26 18 19 11 -0.53 math 0.1482
## 27 19 17 13 0.13 math 0.1360
## 28 20 20 12 0.26 math 0.1344
## 29 21 20 13 0.47 math 0.1303
## 9 22 145 129 0.13 verbal 0.0147
## 10 22 145 129 0.12 math 0.0147
## 11 23 72 129 0.25 verbal 0.0218
## 12 23 72 129 0.06 math 0.0216
## 13 24 71 129 0.31 verbal 0.0221
## 14 24 71 129 0.09 math 0.0219
## 15 25 13 14 0.00 verbal 0.1484
## 16 25 13 14 0.07 math 0.1484
## 30 26 16 17 0.13 verbal 0.1216
## 31 26 16 17 0.48 math 0.1248
# note: this is a subset of the studies in dat.kalaian1996; instead of using
# dummy variable 'x' (where x=0 for the verbal and x=1 for the math subtest),
# we use the 'outcome' variable that has more descriptive labels
# construct variance-covariance matrices assuming rho = 0.7
V <- vcalc(v, cluster=study, type=outcome, rho=0.7, data=dat)
# var-cov matrix for studies 22 to 26
blsplit(V, cluster=dat$study, round, 4)[22:26]
## $`22`
## [,1] [,2]
## [1,] 0.0147 0.0103
## [2,] 0.0103 0.0147
##
## $`23`
## [,1] [,2]
## [1,] 0.0218 0.0152
## [2,] 0.0152 0.0216
##
## $`24`
## [,1] [,2]
## [1,] 0.0221 0.0154
## [2,] 0.0154 0.0219
##
## $`25`
## [,1] [,2]
## [1,] 0.1484 0.1039
## [2,] 0.1039 0.1484
##
## $`26`
## [,1] [,2]
## [1,] 0.1216 0.0862
## [2,] 0.0862 0.1248
# multivariate model with (correlated) random effects for both outcomes
res <- rma.mv(d, V, mods = ~ outcome - 1,
random = ~ outcome | study, struct="UN", data=dat)
print(res, digits=3)
## Multivariate Meta-Analysis Model (k = 31; method: REML)
##
## Variance Components:
##
## outer factor: study (nlvls = 26)
## inner factor: outcome (nlvls = 2)
##
## estim sqrt k.lvl fixed level
## tau^2.1 0.001 0.025 13 no math
## tau^2.2 0.000 0.021 18 no verbal
##
## rho.math rho.vrbl math vrbl
## math 1 - 5
## verbal -1.000 1 no -
##
## Test for Residual Heterogeneity:
## QE(df = 29) = 19.902, p-val = 0.896
##
## Test of Moderators (coefficients 1:2):
## QM(df = 2) = 9.417, p-val = 0.009
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## outcomemath 0.089 0.059 1.511 0.131 -0.026 0.204
## outcomeverbal 0.158 0.051 3.067 0.002 0.057 0.258 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# note: by removing the intercept, the two coefficients directly provide the
# estimated average effect for the two outcomes
# note: both variance components are very close to 0; as in the book, we
# proceed with a model where both variances are constrained to be equal to 0
# multivariate fixed-effects model
res <- rma.mv(d, V, mods = ~ outcome - 1, data=dat)
print(res, digits=3)
## Multivariate Meta-Analysis Model (k = 31; method: REML)
##
## Variance Components: none
##
## Test for Residual Heterogeneity:
## QE(df = 29) = 19.902, p-val = 0.896
##
## Test of Moderators (coefficients 1:2):
## QM(df = 2) = 9.813, p-val = 0.007
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## outcomemath 0.086 0.057 1.507 0.132 -0.026 0.198
## outcomeverbal 0.159 0.051 3.131 0.002 0.059 0.258 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# fit model with varying values of rho
rhos <- c(0, 0.5, 0.6, 0.7, 0.8)
res <- list()
for (i in 1:length(rhos)) {
V <- vcalc(v, cluster=study, type=outcome, rho=rhos[i], data=dat)
res[[i]] <- rma.mv(d, V, mods = ~ outcome - 1, data=dat)
}
# Table 13.2
tab <- data.frame(rho = rhos,
b1 = sapply(res, function(x) coef(x)[1]),
s1 = sapply(res, function(x) x$se[1]),
lcl1 = sapply(res, function(x) x$ci.lb[1]),
ucl1 = sapply(res, function(x) x$ci.ub[1]),
z1 = sapply(res, function(x) x$zval[1]),
b2 = sapply(res, function(x) coef(x)[2]),
s2 = sapply(res, function(x) x$se[2]),
lcl2 = sapply(res, function(x) x$ci.lb[2]),
ucl2 = sapply(res, function(x) x$ci.ub[2]),
z2 = sapply(res, function(x) x$zval[2]))
dfround(tab, digits=c(1,2,3,2,2,2,2,3,2,2,2))
## rho b1 s1 lcl1 ucl1 z1 b2 s2 lcl2 ucl2 z2
## 1 0.0 0.11 0.064 -0.01 0.24 1.77 0.15 0.053 0.05 0.26 2.93
## 2 0.5 0.09 0.060 -0.02 0.21 1.56 0.16 0.052 0.06 0.26 3.03
## 3 0.6 0.09 0.059 -0.03 0.21 1.53 0.16 0.051 0.06 0.26 3.08
## 4 0.7 0.09 0.057 -0.03 0.20 1.51 0.16 0.051 0.06 0.26 3.13
## 5 0.8 0.08 0.055 -0.03 0.19 1.49 0.16 0.050 0.06 0.26 3.20
# note: there are some printing errors in the table in the book;
# also the values given in the Z_2 column in the book are quite off
# robust variance estimation
V <- vcalc(v, cluster=study, type=outcome, rho=0.7, data=dat)
res <- rma.mv(d, V, mods = ~ outcome - 1, data=dat)
robust(res, cluster=study)
## Multivariate Meta-Analysis Model (k = 31; method: REML)
##
## Variance Components: none
##
## Test for Residual Heterogeneity:
## QE(df = 29) = 19.9016, p-val = 0.8960
##
## Number of estimates: 31
## Number of clusters: 26
## Estimates per cluster: 1-2 (mean: 1.19, median: 1)
##
## Test of Moderators (coefficients 1:2):¹
## F(df1 = 2, df2 = 24) = 12.1231, p-val = 0.0002
##
## Model Results:
##
## estimate se¹ tval¹ df¹ pval¹ ci.lb¹ ci.ub¹
## outcomemath 0.0861 0.0383 2.2490 24 0.0340 0.0071 0.1651 *
## outcomeverbal 0.1589 0.0362 4.3950 24 0.0002 0.0843 0.2336 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## 1) results based on cluster-robust inference (var-cov estimator: CR1,
## approx t/F-tests and confidence intervals, df: residual method)
# note: this is not the exact same approach that is described in the book; it
# uses a different weighting scheme and uses k-p for the degrees of freedom
# for the t-statistics; we can make use of the 'effective degrees of freedom'
# method described in the book with the 'clubSandwich' package
# install the 'clubSandwich' package (if it is not already installed) and load it
if (!require(clubSandwich))
install.packages("clubSandwich")
## Multivariate Meta-Analysis Model (k = 31; method: REML)
##
## Variance Components: none
##
## Test for Residual Heterogeneity:
## QE(df = 29) = 19.9016, p-val = 0.8960
##
## Number of estimates: 31
## Number of clusters: 26
## Estimates per cluster: 1-2 (mean: 1.19, median: 1)
##
## Test of Moderators (coefficients 1:2):¹
## F(df1 = 2, df2 = 5.28) = 10.5212, p-val = 0.0144
##
## Model Results:
##
## estimate se¹ tval¹ df¹ pval¹ ci.lb¹ ci.ub¹
## outcomemath 0.0861 0.0395 2.1819 6.67 0.0674 -0.0082 0.1803 .
## outcomeverbal 0.1589 0.0370 4.2915 11.03 0.0013 0.0774 0.2404 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## 1) results based on cluster-robust inference (var-cov estimator: CR2,
## approx t/F-tests and confidence intervals, df: Satterthwaite approx)
# note: but these are still not the same results as given in the book because
# of the different weighting scheme; we can get the same results as given in
# the book with the 'robumeta' package
# install the 'robumeta' package (if it is not already installed) and load it
if (!require(robumeta)) {
install.packages("robumeta")
library(robumeta)
}
## RVE: Correlated Effects Model with Small-Sample Corrections
##
## Model: d ~ outcome - 1
##
## Number of studies = 26
## Number of outcomes = 31 (min = 1 , mean = 1.19 , median = 1 , max = 2 )
## Rho = 0.7
## I.sq = 0
## Tau.sq = 0
##
## Estimate StdErr t-value dfs P(|t|>) 95% CI.L 95% CI.U Sig
## 1 outcomemath 0.114 0.0489 2.32 7.69 0.05018 -0.000115 0.227 *
## 2 outcomeverbal 0.140 0.0338 4.12 13.09 0.00119 0.066450 0.213 ***
## ---
## Signif. codes: < .01 *** < .05 ** < .10 *
## ---
## Note: If df < 4, do not trust the results
# reproduce robumeta results using metafor and clubSandwich
vcalcfun <- function(v, tau2)
diag((tau2 + mean(v)) * length(v), nrow=length(v), ncol=length(v))
V <- lapply(split(dat$v, dat$study), vcalcfun, tau2=0)
res <- rma.mv(d, V, mods = ~ outcome - 1, data=dat)
robust(res, cluster=study, clubSandwich=TRUE)
## Multivariate Meta-Analysis Model (k = 31; method: REML)
##
## Variance Components: none
##
## Test for Residual Heterogeneity:
## QE(df = 29) = 14.6868, p-val = 0.9873
##
## Number of estimates: 31
## Number of clusters: 26
## Estimates per cluster: 1-2 (mean: 1.19, median: 1)
##
## Test of Moderators (coefficients 1:2):¹
## F(df1 = 2, df2 = 10.54) = 11.0146, p-val = 0.0026
##
## Model Results:
##
## estimate se¹ tval¹ df¹ pval¹ ci.lb¹ ci.ub¹
## outcomemath 0.1135 0.0489 2.3198 7.69 0.0502 -0.0001 0.2272 .
## outcomeverbal 0.1395 0.0338 4.1220 13.09 0.0012 0.0664 0.2126 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## 1) results based on cluster-robust inference (var-cov estimator: CR2,
## approx t/F-tests and confidence intervals, df: Satterthwaite approx)
# construct dataset with synthetic effects
agg <- aggregate(dat, by=list(dat$study), function(x) {
if (is.character(x))
paste(unique(x), collapse="/")
else
mean(x)
})
agg$Group.1 <- agg$x <- NULL
agg
## study nt nc d outcome v
## 1 1 28 22 0.220 verbal 0.0817
## 2 2 39 40 0.090 verbal 0.0507
## 3 3 22 17 0.140 verbal 0.1045
## 4 4 48 43 0.140 verbal 0.0442
## 5 5 25 74 -0.010 verbal 0.0535
## 6 6 37 35 0.140 verbal 0.0557
## 7 7 24 70 0.180 verbal 0.0561
## 8 8 16 19 0.010 verbal 0.1151
## 9 9 43 37 0.010 verbal 0.0503
## 10 10 19 13 0.670 verbal 0.1366
## 11 11 16 11 -0.380 verbal 0.1561
## 12 12 20 12 -0.240 verbal 0.1342
## 13 13 39 28 0.290 verbal 0.0620
## 14 14 38 25 0.260 math 0.0669
## 15 15 18 13 -0.410 math 0.1352
## 16 16 19 13 0.080 math 0.1297
## 17 17 37 22 0.300 math 0.0732
## 18 18 19 11 -0.530 math 0.1482
## 19 19 17 13 0.130 math 0.1360
## 20 20 20 12 0.260 math 0.1344
## 21 21 20 13 0.470 math 0.1303
## 22 22 145 129 0.125 verbal/math 0.0147
## 23 23 72 129 0.155 verbal/math 0.0217
## 24 24 71 129 0.200 verbal/math 0.0220
## 25 25 13 14 0.035 verbal/math 0.1484
## 26 26 16 17 0.305 verbal/math 0.1232
# fit standard random-effects model to the dataset with synthetic effects
res <- rma(d, v, data=agg, method="DL")
print(res, digits=3)
## Random-Effects Model (k = 26; tau^2 estimator: DL)
##
## tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.018)
## 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 = 25) = 13.530, p-val = 0.969
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.130 0.048 2.699 0.007 0.036 0.224 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
For these analyses, we could do most of what is described in the chapter using the excellent brms package. However, to fully reproduce all of the analyses conducted, we either need to use WinBUGS (as was done by the chapter authors) or we can use JAGS (which has the advantage that it runs not just under Windows, but also macOS / Mac OS X and Linux). Here, we will make use of the latter, which we can interact with directly from R via the rjags package. Note that JAGS needs to be installed separately (follow the link above for installation instructions).
# install the 'rjags' package (if it is not already installed) and load it
if (!require(rjags)) {
install.packages("rjags")
library(rjags)
}
# Table 14.1: Respiratory Tract Infections Data
dat <- dat.damico2009[c(1:8,16,9:15),-8]
dat$conceal <- 1 - dat$conceal
rownames(dat) <- 1:nrow(dat)
dat <- escalc(measure="OR", ai=xt, n1i=nt, ci=xc, n2i=nc, data=dat, digits=2)
dat
## study year xt nt xc nc conceal yi vi
## 1 Abele-Horn 1997 13 58 23 30 1 -2.43 0.29
## 2 Aerdts 1991 1 28 29 60 0 -3.23 1.10
## 3 Blair 1991 12 161 38 170 0 -1.27 0.12
## 4 Boland 1991 14 32 17 32 1 -0.38 0.25
## 5 Cockerill 1992 4 75 12 75 1 -1.22 0.36
## 6 Finch 1991 4 20 7 24 0 -0.50 0.51
## 7 Jacobs 1992 0 45 4 46 1 -2.27 2.27
## 8 Kerver 1988 5 49 31 47 1 -2.84 0.32
## 9 Krueger 2002 91 265 149 262 0 -0.92 0.03
## 10 Palomar 1997 10 50 25 49 1 -1.43 0.21
## 11 Rocha 1992 7 47 25 54 1 -1.59 0.24
## 12 Sanchez-Garcia 1992 32 131 60 140 0 -0.84 0.07
## 13 Stoutenbeek 2007 62 201 100 200 0 -0.81 0.04
## 14 Ulrich 1989 7 55 26 57 0 -1.75 0.23
## 15 Verwaest 1997 22 193 40 185 0 -0.76 0.08
## 16 Winter 1992 3 91 17 92 0 -1.89 0.42
# note: including 'conceal' variable (coded 0 = allocation concealment
# adequate, 1 = allocation concealment inadequate)
# Bayesian equal-effects model
k <- length(dat$yi)
jags.data <- list(yi=dat$yi, vi=dat$vi, k=k)
fe.model <- "model {
for (i in 1:k) {
yi[i] ~ dnorm(theta, 1/vi[i])
}
theta ~ dnorm(0, 0.0001)
}"
inits <- list(theta=0, .RNG.name="base::Mersenne-Twister")
inits <- list(inits, inits, inits)
inits[[1]]$.RNG.seed <- 12341
inits[[2]]$.RNG.seed <- 12342
inits[[3]]$.RNG.seed <- 12343
# note: usually 'inits <- list(theta=0)' would have been sufficient; however,
# for full reproducibility, I am setting the random number generator and
# the seed for each chain (the same will be done below for other models)
model <- jags.model(textConnection(fe.model), inits=inits, data=jags.data, n.chains=3, quiet=TRUE)
update(model, n.iter=10000, progress.bar="none")
bma.ee <- coda.samples(model, variable.names=c("theta"), n.iter=100000, progress.bar="none")
dic.ee <- dic.samples(model, n.iter=100000, type="pD", progress.bar="none")
bma.ee <- summary(bma.ee)$quantiles[c(3,1,5)]
bma.ee <- rbind(theta=bma.ee, exp.theta = exp(bma.ee))
round(bma.ee, digits=2)
## 50% 2.5% 97.5%
## theta -1.09 -1.27 -0.91
## exp.theta 0.34 0.28 0.40
# compare with results from a non-Bayesian analysis
res.ee <- rma(yi, vi, data=dat, method="EE", digits=2)
predict(res.ee, transf=exp)
## pred ci.lb ci.ub
## 0.34 0.28 0.40
# Bayesian random-effects model with uniform(0,2) prior for tau
re.model <- "model {
for (i in 1:k) {
yi[i] ~ dnorm(theta[i], 1/vi[i])
theta[i] ~ dnorm(mu, 1/tau2)
}
mu ~ dnorm(0, 0.0001)
tau ~ dunif(0, 2)
tau2 <- tau^2
theta.2 <- theta[2]
theta.3 <- theta[3]
theta.new ~ dnorm(mu, 1/tau2)
}"
inits <- list(theta=rep(0,k), theta.new=0, mu=0, tau=1, .RNG.name="base::Mersenne-Twister")
inits <- list(inits, inits, inits)
inits[[1]]$.RNG.seed <- 12341
inits[[2]]$.RNG.seed <- 12342
inits[[3]]$.RNG.seed <- 12343
model <- jags.model(textConnection(re.model), inits=inits, data=jags.data, n.chains=3, quiet=TRUE)
update(model, n.iter=10000, progress.bar="none")
bma.re <- coda.samples(model, variable.names=c("mu","tau2","theta.2","theta.3","theta.new"),
n.iter=100000, progress.bar="none")
dic.re <- dic.samples(model, n.iter=100000, type="pD", progress.bar="none")
bma.re <- summary(bma.re)$quantiles[,c(3,1,5)]
bma.re <- rbind(mu=bma.re[1,], exp.mu=exp(bma.re[1,]), tau2=bma.re[2,],
bma.re[3:5,], exp.theta.new=exp(bma.re[5,]))
round(bma.re[1:3,], digits=2)
## 50% 2.5% 97.5%
## mu -1.29 -1.73 -0.94
## exp.mu 0.28 0.18 0.39
## tau2 0.27 0.01 1.09
# compare with results from a non-Bayesian analysis
res.re <- rma(yi, vi, data=dat, method="DL", digits=2)
pred.re <- predict(res.re, transf=exp)
pred.re
## pred ci.lb ci.ub pi.lb pi.ub
## 0.28 0.20 0.38 0.12 0.67
## estimate ci.lb ci.ub
## 0.18 0.04 1.20
# Bayesian random-effects model with half-normal(0,0.5^2) prior for tau
re.model <- "model {
for (i in 1:k) {
yi[i] ~ dnorm(theta[i], 1/vi[i])
theta[i] ~ dnorm(mu, 1/tau2)
}
mu ~ dnorm(0, 0.0001)
tau ~ dnorm(0, 1/0.5^2) T(0,)
tau2 <- tau^2
}"
inits <- list(theta=rep(0,k), mu=0, tau=1, .RNG.name="base::Mersenne-Twister")
inits <- list(inits, inits, inits)
inits[[1]]$.RNG.seed <- 12341
inits[[2]]$.RNG.seed <- 12342
inits[[3]]$.RNG.seed <- 12343
model <- jags.model(textConnection(re.model), inits=inits, data=jags.data, n.chains=3, quiet=TRUE)
update(model, n.iter=10000, progress.bar="none")
bma.hn <- coda.samples(model, variable.names=c("mu","tau2"), n.iter=100000, progress.bar="none")
bma.hn <- summary(bma.hn)$quantiles[,c(3,1,5)]
bma.hn <- rbind(mu=bma.hn[1,], exp.mu=exp(bma.hn[1,]), tau2=bma.hn[2,])
# Bayesian random-effects model with gamma(0.001,0.001) prior for 1/tau^2
re.model <- "model {
for (i in 1:k) {
yi[i] ~ dnorm(theta[i], 1/vi[i])
theta[i] ~ dnorm(mu, prec)
}
mu ~ dnorm(0, 0.0001)
prec ~ dgamma(0.001, 0.001)
tau2 <- 1/prec
}"
inits <- list(theta=rep(0,k), mu=0, prec=1, .RNG.name="base::Mersenne-Twister")
inits <- list(inits, inits, inits)
inits[[1]]$.RNG.seed <- 12341
inits[[2]]$.RNG.seed <- 12342
inits[[3]]$.RNG.seed <- 12343
model <- jags.model(textConnection(re.model), inits=inits, data=jags.data, n.chains=3, quiet=TRUE)
update(model, n.iter=10000, progress.bar="none")
bma.ig <- coda.samples(model, variable.names=c("mu","tau2"), n.iter=100000, progress.bar="none")
bma.ig <- summary(bma.ig)$quantiles[,c(3,1,5)]
bma.ig <- rbind(mu=bma.ig[1,], exp.mu=exp(bma.ig[1,]), tau2=bma.ig[2,])
# Table 14.2
tab <- rbind(
c(or = pred.re$pred, ci.lb = pred.re$ci.lb, ci.ub = pred.re$ci.ub,
tau2 = conf.re$random[1,1], ci.lb = conf.re$random[1,2], ci.ub = conf.re$random[1,3]),
c(bma.re[2,], bma.re[3,]), c(bma.hn[2,], bma.hn[3,]), c(bma.ig[2,], bma.ig[3,]))
rownames(tab) <- c("Frequentist RE model",
"Bayesian RE model, uniform(0,2) prior for tau",
"Bayesian RE model, half-normal(0,0.5^2) prior for tau",
"Bayesian RE model, gamma(0.001,0.001) prior for 1/tau^2")
round(tab, digits=2)
## or ci.lb ci.ub tau2 ci.lb ci.ub
## Frequentist RE model 0.28 0.20 0.38 0.18 0.04 1.20
## Bayesian RE model, uniform(0,2) prior for tau 0.28 0.18 0.39 0.27 0.01 1.09
## Bayesian RE model, half-normal(0,0.5^2) prior for tau 0.28 0.19 0.39 0.19 0.01 0.70
## Bayesian RE model, gamma(0.001,0.001) prior for 1/tau^2 0.29 0.19 0.39 0.16 0.00 0.82
# Bayesian random-effects meta-regression model with uniform(0,2) prior for tau
jags.data <- list(yi=dat$yi, vi=dat$vi, xi=dat$conceal, k=k)
me.model <- "model {
for (i in 1:k) {
yi[i] ~ dnorm(theta[i], 1/vi[i])
theta[i] ~ dnorm(beta0 + beta1 * xi[i], 1/tau2)
}
beta0 ~ dnorm(0, 0.0001)
beta1 ~ dnorm(0, 0.0001)
tau ~ dunif(0, 2)
tau2 <- tau^2
betasum <- beta0 + beta1
}"
inits <- list(theta=rep(0,k), beta0=0, beta1=1, tau=1, .RNG.name="base::Mersenne-Twister")
inits <- list(inits, inits, inits)
inits[[1]]$.RNG.seed <- 12341
inits[[2]]$.RNG.seed <- 12342
inits[[3]]$.RNG.seed <- 12343
model <- jags.model(textConnection(me.model), inits=inits, data=jags.data, n.chains=3, quiet=TRUE)
update(model, n.iter=10000, progress.bar="none")
bma.me <- coda.samples(model, variable.names=c("beta0","beta1","tau2","betasum"),
n.iter=100000, progress.bar="none")
dic.me <- dic.samples(model, n.iter=100000, type="pD", progress.bar="none")
bma.me <- summary(bma.me)$quantiles[,c(3,1,5)]
bma.me <- rbind(bma.me[c(1,2,4),], exp.beta0=exp(bma.me[1,]), exp.betasum=exp(bma.me[3,]))
round(bma.me, digits=2)
## 50% 2.5% 97.5%
## beta0 -1.05 -1.55 -0.71
## beta1 -0.58 -1.22 0.16
## tau2 0.12 0.00 0.87
## exp.beta0 0.35 0.21 0.49
## exp.betasum 0.19 0.11 0.33
# compare with results from a non-Bayesian analysis
res.me <- rma(yi, vi, mods = ~ conceal, data=dat, method="DL", digits=2)
predict(res.me, newmods=c(0,1), transf=exp)
## pred ci.lb ci.ub pi.lb pi.ub
## 1 0.34 0.25 0.48 0.16 0.74
## 2 0.20 0.12 0.32 0.08 0.45
## estimate ci.lb ci.ub
## 0.12 0.00 1.23
# Table 14.3
tab <- rbind(c(sum(dic.ee$deviance), sum(dic.ee[[2]]), sum(dic.ee$deviance) + sum(dic.ee[[2]])),
c(sum(dic.re$deviance), sum(dic.re[[2]]), sum(dic.re$deviance) + sum(dic.re[[2]])),
c(sum(dic.me$deviance), sum(dic.me[[2]]), sum(dic.me$deviance) + sum(dic.me[[2]])))
colnames(tab) <- c("mean(D)", "p_D", "DIC")
rownames(tab) <- c("Equal-effect meta-analysis",
"Random-effects meta-analysis",
"Random-effects meta-regression")
round(tab, digits=1)
## mean(D) p_D DIC
## Equal-effect meta-analysis 40.3 1.0 41.3
## Random-effects meta-analysis 24.4 9.7 34.1
## Random-effects meta-regression 26.5 8.3 34.7
## 50% 2.5% 97.5%
## theta.2 -1.62 -3.01 -0.77
## theta.3 -1.26 -1.86 -0.71
## 50% 2.5% 97.5%
## theta.new -1.27 -2.61 -0.10
## exp.theta.new 0.28 0.07 0.91
# Bayesian random-effects model for binary data
jags.data <- list(xt=dat$xt, nt=dat$nt, xc=dat$xc, nc=dat$nc, k=k)
re.model <- "model {
for (i in 1:k) {
xc[i] ~ dbin(pc[i], nc[i])
xt[i] ~ dbin(pt[i], nt[i])
logit(pc[i]) <- alpha[i] - theta[i]/2
logit(pt[i]) <- alpha[i] + theta[i]/2
theta[i] ~ dnorm(mu, 1/tau2)
alpha[i] ~ dnorm(0, 0.0001)
}
mu ~ dnorm(0, 0.0001)
tau ~ dunif(0, 2)
tau2 <- tau^2
}"
inits <- list(theta=rep(0,k), alpha=rep(0,k), mu=0, tau=1, .RNG.name="base::Mersenne-Twister")
inits <- list(inits, inits, inits)
inits[[1]]$.RNG.seed <- 12341
inits[[2]]$.RNG.seed <- 12342
inits[[3]]$.RNG.seed <- 12343
model <- jags.model(textConnection(re.model), inits=inits, data=jags.data, n.chains=3, quiet=TRUE)
update(model, n.iter=10000, progress.bar="none")
bma.re <- coda.samples(model, variable.names=c("mu","tau2"), n.iter=100000, progress.bar="none")
bma.re <- summary(bma.re)$quantiles[,c(3,1,5)]
bma.re <- rbind(mu=bma.re[1,], exp.mu=exp(bma.re[1,]), tau2=bma.re[2,])
round(bma.re, digits=2)
## 50% 2.5% 97.5%
## mu -1.40 -1.88 -1.01
## exp.mu 0.25 0.15 0.36
## tau2 0.40 0.08 1.44
# Bayesian random-effects model with log-normal(-2.49,1.52^2) prior for tau^2
jags.data <- list(yi=dat$yi, vi=dat$vi, k=k)
re.model <- "model {
for (i in 1:k) {
yi[i] ~ dnorm(theta[i], 1/vi[i])
theta[i] ~ dnorm(mu, 1/tau2)
}
mu ~ dnorm(0, 1.0E-4)
tau2 ~ dlnorm(-2.49, 1/1.52^2)
}"
inits <- list(theta=rep(0,k), mu=0, tau2=1, .RNG.name="base::Mersenne-Twister")
inits <- list(inits, inits, inits)
inits[[1]]$.RNG.seed <- 12341
inits[[2]]$.RNG.seed <- 12342
inits[[3]]$.RNG.seed <- 12343
model <- jags.model(textConnection(re.model), inits=inits, data=jags.data, n.chains=3, quiet=TRUE)
update(model, n.iter=10000, progress.bar="none")
bma.ip <- coda.samples(model, variable.names=c("mu","tau2"), n.iter=100000, progress.bar="none")
bma.ip <- summary(bma.ip)$quantiles[,c(3,1,5)]
bma.ip <- rbind(mu=bma.ip[1,], exp.mu=exp(bma.ip[1,]), tau2=bma.ip[2,])
round(bma.ip, digits=2)
## 50% 2.5% 97.5%
## mu -1.25 -1.63 -0.95
## exp.mu 0.29 0.20 0.39
## tau2 0.15 0.01 0.62
# Bayesian random-effects meta-regression model with log-normal(-2.49,1.52^2) prior for tau^2
jags.data <- list(yi=dat$yi, vi=dat$vi, xi=dat$conceal, k=k)
me.model <- "model {
for (i in 1:k) {
yi[i] ~ dnorm(theta[i], 1/vi[i])
theta[i] ~ dnorm(beta0 + beta1 * xi[i], 1/tau2)
}
beta0 ~ dnorm(0, 0.0001)
beta1 ~ dnorm(0, 0.0001)
tau2 ~ dlnorm(-2.49, 1/1.52^2)
betasum <- beta0 + beta1
}"
inits <- list(theta=rep(0,k), beta0=0, beta1=1, tau2=1, .RNG.name="base::Mersenne-Twister")
inits <- list(inits, inits, inits)
inits[[1]]$.RNG.seed <- 12341
inits[[2]]$.RNG.seed <- 12342
inits[[3]]$.RNG.seed <- 12343
model <- jags.model(textConnection(me.model), inits=inits, data=jags.data, n.chains=3, quiet=TRUE)
update(model, n.iter=10000, progress.bar="none")
bma.me <- coda.samples(model, variable.names=c("beta0","beta1","tau2","betasum"),
n.iter=100000, progress.bar="none")
bma.me <- summary(bma.me)$quantiles[,c(3,1,5)]
bma.me <- rbind(bma.me[c(1,2,4),], exp.beta0=exp(bma.me[1,]), exp.betasum=exp(bma.me[3,]))
round(bma.me, digits=2)
## 50% 2.5% 97.5%
## beta0 -1.03 -1.42 -0.74
## beta1 -0.60 -1.16 0.02
## tau2 0.07 0.00 0.44
## exp.beta0 0.36 0.24 0.48
## exp.betasum 0.20 0.12 0.32
# Table 14.4: Recurrence of Violence Data
dat <- read.table(header=TRUE, text = "
study year xt nt xc nc
'Bronx' 2005 20 202 11 218
'Brooklyn' 2000 13 129 100 386
'Broward' 2000 52 216 45 188
'San Diego Navy' 2000 63 218 75 214")
dat <- escalc(measure="OR", ai=xt, n1i=nt, ci=xc, n2i=nc, data=dat, digits=2)
dat
## study year xt nt xc nc yi vi
## 1 Bronx 2005 20 202 11 218 0.73 0.15
## 2 Brooklyn 2000 13 129 100 386 -1.14 0.10
## 3 Broward 2000 52 216 45 188 0.01 0.05
## 4 San Diego Navy 2000 63 218 75 214 -0.28 0.04
Using the same code as above, one can repeat all analyses with this
dataset. Doing so yields the following results (code now shown; see
cooper2019.rmd
file for the code):
## or ci.lb ci.ub tau2 ci.lb ci.ub
## Frequentist RE model 0.82 0.45 1.51 0.31 0.07 8.15
## Bayesian RE model, uniform(0,2) prior for tau 0.82 0.28 2.54 0.74 0.06 3.51
## Bayesian RE model, half-normal(0,0.5^2) prior for tau 0.82 0.41 1.67 0.30 0.02 1.28
## Bayesian RE model, gamma(0.001,0.001) prior for 1/tau^2 0.82 0.29 2.36 0.42 0.01 6.05
## Bayesian RE model, log-normal(-2.01,1.64^2) prior for tau^2 0.82 0.39 1.74 0.29 0.03 1.93
# note: the prior for the last model was log-normal(-2.01,1.64^2) (see page 311)
# and not log-normal(-3.95,1.79^2) as stated in the table
# Table 17.1: Oral Anticoagulant Therapy Data
dat <- dat.anand1999
# reorder rows
dat <- dat[c(1:21,29,22:28,30:34),]
rownames(dat) <- 1:nrow(dat)
# add 'age' variable
dat$age <- 1999 - dat$year
# compute log odds ratios and corresponding sampling variances
dat <- escalc(measure="OR", ai=ai, n1i=n1i, ci=ci, n2i=n2i, data=dat, to="all", digits=2)
dat$yi <- round(dat$yi, 2)
dat$vi <- round(dat$vi, 2)
# flip signs of some of the log odds ratios (see note below)
sel <- c(4:6,11:14,16,19,20,22,24:27,29,30,32,33)
dat$yi[sel] <- -1 * dat$yi[sel]
# code intensity for study 22 as 'moderate' (see note below)
dat$intensity[22] <- "moderate"
# turn intensity into a factor with the desired order of levels
dat$intensity <- factor(dat$intensity, levels=c("low", "high", "moderate"))
# add mcar and mar dummy variables
dat$mcar <- c(0,1,0,1,1,1,1,0,1,0,0,1,0,1,1,1,1,0,1,0,1,0,1,1,1,1,1,1,1,1,1,0,1,1)
dat$mar <- c(0,0,1,0,0,1,0,1,0,0,0,0,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1)
# select (order of) columns
dat <- dat[c(11,12,3,2,10,13,14)]
dat
## yi vi intensity year age mcar mar
## 1 3.02 2.21 high 1960 39 0 0
## 2 -1.84 0.80 high 1960 39 1 0
## 3 0.24 0.16 high 1961 38 0 1
## 4 0.15 0.07 high 1961 38 1 0
## 5 0.47 0.07 high 1964 35 1 0
## 6 0.38 0.28 high 1964 35 1 1
## 7 -0.38 0.18 high 1966 33 1 0
## 8 -0.47 0.22 high 1967 32 0 1
## 9 -0.25 0.07 high 1967 32 1 0
## 10 0.22 0.10 high 1969 30 0 0
## 11 0.84 0.08 high 1969 30 0 0
## 12 0.35 0.04 high 1980 19 1 0
## 13 0.33 0.02 high 1990 9 0 1
## 14 0.12 0.01 high 1994 5 1 1
## 15 -1.81 0.82 high 1969 30 1 1
## 16 0.43 0.02 high 1974 25 1 1
## 17 0.18 0.06 high 1980 19 1 1
## 18 0.39 0.07 high 1980 19 0 1
## 19 0.90 0.41 high 1993 6 1 1
## 20 0.65 0.29 high 1996 3 0 0
## 21 1.42 2.74 high 1990 9 1 1
## 22 0.14 4.06 moderate 1990 9 0 1
## 23 0.04 1.36 moderate 1982 17 1 1
## 24 0.35 0.35 moderate 1981 18 1 1
## 25 0.08 0.03 moderate 1982 17 1 1
## 26 0.06 0.03 moderate 1969 30 1 1
## 27 0.43 0.07 moderate 1964 35 1 1
## 28 -1.16 0.93 moderate 1986 13 1 1
## 29 0.75 0.98 moderate 1982 17 1 1
## 30 0.81 0.60 moderate 1998 1 1 1
## 31 0.04 0.82 moderate 1994 5 1 1
## 32 0.35 0.70 low 1998 1 0 1
## 33 0.34 0.06 low 1997 2 1 1
## 34 0.15 0.02 low 1997 2 1 1
# note: compared to the original dataset, some of the log odds ratios have
# flipped signs; also, intensity is coded as 'moderate' for study 22, but
# this was actually 'high' in the original dataset; all of these differences
# are introduced above
# fit model to full data
res.full <- rma(yi, vi, mods = ~ intensity + age, data=dat)
# fit model to MCAR data
dat$age[dat$mcar == 0] <- NA
res.mcar <- rma(yi, vi, mods = ~ intensity + age, data=dat)
# fit model to MAR data
dat$age <- 1999 - dat$year
dat$age[dat$mar == 0] <- NA
res.mar <- rma(yi, vi, mods = ~ intensity + age, data=dat)
# Table 17.2
tab <- cbind(full = coef(summary(res.full))[1:2],
mcar = coef(summary(res.mcar))[1:2],
mar = coef(summary(res.mar))[1:2])
round(tab, digits=3)
## full.estimate full.se mcar.estimate mcar.se mar.estimate mar.se
## intrcpt 0.213 0.142 0.214 0.153 0.193 0.122
## intensityhigh 0.060 0.180 0.014 0.209 -0.008 0.152
## intensitymoderate -0.033 0.216 -0.011 0.239 -0.149 0.206
## age -0.001 0.005 -0.002 0.006 0.004 0.006
# fit model to MCAR data with mean imputation
dat$age <- 1999 - dat$year
dat$age[dat$mcar == 0] <- NA
dat$age <- replmiss(dat$age, mean(dat$age, na.rm=TRUE))
res.mcar <- rma(yi, vi, mods = ~ intensity + age, data=dat)
# fit model to MAR data with mean imputation
dat$age <- 1999 - dat$year
dat$age[dat$mar == 0] <- NA
dat$age <- replmiss(dat$age, mean(dat$age, na.rm=TRUE))
res.mar <- rma(yi, vi, mods = ~ intensity + age, data=dat)
# Table 17.3
tab <- cbind(mcar = coef(summary(res.mcar))[1:2],
mar = coef(summary(res.mar))[1:2])
round(tab, digits=3)
## mcar.estimate mcar.se mar.estimate mar.se
## intrcpt 0.215 0.144 0.197 0.128
## intensityhigh 0.063 0.189 0.001 0.160
## intensitymoderate -0.029 0.226 -0.132 0.216
## age -0.001 0.006 0.003 0.006
# fit model to MCAR data with regression imputation
dat$age <- 1999 - dat$year
dat$age[dat$mcar == 0] <- NA
fit <- lm(age ~ intensity + yi, data=dat)
dat$age <- replmiss(dat$age, predict(fit, newdata=dat))
res.mcar <- rma(yi, vi, mods = ~ intensity + age, data=dat)
# fit model to MAR data with regression imputation
dat$age <- 1999 - dat$year
dat$age[dat$mar == 0] <- NA
fit <- lm(age ~ intensity + yi, data=dat)
dat$age <- replmiss(dat$age, predict(fit, newdata=dat))
res.mar <- rma(yi, vi, mods = ~ intensity + age, data=dat)
# Table 17.4
tab <- cbind(mcar = coef(summary(res.mcar))[1:2],
mar = coef(summary(res.mar))[1:2])
round(tab, digits=3)
## mcar.estimate mcar.se mar.estimate mar.se
## intrcpt 0.221 0.151 0.220 0.149
## intensityhigh 0.098 0.204 0.086 0.192
## intensitymoderate 0.012 0.238 0.013 0.243
## age -0.003 0.006 -0.003 0.007
# install the 'mice' package (if it is not already installed) and load it
if (!require(mice)) {
install.packages("mice")
library(mice)
}
# multiple imputation using the MAR data
dat$age <- 1999 - dat$year
dat$age[dat$mar == 0] <- NA
predMatrix <- make.predictorMatrix(dat)
predMatrix[,] <- 0
predMatrix["age", c("yi", "intensity")] <- 1
predMatrix
## yi vi intensity year age mcar mar
## yi 0 0 0 0 0 0 0
## vi 0 0 0 0 0 0 0
## intensity 0 0 0 0 0 0 0
## year 0 0 0 0 0 0 0
## age 1 0 1 0 0 0 0
## mcar 0 0 0 0 0 0 0
## mar 0 0 0 0 0 0 0
## yi vi intensity year age mcar mar
## "" "" "" "" "norm" "" ""
imp <- mice(dat, print=FALSE, m=100, predictorMatrix=predMatrix, method=impMethod, seed=1234)
fit <- with(imp, rma(yi, vi, mods = ~ intensity + age))
pool <- summary(pool(fit))
# Table 17.5
pool[-1] <- round(pool[-1], digits=3)
pool[1:3]
## term estimate std.error
## 1 intercept 0.214 0.144
## 2 intensityhigh 0.068 0.186
## 3 intensitymoderate -0.018 0.243
## 4 age -0.002 0.007
# note: results in the table are different; this could be due to several reasons:
#
# 1) multiple imputation involves an element of randomness (above the seed of the
# random number generator is set to some arbitrary value for reproducibility,
# but a different seed would generate (slightly) different results)
# 2) the chapter author used the Amelia package instead of the mice package for
# the analyses and there could be differences in how these packages implement
# the multiple imputation algorithms
#
# the coefficients are actually not all that different, but the SEs shown in the
# chapter are much larger than the ones obtained above
#
# note: above m=100 imputations were used (instead of 10 as was done in the chapter)
# to obtain more stable results
#
# let's try this with the Amelia package
# install the 'Amelia' package (if it is not already installed) and load it
if (!require(Amelia)) {
install.packages("Amelia")
library(Amelia)
}
# multiple imputation using the MAR data
set.seed(1234)
invisible(capture.output(imp <- amelia(dat, m=100, idvars=c(2,4,6,7), noms=3, p2s=0)))
# the invisible(capture.output()) part is used here to suppress some output that
# is generated by the amelia() function, but in general this can be left out
fit <- lapply(imp$imputations,
function(x) if (length(x) == 1L) NULL else rma(yi, vi, mods = ~ intensity + age, data=x))
fit <- fit[!sapply(fit, is.null)]
b <- sapply(fit, function(x) coef(x))
se <- sapply(fit, function(x) x$se)
pool2 <- mi.meld(b, se, byrow=FALSE)
pool2 <- data.frame(estimate=pool2$q.mi[1,], se=pool2$se.mi[1,])
round(pool2, digits=3)
## estimate se
## intrcpt 0.214 0.144
## intensityhigh 0.068 0.190
## intensitymoderate -0.017 0.248
## age -0.002 0.007
# note: these results are very similar to the ones obtained above using the mice
# package and the SEs again differ substantially from those shown in Table 17.5
Part of the code for this chapter (by Jack Vevea, Kathleen Coburn, and Alexander Sutton) is adapted from the chapter itself (see section 18.6).
Below is the code for the analysis of the irritable bowl syndrome (IBS) data.
# compute log risk ratios and corresponding sampling variances
dat <- escalc(measure="RR", ai=x.a, n1i=n.a, ci=x.p, n2i=n.p,
data=dat.dorn2007, slab=paste(study, ", ", year))
# fit random-effects model using ML estimation
res <- rma(yi, vi, data=dat, digits=2, method="ML")
res
## Random-Effects Model (k = 19; tau^2 estimator: ML)
##
## tau^2 (estimated amount of total heterogeneity): 0.15 (SE = 0.07)
## tau (square root of estimated tau^2 value): 0.39
## I^2 (total heterogeneity / total variability): 72.90%
## H^2 (total variability / sampling variability): 3.69
##
## Test for Heterogeneity:
## Q(df = 18) = 62.40, p-val < .01
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.42 0.11 3.77 <.01 0.20 0.63 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# funnel plot
par(mar=c(5,4,2,2))
funnel(res, atransf=exp, at=log(c(0.25, 0.5, 1, 2, 4, 8)),
ylim=c(0,0.6), steps=7, las=1)
# cumulative meta-analysis (in order of precision) and corresponding forest plot
sav <- cumul(res, order=vi)
par(mar=c(5,4,2,2))
forest(sav, header=TRUE, xlim=c(-1,1.5))
An alternative way of plotting the results from such a cumulative meta-analysis is to plot the pooled estimates (in this case, the estimated average log risk ratios) against the estimates of the amount of heterogeneity as studies are added. The color gradient of the points/lines indicates the order of the cumulative results (light gray at the beginning, dark gray at the end).
# plot cumulative meta-analysis results
par(mar=c(5,4,2,2))
plot(sav, xlim=c(0.1,0.5), ylim=c(0,0.16), las=1)
# trim and fill analysis
tf.l.l0 <- trimfill(res, side="left", estimator="L0")
tf.l.r0 <- trimfill(res, side="left", estimator="R0")
tf.r.l0 <- trimfill(res, side="right", estimator="L0")
tf.r.r0 <- trimfill(res, side="right", estimator="R0")
# show results for side="left" and estimator="R0"
tf.l.r0
## Estimated number of missing studies on the left side: 1 (SE = 2.00)
## Test of H0: no missing studies on the left side: p-val = 0.25
##
## Random-Effects Model (k = 20; tau^2 estimator: ML)
##
## tau^2 (estimated amount of total heterogeneity): 0.18 (SE = 0.08)
## tau (square root of estimated tau^2 value): 0.43
## I^2 (total heterogeneity / total variability): 75.56%
## H^2 (total variability / sampling variability): 4.09
##
## Test for Heterogeneity:
## Q(df = 19) = 69.35, p-val < .01
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.38 0.12 3.26 <.01 0.15 0.61 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Egger regression using a weighted regression model with the standard error as predictor
regtest(res, model="lm")
## Regression Test for Funnel Plot Asymmetry
##
## Model: weighted regression with multiplicative dispersion
## Predictor: standard error
##
## Test for Funnel Plot Asymmetry: t = 2.22, df = 17, p = 0.04
## Limit Estimate (as sei -> 0): b = -0.23 (CI: -0.78, 0.32)
## Regression Test for Funnel Plot Asymmetry
##
## Model: mixed-effects meta-regression model
## Predictor: standard error
##
## Test for Funnel Plot Asymmetry: z = 2.64, p < .01
## Limit Estimate (as sei -> 0): b = -0.23 (CI: -0.72, 0.27)
# The PET-PEESE procedure is this: Start with Egger regression using a weighted regression
# model with the standard error as predictor. If the limit estimate (i.e., intercept) is not
# significant, the adjusted estimate is the limit estimate. On the other hand, if the limit
# estimate is significant, repeat the Egger regression but using the sampling variances as
# predictor. The adjusted estimate is then the limit estimate from that model.
regtest(res, model="lm", ret.fit=TRUE)
## Regression Test for Funnel Plot Asymmetry
##
## Model: weighted regression with multiplicative dispersion
## Predictor: standard error
##
## Call:
## lm(formula = yi ~ X - 1, weights = 1/vi)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -2.5967 -1.6355 0.6988 1.4555 1.9691
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## Xintrcpt -0.2313 0.2612 -0.885 0.3883
## Xsei 2.4617 1.1087 2.220 0.0403 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.687 on 17 degrees of freedom
## Multiple R-squared: 0.4949, Adjusted R-squared: 0.4355
## F-statistic: 8.329 on 2 and 17 DF, p-value: 0.00301
##
## Test for Funnel Plot Asymmetry: t = 2.22, df = 17, p = 0.04
## Limit Estimate (as sei -> 0): b = -0.23 (CI: -0.78, 0.32)
# note: since the intercept is not significant (p = .39), the adjusted estimate is the limit
# estimate; but just for illustration purposes, if the intercept had been significant, then
# the adjusted estimate would be the limit estimate from the following model
regtest(res, model="lm", predictor="vi")
## Regression Test for Funnel Plot Asymmetry
##
## Model: weighted regression with multiplicative dispersion
## Predictor: sampling variance
##
## Test for Funnel Plot Asymmetry: t = 2.13, df = 17, p = 0.05
## Limit Estimate (as vi -> 0): b = 0.09 (CI: -0.20, 0.38)
## Rank Correlation Test for Funnel Plot Asymmetry
##
## Kendall's tau = 0.26, p = 0.12
For a p-curve analysis, one needs to use the online app at http://p-curve.com/.
# install the 'puniform' package (if it is not already installed) and load it
if (!require(puniform)) {
install.packages("puniform")
library(puniform)
}
##
## Method: P
##
## Effect size estimation p-uniform
##
## est ci.lb ci.ub L.0 pval ksig
## 0.6439 0.4067 0.9005 -4.03 <.001 10
##
## ===
##
## Publication bias test p-uniform
##
## L.pb pval
## -2.6085 0.9955
##
## ===
##
## Fixed-effect meta-analysis
##
## est.fe se.fe zval.fe pval.fe ci.lb.fe ci.ub.fe Qstat Qpval
## 0.3123 0.0541 5.777 <.001 0.2063 0.4182 62.4013 <.001
## Test of Excess Significance
##
## Observed Number of Significant Findings: 10 (out of 19)
## Expected Number of Significant Findings: 4.9233
## Observed Number / Expected Number: 2.0311
##
## Estimated Power of Tests (based on theta = 0.3123)
##
## min q1 median q3 max
## 0.0894 0.1270 0.2118 0.3537 0.5397
##
## Test of Excess Significance: p = 0.0078 (exact test)
## Limit Estimate (theta_lim): 0.3971 (where p = 0.1)
# install the 'selectMeta' package (if it is not already installed) and load it
if (!require(selectMeta)) {
install.packages("selectMeta")
library(selectMeta)
}
# Dear and Begg weight-function model
sav <- DearBegg(dat$yi, sqrt(dat$vi), trace=FALSE)
par(mar=c(5,4,2,2))
plot(0, 0, type="n", xlim=c(0, 1), ylim=c(0, 1), xlab="p-Values", ylab="Estimated Weight Function")
ps <- seq(0, 1, by = 0.01)
pval <- summary(dat)$pval
rug(pval, lwd=3)
weightLine(pval, w=sav$w, col0="black", lwd0=1, lty0="solid")
# Rufibach weight-function model
sav <- DearBeggMonotone(dat$yi, sqrt(dat$vi), trace=FALSE)
par(mar=c(5,4,2,2))
plot(0, 0, type="n", xlim=c(0, 1), ylim=c(0, 1), xlab="p-Values", ylab="Estimated Weight Function")
ps <- seq(0, 1, by = 0.01)
pval <- summary(dat)$pval
rug(pval, lwd=3)
weightLine(pval, w=sav$w, col0="black", lwd0=1, lty0="solid")
# Vevea and Hedges (1995) model (with a single cutpoint at p=0.025)
selmodel(res, type="stepfun", steps=0.025)
## Random-Effects Model (k = 19; tau^2 estimator: ML)
##
## tau^2 (estimated amount of total heterogeneity): 0.09 (SE = 0.07)
## tau (square root of estimated tau^2 value): 0.31
##
## Test for Heterogeneity:
## LRT(df = 1) = 6.10, p-val = 0.01
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.16 0.15 1.01 0.31 -0.14 0.45
##
## Test for Selection Model Parameters:
## LRT(df = 1) = 3.59, p-val = 0.06
##
## Selection Model Results:
##
## k estimate se zval pval ci.lb ci.ub
## 0 < p <= 0.025 10 1.00 --- --- --- --- ---
## 0.025 < p <= 1 9 0.19 0.18 -4.57 <.01 0.00 0.54 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Table 18.1: Sample Selection Patterns for the Vevea and Woods Method
ssp <- data.frame(
steps = c(0.005, 0.01, 0.05, 0.10, 0.25, 0.35, 0.50, 0.65, 0.75, 0.90, 0.95, 0.99, 0.995, 1),
weights.mod.1 = c(1, 0.99, 0.95, 0.80, 0.75, 0.65, 0.60, 0.55, 0.50, 0.50, 0.50, 0.50, 0.50, 0.50),
weights.sev.1 = c(1, 0.99, 0.90, 0.75, 0.60, 0.50, 0.40, 0.35, 0.30, 0.25, 0.10, 0.10, 0.10, 0.10),
weights.mod.2 = c(1, 0.99, 0.95, 0.90, 0.80, 0.75, 0.60, 0.60, 0.75, 0.80, 0.90, 0.95, 0.99, 1.00),
weights.sev.2 = c(1, 0.99, 0.90, 0.75, 0.60, 0.50, 0.25, 0.25, 0.50, 0.60, 0.75, 0.90, 0.99, 1.00))
ssp
## steps weights.mod.1 weights.sev.1 weights.mod.2 weights.sev.2
## 1 0.005 1.00 1.00 1.00 1.00
## 2 0.010 0.99 0.99 0.99 0.99
## 3 0.050 0.95 0.90 0.95 0.90
## 4 0.100 0.80 0.75 0.90 0.75
## 5 0.250 0.75 0.60 0.80 0.60
## 6 0.350 0.65 0.50 0.75 0.50
## 7 0.500 0.60 0.40 0.60 0.25
## 8 0.650 0.55 0.35 0.60 0.25
## 9 0.750 0.50 0.30 0.75 0.50
## 10 0.900 0.50 0.25 0.80 0.60
## 11 0.950 0.50 0.10 0.90 0.75
## 12 0.990 0.50 0.10 0.95 0.90
## 13 0.995 0.50 0.10 0.99 0.99
## 14 1.000 0.50 0.10 1.00 1.00
# Vevea and Wood (2005) model using moderate one-tailed selection
selmodel(res, type="stepfun", steps=ssp$steps, delta=ssp$weights.mod.1)
## Random-Effects Model (k = 19; tau^2 estimator: ML)
##
## tau^2 (estimated amount of total heterogeneity): 0.16 (SE = 0.09)
## tau (square root of estimated tau^2 value): 0.40
##
## Test for Heterogeneity:
## LRT(df = 1) = 17.74, p-val < .01
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.31 0.12 2.61 <.01 0.08 0.54 **
##
## Test for Selection Model Parameters:
## LRT(df = 0) = 0.00, p-val = NA
##
## Selection Model Results:
##
## k estimate se zval pval ci.lb ci.ub
## 0 < p <= 0.005 7 1.00 --- --- --- --- ---
## 0.005 < p <= 0.01 2 0.99 --- --- --- --- ---
## 0.01 < p <= 0.05 2 0.95 --- --- --- --- ---
## 0.05 < p <= 0.1 0 0.80 --- --- --- --- ---
## 0.1 < p <= 0.25 2 0.75 --- --- --- --- ---
## 0.25 < p <= 0.35 0 0.65 --- --- --- --- ---
## 0.35 < p <= 0.5 0 0.60 --- --- --- --- ---
## 0.5 < p <= 0.65 2 0.55 --- --- --- --- ---
## 0.65 < p <= 0.75 1 0.50 --- --- --- --- ---
## 0.75 < p <= 0.9 3 0.50 --- --- --- --- ---
## 0.9 < p <= 0.95 0 0.50 --- --- --- --- ---
## 0.95 < p <= 0.99 0 0.50 --- --- --- --- ---
## 0.99 < p <= 0.995 0 0.50 --- --- --- --- ---
## 0.995 < p <= 1 0 0.50 --- --- --- --- ---
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Vevea and Wood (2005) model using severe one-tailed selection
selmodel(res, type="stepfun", steps=ssp$steps, delta=ssp$weights.sev.1)
## Random-Effects Model (k = 19; tau^2 estimator: ML)
##
## tau^2 (estimated amount of total heterogeneity): 0.23 (SE = 0.14)
## tau (square root of estimated tau^2 value): 0.48
##
## Test for Heterogeneity:
## LRT(df = 1) = 19.45, p-val < .01
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.11 0.18 0.57 0.57 -0.26 0.47
##
## Test for Selection Model Parameters:
## LRT(df = 0) = 0.66, p-val = NA
##
## Selection Model Results:
##
## k estimate se zval pval ci.lb ci.ub
## 0 < p <= 0.005 7 1.00 --- --- --- --- ---
## 0.005 < p <= 0.01 2 0.99 --- --- --- --- ---
## 0.01 < p <= 0.05 2 0.90 --- --- --- --- ---
## 0.05 < p <= 0.1 0 0.75 --- --- --- --- ---
## 0.1 < p <= 0.25 2 0.60 --- --- --- --- ---
## 0.25 < p <= 0.35 0 0.50 --- --- --- --- ---
## 0.35 < p <= 0.5 0 0.40 --- --- --- --- ---
## 0.5 < p <= 0.65 2 0.35 --- --- --- --- ---
## 0.65 < p <= 0.75 1 0.30 --- --- --- --- ---
## 0.75 < p <= 0.9 3 0.25 --- --- --- --- ---
## 0.9 < p <= 0.95 0 0.10 --- --- --- --- ---
## 0.95 < p <= 0.99 0 0.10 --- --- --- --- ---
## 0.99 < p <= 0.995 0 0.10 --- --- --- --- ---
## 0.995 < p <= 1 0 0.10 --- --- --- --- ---
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Vevea and Wood (2005) model using moderate two-tailed selection
selmodel(res, type="stepfun", steps=ssp$steps, delta=ssp$weights.mod.2)
## Random-Effects Model (k = 19; tau^2 estimator: ML)
##
## tau^2 (estimated amount of total heterogeneity): 0.14 (SE = 0.07)
## tau (square root of estimated tau^2 value): 0.37
##
## Test for Heterogeneity:
## LRT(df = 1) = 16.73, p-val < .01
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.37 0.11 3.52 <.01 0.16 0.58 ***
##
## Test for Selection Model Parameters:
## LRT(df = 0) = 0.24, p-val = NA
##
## Selection Model Results:
##
## k estimate se zval pval ci.lb ci.ub
## 0 < p <= 0.005 7 1.00 --- --- --- --- ---
## 0.005 < p <= 0.01 2 0.99 --- --- --- --- ---
## 0.01 < p <= 0.05 2 0.95 --- --- --- --- ---
## 0.05 < p <= 0.1 0 0.90 --- --- --- --- ---
## 0.1 < p <= 0.25 2 0.80 --- --- --- --- ---
## 0.25 < p <= 0.35 0 0.75 --- --- --- --- ---
## 0.35 < p <= 0.5 0 0.60 --- --- --- --- ---
## 0.5 < p <= 0.65 2 0.60 --- --- --- --- ---
## 0.65 < p <= 0.75 1 0.75 --- --- --- --- ---
## 0.75 < p <= 0.9 3 0.80 --- --- --- --- ---
## 0.9 < p <= 0.95 0 0.90 --- --- --- --- ---
## 0.95 < p <= 0.99 0 0.95 --- --- --- --- ---
## 0.99 < p <= 0.995 0 0.99 --- --- --- --- ---
## 0.995 < p <= 1 0 1.00 --- --- --- --- ---
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Vevea and Wood (2005) model using severe two-tailed selection
selmodel(res, type="stepfun", steps=ssp$steps, delta=ssp$weights.sev.2)
## Random-Effects Model (k = 19; tau^2 estimator: ML)
##
## tau^2 (estimated amount of total heterogeneity): 0.11 (SE = 0.06)
## tau (square root of estimated tau^2 value): 0.34
##
## Test for Heterogeneity:
## LRT(df = 1) = 14.29, p-val < .01
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.32 0.10 3.30 <.01 0.13 0.51 ***
##
## Test for Selection Model Parameters:
## LRT(df = 0) = 0.00, p-val = NA
##
## Selection Model Results:
##
## k estimate se zval pval ci.lb ci.ub
## 0 < p <= 0.005 7 1.00 --- --- --- --- ---
## 0.005 < p <= 0.01 2 0.99 --- --- --- --- ---
## 0.01 < p <= 0.05 2 0.90 --- --- --- --- ---
## 0.05 < p <= 0.1 0 0.75 --- --- --- --- ---
## 0.1 < p <= 0.25 2 0.60 --- --- --- --- ---
## 0.25 < p <= 0.35 0 0.50 --- --- --- --- ---
## 0.35 < p <= 0.5 0 0.25 --- --- --- --- ---
## 0.5 < p <= 0.65 2 0.25 --- --- --- --- ---
## 0.65 < p <= 0.75 1 0.50 --- --- --- --- ---
## 0.75 < p <= 0.9 3 0.60 --- --- --- --- ---
## 0.9 < p <= 0.95 0 0.75 --- --- --- --- ---
## 0.95 < p <= 0.99 0 0.90 --- --- --- --- ---
## 0.99 < p <= 0.995 0 0.99 --- --- --- --- ---
## 0.995 < p <= 1 0 1.00 --- --- --- --- ---
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# install the 'metasens' package (if it is not already installed) and load it
if (!require(metasens)) {
install.packages("metasens")
library(metasens)
}
# Copas and Shi (2001) model
res <- metagen(yi, sqrt(vi), method.tau="ML", data=dat)
sav <- copas(res)
plot(sav)
## Copas selection model analysis
##
## p.publ 95%-CI tau^2 tau p.trt p.rsb N
## 1.0000 0.4176 [ 0.1966; 0.6386] 0.1534 0.3917 0.0002 0.0144 0
## 0.9612 0.3999 [ 0.1903; 0.6095] 0.1391 0.3729 0.0002 0.0225 0
## 0.8468 0.3501 [ 0.1535; 0.5466] 0.1197 0.3460 0.0005 0.0541 2
## 0.7417 0.3004 [ 0.1295; 0.4714] 0.1117 0.3342 0.0006 0.0997 3
## 0.6475 0.2506 [-0.1802; 0.6813] 0.1097 0.3313 0.2542 0.1504 6
## 0.5432 0.1925 [-0.0460; 0.4310] 0.1115 0.3339 0.1137 0.1957 9
## 0.4352 0.1511 [-0.1019; 0.4041] 0.1083 0.3291 0.2418 0.2460 15
## 0.3288 0.1026 [-0.1695; 0.3746] 0.1047 0.3236 0.4599 0.3191 24
##
## Adjusted estimate 0.2506 [-0.1802; 0.6813] 0.1097 0.3313 0.2542 0.1504 6
## Unadjusted estimate 0.4176 [ 0.2003; 0.6348] 0.1534 0.3917 0.0002
##
## Significance level for test of residual selection bias: 0.1
##
## min max
## range of gamma0: -0.7558 2.0000
## range of gamma1: 0.0000 0.2723
##
## Largest standard error (SE): 0.5383
##
## Range of probability publishing trial with largest SE:
## min max
## 0.2249 0.9939
##
## Calculation of orthogonal line:
##
## level nobs adj.r.square slope se.slope
## -0.05 5 0.9999830 -1.362614 0.002809321
## 0.00 9 0.9998969 -1.440954 0.005173595
## 0.05 14 0.9994926 -1.513392 0.009457356
## 0.10 19 0.9925764 -1.774954 0.036173183
## 0.15 24 0.9886879 -2.100479 0.046836982
## 0.20 28 0.9914366 -2.508877 0.044866233
## 0.25 26 0.9997682 -2.980742 0.009078334
## 0.30 26 0.9995475 -3.081597 0.013113046
## 0.35 26 0.9993143 -3.196848 0.016748454
## 0.40 24 0.9990672 -3.426336 0.021830322
##
## Legend:
## p.publ - Probability of publishing study with largest SE
## p.trt - P-value for test of overall treatment effect
## p.rsb - P-value for test of residual selection bias
## N - Estimated number of unpublished studies
## Result of limit meta-analysis:
##
## Random effects model 95%-CI z pval
## Adjusted estimate 0.0854 [-0.2800; 0.4508] 0.46 0.6469
## Unadjusted estimate 0.4176 [ 0.2003; 0.6348] 3.77 0.0002
##
## Quantifying heterogeneity:
## tau^2 = 0.1534; I^2 = 71.2% [54.1%; 81.9%]; G^2 = 99.9%
##
## Test of heterogeneity:
## Q d.f. p-value
## 62.40 18 < 0.0001
##
## Test of small-study effects:
## Q-Q' d.f. p-value
## 14.03 1 0.0002
##
## Test of residual heterogeneity beyond small-study effects:
## Q' d.f. p-value
## 48.37 17 < 0.0001
##
## Details on adjustment method:
## - expectation (beta0)
## [1] 0.5537899
## [1] 0.6914625
## [1] 0.8413447
## [1] 0.583998
## [1] 0.5741702
# note: this is the proportion in the first group scoring above the median; the proportion
# in the second group scoring above the median is simply 1 minus the computed value:
1 - transf.dtobesd(0.3)
## [1] 0.4258298
# converting a 2x2 table to a correlation (phi coefficient)
escalc(measure="PHI", ai=100, bi=50, ci=60, di=90)
## yi vi
## 1 0.2673 0.0031
# converting a 2x2 table to an odds ratio
summary(escalc(measure="OR", ai=100, bi=50, ci=60, di=90), transf=exp)
## yi ci.lb ci.ub
## 1 3.0000 1.8729 4.8053
# converting a 2x2 table to a risk ratio
summary(escalc(measure="RR", ai=100, bi=50, ci=60, di=90), transf=exp)
## yi ci.lb ci.ub
## 1 1.6667 1.3291 2.0900
## yi vi
## 1 0.2667 0.0031
# converting a 2x2 table to a number needed to treat
summary(escalc(measure="RD", ai=100, bi=50, ci=60, di=90), transf=function(x) 1/x)
## yi ci.lb ci.ub
## 1 3.7500 2.6634 6.3344
# using effect-size translations with confidence intervals (after a meta-analysis)
# note: using some made-up data that yield the same results as shown on page 443
# example with standardized mean differences
dat <- data.frame(di = c(0.03, 0.93, 0.66),
vi = c(0.072, 0.063, 0.047))
res <- rma(di, vi, data=dat, digits=2)
res
## Random-Effects Model (k = 3; tau^2 estimator: REML)
##
## tau^2 (estimated amount of total heterogeneity): 0.14 (SE = 0.20)
## tau (square root of estimated tau^2 value): 0.37
## I^2 (total heterogeneity / total variability): 69.70%
## H^2 (total variability / sampling variability): 3.30
##
## Test for Heterogeneity:
## Q(df = 2) = 6.28, p-val = 0.04
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.55 0.26 2.15 0.03 0.05 1.05 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## pred ci.lb ci.ub pi.lb pi.ub
## 0.71 0.52 0.85 0.37 0.92
# example with (log) odds ratios
dat <- data.frame(logori = c(0.40, 0.74, 0.99, 0.84, 0.75),
vi = c(0.111, 0.093, 0.165, 0.128, 0.114))
res <- rma(logori, vi, data=dat, digits=3)
res
## Random-Effects Model (k = 5; tau^2 estimator: REML)
##
## tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.083)
## 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 = 4) = 1.488, p-val = 0.829
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.724 0.154 4.715 <.001 0.423 1.025 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## pred ci.lb ci.ub pi.lb pi.ub
## 2.06 1.53 2.79 1.53 2.79
## pred ci.lb ci.ub pi.lb pi.ub
## 0.10 0.05 0.16 0.05 0.16
See also the R functions provided by the chapter authors (Jeffrey Valentine, Ariel Aloe, and Sandra Wilson) at the end of the chapter. The code is also available on the publisher’s website for the book (see here).
## yi vi
## 1 10 562.5
## 2 20 562.5
## 3 30 562.5
## 4 40 562.5
## 5 50 562.5
## 6 60 562.5
## 7 70 562.5
## 8 80 562.5
## 9 90 562.5
## Random-Effects Model (k = 9; tau^2 estimator: DL)
##
## tau^2 (estimated amount of total heterogeneity): 187.5000 (SE = 375.0000)
## tau (square root of estimated tau^2 value): 13.6931
## I^2 (total heterogeneity / total variability): 25.00%
## H^2 (total variability / sampling variability): 1.33
##
## Test for Heterogeneity:
## Q(df = 8) = 10.6667, p-val = 0.2213
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 50.0000 9.1287 5.4772 <.0001 32.1081 67.8919 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# approximate 95% prediction interval (based on equation 20.6)
predict(res, pi.type="simple", digits=0)
## pred se ci.lb ci.ub pi.lb pi.ub
## 50 9 32 68 23 77
# approximate 95% prediction interval (based on equation 20.7)
predict(res, pi.type="riley", digits=0)
## pred se ci.lb ci.ub pi.lb pi.ub
## 50 9 32 68 11 89
Note: 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"
, equation 20.7
is used.
# data for Figure 20.3
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
# fit equal- and random-effects models
res.ee <- rma(yi, vi, data=dat, method="EE")
res.re <- rma(yi, vi, data=dat, method="DL")
# show random-effects model results
print(res.re, digits=3)
## Random-Effects Model (k = 6; tau^2 estimator: DL)
##
## tau^2 (estimated amount of total heterogeneity): 0.037 (SE = 0.042)
## tau (square root of estimated tau^2 value): 0.193
## I^2 (total heterogeneity / total variability): 58.34%
## H^2 (total variability / sampling variability): 2.40
##
## Test for Heterogeneity:
## Q(df = 5) = 12.003, p-val = 0.035
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.358 0.105 3.404 <.001 0.152 0.565 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Figure 20.3: Meta-Analysis Showing Relative Weights for the Equal- and Random-Effects Models
wi <- fmtx(cbind(weights(res.ee), weights(res.re)), digits=2)
par(mar=c(5,4,2,2))
forest(res.ee, header=TRUE, xlim=c(-1.5,4), ylim=c(-2.5,9), psize=1, efac=c(0,1),
ilab=wi, ilab.xpos=c(1.7,2.4), ilab.pos=2, digits=c(3,2))
addpoly(res.re)
text(1.45, 8, "Weight\n(Equal)", font=2)
text(2.15, 8, "Weight\n(Random)", font=2)
This documented is licensed under the following license: CC Attribution-Noncommercial-Share Alike 4.0 International.