R Code Corresponding to the Book Publication Bias in Meta-Analysis by Rothstein et al. (2005)
Wolfgang Viechtbauer
Maastricht University
Maastricht University
01 October, 2023
General Notes / Setup
The book Publication Bias in Meta-Analysis: Prevention,
Assessment and Adjustments by Rothstein et al. (2005) provides a
very comprehensive treatment of the topic of publication bias. 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 a few 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:
- Results are only reproduced for chapters containing worked examples.
- Occasionally, there are some minor discrepancies between the results shown in the book and those obtained below. These can result from using different software packages that implement methods in slightly different ways, due to intermittent rounding or using a different rounding scheme, or due to chance when the analyses involve some stochastic process. Minor discrepancies will (usually) not be commented on. However, where discrepancies are more substantial, they will be noted (and the reasons for them if they are known).
- 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.
Appendix A: Data Sets
We will actually start with Appendix A, which provides the three datasets used throughout the book for illustrative purposes.
Dataset 1
# data for the meta-analysis on teacher expectancy effects
dat1 <- dat.raudenbush1985
dat1$vi <- c(0.126, 0.147, 0.167, 0.397, 0.371, 0.103, 0.103, 0.221, 0.165, 0.260,
0.307, 0.223, 0.289, 0.291, 0.159, 0.167, 0.139, 0.094, 0.174)^2
dat1$weeks <- ifelse(dat1$weeks > 1, 1, 0) # dummy-code weeks variable into 'high' vs 'low'
dat1## study author year weeks setting tester n1i n2i yi vi
## 1 1 Rosenthal et al. 1974 1 group aware 77 339 0.0300 0.0159
## 2 2 Conn et al. 1968 1 group aware 60 198 0.1200 0.0216
## 3 3 Jose & Cody 1971 1 group aware 72 72 -0.1400 0.0279
## 4 4 Pellegrini & Hicks 1972 0 group aware 11 22 1.1800 0.1576
## 5 5 Pellegrini & Hicks 1972 0 group blind 11 22 0.2600 0.1376
## 6 6 Evans & Rosenthal 1969 1 group aware 129 348 -0.0600 0.0106
## 7 7 Fielder et al. 1971 1 group blind 110 636 -0.0200 0.0106
## 8 8 Claiborn 1969 1 group aware 26 99 -0.3200 0.0488
## 9 9 Kester 1969 0 group aware 75 74 0.2700 0.0272
## 10 10 Maxwell 1970 0 indiv blind 32 32 0.8000 0.0676
## 11 11 Carter 1970 0 group blind 22 22 0.5400 0.0942
## 12 12 Flowers 1966 0 group blind 43 38 0.1800 0.0497
## 13 13 Keshock 1970 0 indiv blind 24 24 -0.0200 0.0835
## 14 14 Henrikson 1970 1 indiv blind 19 32 0.2300 0.0847
## 15 15 Fine 1972 1 group aware 80 79 -0.1800 0.0253
## 16 16 Grieger 1970 1 group blind 72 72 -0.0600 0.0279
## 17 17 Rosenthal & Jacobson 1968 0 group aware 65 255 0.3000 0.0193
## 18 18 Fleming & Anttonen 1971 1 group blind 233 224 0.0700 0.0088
## 19 19 Ginsburg 1970 1 group aware 65 67 -0.0700 0.0303
# note: 'yi' are the standardized mean differences and 'vi' the corresponding sampling variances;
# the sampling variances in 'dat.raudenbush1985' were computed in a slightly different way than in
# the dataset included in the book, so to make sure we can obtain the same results as provided in
# the book, we just overwrite 'vi' with the squared standard errors given in Table A.1# equal-effects model for studies where teachers had more than one week of contact with the students
res.h <- rma(yi, vi, data=dat1, subset=weeks==1, method="EE", digits=2)
res.h## Equal-Effects Model (k = 11)
##
## I^2 (total heterogeneity / total variability): 0.00%
## H^2 (total variability / sampling variability): 0.64
##
## Test for Heterogeneity:
## Q(df = 10) = 6.38, p-val = 0.78
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.02 0.04 -0.52 0.60 -0.10 0.06
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# equal-effects model for studies where teachers had a week or less of contact with the students
res.l <- rma(yi, vi, data=dat1, subset=weeks==0, method="EE", digits=2)
res.l## Equal-Effects Model (k = 8)
##
## I^2 (total heterogeneity / total variability): 32.65%
## H^2 (total variability / sampling variability): 1.48
##
## Test for Heterogeneity:
## Q(df = 7) = 10.39, p-val = 0.17
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.35 0.08 4.41 <.01 0.19 0.50 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Figure A.1
par(mar=c(4,4,2,2))
tmp <- rbind(dat1[dat1$weeks == 1,], dat1[dat1$weeks == 0,])
tmp$weeks <- ifelse(tmp$weeks == 1, "High", "Low")
with(tmp, forest(yi, vi, slab=paste0(author, ", ", year), psize=1, efac=c(0,1),
xlim=c(-5, 3.5), ylim=c(-1,25), header=TRUE, at=seq(-1,1.5,by=0.5),
rows=c(22:12, 9:2), ilab=weeks, ilab.xpos=-1.2, ilab.pos=2))
text(-1.6, 24, "Prior Contact", font=2)
addpoly(res.h, row=11, mlab="High", col="white", font=c(sans=2))
addpoly(res.l, row= 1, mlab="Low", col="white", font=c(sans=2))
res <- rma(yi, vi, data=tmp, method="EE")
addpoly(res, row=-1, mlab="Overall", font=c(sans=2))
abline(h=0)
Dataset 2
# data for the meta-analysis on the effect of environmental tobacco smoke on lung cancer risk
dat2 <- dat.hackshaw1998
head(dat2, 10) # show the first 10 rows## 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
# note: 'yi' are the log odds ratios and 'vi' the corresponding sampling variances; the studies in
# 'dat.hackshaw1998' are ordered by their publication year and hence are not in the same order as in
# Figure A.2, but this is (with some minor rounding discrepancies) the same dataset; also note that
# the 'yi' values are in some chapters referred to as log risk ratios# random-effects model
res <- rma(yi, vi, data=dat2, method="DL", digits=2, slab=paste0(author, ", ", year))
predict(res, transf=exp)## pred ci.lb ci.ub pi.lb pi.ub
## 1.24 1.13 1.36 0.94 1.63
# Figure A.2
par(mar=c(4,4,2,2))
wi <- fmtx(weights(res), digits=2)
forest(res, atransf=exp, psize=1, header=TRUE, xlim=c(-4.5,7), ylim=c(-0.6,40),
at=log(c(0.2, 0.5, 1, 2, 5, 10)), efac=c(0,1),
ilab=wi, ilab.xpos=3.8, ilab.pos=2)
text(3.3, 39, "Weight", font=2)
Dataset 3
# data for the meta-analysis on employment interview scores and job performance
dat3 <- dat.mcdaniel1994
dat3r <- escalc(measure="COR", ri=ri, ni=ni, data=dat3)
dat3z <- escalc(measure="ZCOR", ri=ri, ni=ni, data=dat3)
head(dat3r, 10) # show the first 10 rows## study ni ri type struct yi vi
## 1 1 123 0.00 j s 0.0000 0.0082
## 2 2 95 0.06 p u 0.0600 0.0106
## 3 3 69 0.36 j s 0.3600 0.0111
## 4 4 1832 0.15 j s 0.1500 0.0005
## 5 5 78 0.14 j s 0.1400 0.0125
## 6 6 329 0.06 j s 0.0600 0.0030
## 7 7 153 0.09 j s 0.0900 0.0065
## 8 8 29 0.40 j s 0.4000 0.0252
## 9 9 29 0.39 s s 0.3900 0.0257
## 10 10 157 0.14 s s 0.1400 0.0062
# note: in 'dat3r', 'yi' are the raw correlation coefficients and 'vi' the corresponding sampling
# variances; in 'dat3z', 'yi' are the r-to-z transformed correlation coefficients and 'vi' the
# corresponding sampling variances# random-effects model
res <- rma(yi, vi, data=dat3z, digits=2, method="DL")
predict(res, transf=transf.ztor)## pred ci.lb ci.ub pi.lb pi.ub
## 0.23 0.20 0.26 -0.09 0.51
# fit random-effects models for the various subgroups
res.s <- rma(yi, vi, method="DL", data=dat3z, subset=struct=="s")
res.u <- rma(yi, vi, method="DL", data=dat3z, subset=struct=="u")
res.j <- rma(yi, vi, method="DL", data=dat3z, subset=type=="j" | type=="s")
res.p <- rma(yi, vi, method="DL", data=dat3z, subset=type=="p")# Figure A.3
par(mar=c(4,4,2,2))
yi <- c(coef(res), coef(res.s), coef(res.u), coef(res.j), coef(res.p))
vi <- c(vcov(res), vcov(res.s), vcov(res.u), vcov(res.j), vcov(res.p))
slab <- c("Overall", "Structured", "Unstructured", "Job-Related", "Psychological")
forest(yi, vi, slab=slab, xlim=c(-0.2,0.65), ylim=c(0,10),
atransf=transf.ztor, at=transf.rtoz(c(0,0.1,0.2,0.3,0.4)),
xlab="Correlation Coefficient", header=c("Subgroup", "Correlation [95% CI]"),
refline=NA, efac=0, psize=1, rows=c(7, 5:4, 2:1))
5) The Funnel Plot
# data for the meta-analysis on magnesium treatment in the prevention of death following MI
dat <- dat.egger2001
# calculate log odds ratios and corresponding sampling variances
dat <- escalc(measure="OR", ai=ai, n1i=n1i, ci=ci, n2i=n2i, data=dat)
dat## id study year ai n1i ci n2i yi vi
## 1 1 Morton 1984 1 40 2 36 -0.8303 1.5551
## 2 2 Rasmussen 1986 9 135 23 135 -1.0561 0.1715
## 3 3 Smith 1986 2 200 7 200 -1.2783 0.6531
## 4 4 Abraham 1987 1 48 1 46 -0.0435 2.0435
## 5 5 Feldstedt 1988 10 150 8 148 0.2231 0.2393
## 6 6 Shechter 1989 1 59 9 56 -2.4075 1.1496
## 7 7 Ceremuzynski 1989 1 25 3 23 -1.2809 1.4250
## 8 8 Bertschat 1989 0 22 1 21 -1.1917 2.7599
## 9 9 Singh 1990 6 76 11 75 -0.6957 0.2875
## 10 10 Pereira 1990 1 27 7 27 -2.2083 1.2313
## 11 11 Shechter 1991 2 89 12 80 -2.0382 0.6095
## 12 12 Golf 1991 5 23 13 33 -0.8502 0.3825
## 13 13 Thogersen 1991 4 130 8 122 -0.7932 0.3917
## 14 14 LIMIT-2 1992 90 1159 118 1157 -0.2993 0.0215
## 15 15 Shechter 1995 4 107 17 108 -1.5708 0.3295
## 16 16 ISIS-4 1995 2216 29011 2103 29039 0.0576 0.0010
## Equal-Effects Model (k = 16)
##
## I^2 (total heterogeneity / total variability): 68.13%
## H^2 (total variability / sampling variability): 3.14
##
## Test for Heterogeneity:
## Q(df = 15) = 47.059, p-val < .001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.015 0.031 0.484 0.629 -0.045 0.075
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Figure 5.6: Funnel plots for the magnesium trials with different vertical axes
par(mfrow=c(3,2), mar=c(5,6,2,2), mgp=c(4,1,0), las=1)
funnel(res, yaxis="sei", xlim=c(-4,4), ylim=c(0,2), steps=3, back=NA, lty=1, xlab="")
funnel(res, yaxis="seinv", xlim=c(-4,4), ylim=c(1e-6,32), back=NA, lty=1, xlab="")
funnel(res, yaxis="vi", xlim=c(-4,4), ylim=c(0,3), steps=4, back=NA, lty=1, xlab="")
funnel(res, yaxis="vinv", xlim=c(-4,4), ylim=c(1e-6,1000), steps=3, back=NA, lty=1, xlab="")
funnel(res, yaxis="ni", xlim=c(-4,4), ylim=c(0,60000), steps=4, back=NA, lty=1, xlab="")
mtext("Log odds ratio", side=1, line=3, cex=0.7)
funnel(res, yaxis="lni", xlim=c(-4,4), ylim=c(2,12), steps=6, back=NA, lty=1, xlab="",
yaxt="n", ylab="Sample Size (log scale)")
axis(side=2, at=log(c(10, 100, 1000, 10000, 100000)),
labels=c("10", "100", "1000", "10000", "100000"))
mtext("Log odds ratio", side=1, line=3, cex=0.7)
## 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
# Figure 5.7
par(mar=c(5,4,2,2))
funnel(res, atransf=exp, at=log(c(0.25,0.50,1,2,4,8)), ylim=c(0,0.8), las=1, digits=c(2,1),
back=NA, shade=NA, hlines="lightgray")
## Equal-Effects Model (k = 19)
##
## I^2 (total heterogeneity / total variability): 47.12%
## H^2 (total variability / sampling variability): 1.89
##
## Test for Heterogeneity:
## Q(df = 18) = 34.0384, p-val = 0.0125
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.0578 0.0366 1.5786 0.1144 -0.0140 0.1295
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Figure 5.8a
par(mar=c(5,4,2,2))
funnel(res, at=seq(-1.2,1.2,by=0.4), ylim=c(0,0.4), las=1, digits=1,
back=NA, shade=NA, hlines="lightgray")
# Figure 5.8b
par(mar=c(5,4,2,2))
funnel(res, at=seq(-1.2,1.2,by=0.4), ylim=c(0,0.4), las=1, digits=c(1,1),
back=NA, shade=NA, hlines="lightgray", pch=ifelse(weeks == 1, 19, 17))
legend("topleft", inset=0.02, bg="white", pch=c(19,17), legend=c("High contact", "Low contact"))
## Equal-Effects Model (k = 160)
##
## I^2 (total heterogeneity / total variability): 97.41%
## H^2 (total variability / sampling variability): 38.60
##
## Test for Heterogeneity:
## Q(df = 159) = 6136.9194, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.4713 0.0048 97.9829 <.0001 0.4619 0.4808 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Figure 5.8
funnel(res, at=seq(-1,1,by=0.25), ylim=c(0,0.3), steps=4, las=1, digits=c(2,1),
back=NA, shade=NA, hlines="lightgray")
## Equal-Effects Model (k = 160)
##
## I^2 (total heterogeneity / total variability): 79.87%
## H^2 (total variability / sampling variability): 4.97
##
## Test for Heterogeneity:
## Q(df = 159) = 789.7321, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.2105 0.0064 33.1196 <.0001 0.1980 0.2229 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Figure 5.10a
par(mar=c(5,4,2,2))
funnel(res, at=seq(-1,3,by=0.5), ylim=c(0,0.6), steps=7, las=1, digits=1,
back=NA, shade=NA, hlines="lightgray")
# Figure 5.10b
par(mar=c(5,4,2,2))
sav <- funnel(res, at=seq(-1,3,by=0.5), ylim=c(0,0.6), steps=7, las=1, digits=1,
back=NA, shade=NA, hlines="lightgray", pch=ifelse(struct == "u", 2, 19))
with(subset(sav, dat3z$struct == "u"), points(x, y, pch=24, bg="white"))
legend("topright", inset=0.02, bg="white", pch=c(2,19), legend=c("Unstructured", "Structured"))
6) Regression Methods to Detect Publication and Other Bias in Meta-Analysis
# data for the meta-analysis on magnesium treatment in the prevention of death following MI
dat <- dat.egger2001
# calculate log odds ratios and corresponding sampling variances
dat <- escalc(measure="OR", ai=ai, n1i=n1i, ci=ci, n2i=n2i, data=dat)# rank correlation test excluding ISIS-4 study (and the study by Bertschat et al., 1989)
ranktest(yi, vi, data=dat, subset=-c(8,16), digits=3, exact=FALSE)## Rank Correlation Test for Funnel Plot Asymmetry
##
## Kendall's tau = -0.077, p = 0.702
# rank correlation test (excluding the study by Bertschat et al., 1989)
ranktest(yi, vi, data=dat, subset=-8, digits=3, exact=FALSE)## Rank Correlation Test for Funnel Plot Asymmetry
##
## Kendall's tau = 0.105, p = 0.586
# note: by default, ranktest() computes exact p-values (if k < 50 and there are no ties);
# using exact=FALSE, we can force the use of the normalized test described in the book# regression test excluding ISIS-4 study (and the study by Bertschat et al., 1989)
reg <- regtest(yi, vi, data=dat, subset=-c(8,16), model="lm", digits=3)
reg## Regression Test for Funnel Plot Asymmetry
##
## Model: weighted regression with multiplicative dispersion
## Predictor: standard error
##
## Test for Funnel Plot Asymmetry: t = -3.136, df = 12, p = 0.009
## Limit Estimate (as sei -> 0): b = -0.130 (CI: -0.511, 0.252)
## coef 2.5 % 97.5 %
## -1.27 -2.15 -0.39
# rank correlation test (excluding the study by Bertschat et al., 1989)
reg <- regtest(yi, vi, data=dat, subset=-8, model="lm", digits=3)
reg## Regression Test for Funnel Plot Asymmetry
##
## Model: weighted regression with multiplicative dispersion
## Predictor: standard error
##
## Test for Funnel Plot Asymmetry: t = -5.744, df = 13, p < .001
## Limit Estimate (as sei -> 0): b = 0.103 (CI: 0.029, 0.177)
## coef 2.5 % 97.5 %
## -1.67 -2.29 -1.04
# Figure 6.1a
res <- rma(yi, vi, data=dat, method="EE", subset=-c(8,16))
par(mar=c(5,4,2,2))
funnel(res, atransf=exp, at=log(c(0.05, 0.1, 0.25, 0.5, 1, 2, 4, 8, 16)), xlim=log(c(0.04,16)),
ylim=c(0,1.5), steps=4, las=1, digits=list(2L,1), back=NA, shade=NA, hlines="lightgray",
lty=5, lwd=2)
reg <- regtest(res, model="lm")
se <- seq(0.13, 1.6, length=100)
lines(coef(reg$fit)[1] + coef(reg$fit)[2]*se, se, lwd=2, lty=3)
# Figure 6.1b
res <- rma(yi, vi, data=dat, method="EE", subset=-8)
par(mar=c(5,4,2,2))
funnel(res, atransf=exp, at=log(c(0.05, 0.1, 0.25, 0.5, 1, 2, 4, 8, 16)), xlim=log(c(0.04,16)),
ylim=c(0,1.5), steps=4, las=1, digits=list(2L,1), back=NA, shade=NA, hlines="lightgray",
lty=5, lwd=2)
reg <- regtest(res, model="lm")
se <- seq(0.02, 1.6, length=100)
lines(coef(reg$fit)[1] + coef(reg$fit)[2]*se, se, lwd=2, lty=3)
Unfortunately, I do not have access to the dataset by Linde et al. (1997) and hence cannot reproduce the results in Table 6.1. To illustrate the same principle (i.e., extending the model by including other study characteristics as predictors and using a random/mixed-effects meta-regression model for the analysis and the regression test), I will continue to use the data from the magnesium treatment meta-analysis and include study year as a predictor.
# fit (mixed-effects) meta-regression model with study year as predictor
res <- rma(yi, vi, mods = ~ year, data=dat, subset=-c(8,16))
res## Mixed-Effects Model (k = 14; tau^2 estimator: REML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.2470 (SE = 0.2628)
## tau (square root of estimated tau^2 value): 0.4970
## I^2 (residual heterogeneity / unaccounted variability): 41.70%
## H^2 (unaccounted variability / sampling variability): 1.72
## R^2 (amount of heterogeneity accounted for): 0.00%
##
## Test for Residual Heterogeneity:
## QE(df = 12) = 20.0453, p-val = 0.0662
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 0.1004, p-val = 0.7514
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 50.3481 161.7232 0.3113 0.7556 -266.6236 367.3197
## year -0.0258 0.0813 -0.3169 0.7514 -0.1850 0.1335
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Regression Test for Funnel Plot Asymmetry
##
## Model: mixed-effects meta-regression model
## Predictor: standard error
##
## Mixed-Effects Model (k = 14; tau^2 estimator: REML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.0690 (SE = 0.1875)
## tau (square root of estimated tau^2 value): 0.2626
## I^2 (residual heterogeneity / unaccounted variability): 14.10%
## H^2 (unaccounted variability / sampling variability): 1.16
## R^2 (amount of heterogeneity accounted for): 64.69%
##
## Test for Residual Heterogeneity:
## QE(df = 11) = 10.6881, p-val = 0.4698
##
## Test of Moderators (coefficients 2:3):
## QM(df = 2) = 5.7614, p-val = 0.0561
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 127.7601 139.1082 0.9184 0.3584 -144.8870 400.4072
## year -0.0642 0.0698 -0.9197 0.3577 -0.2011 0.0727
## sei -1.4197 0.5915 -2.4002 0.0164 -2.5790 -0.2604 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Test for Funnel Plot Asymmetry: z = -2.4002, p = 0.0164
# Table 6.2 (rank correlation and regression tests for the three example meta-analyses)
ran1 <- ranktest(yi, vi, data=dat1, exact=FALSE)
reg1 <- regtest(yi, vi, data=dat1, model="lm")
ran2 <- ranktest(yi, vi, data=dat2, exact=FALSE)
reg2 <- regtest(yi, vi, data=dat2, model="lm")
ran3 <- ranktest(yi, vi, data=dat3z, exact=FALSE)
reg3 <- regtest(yi, vi, data=dat3z, model="lm")
tab <- data.frame(
rank_pval = c(ran2$pval, ran1$pval, ran3$pval),
rank_cor = c(ran2$tau, ran1$tau , ran3$tau),
reg_pval = c(reg2$pval, reg1$pval, reg3$pval),
reg_coef = c(coef(reg2$fit)[2], coef(reg1$fit)[2], coef(reg3$fit)[2]),
reg_coef_lb = c(confint(reg2$fit)[2,1], confint(reg1$fit)[2,1], confint(reg3$fit)[2,1]),
reg_coef_ub = c(confint(reg2$fit)[2,2], confint(reg1$fit)[2,2], confint(reg3$fit)[2,2]))
rownames(tab) <- c("dataset 1", "dataset 2", "dataset 3")
dfround(tab, c(3,2,3,2,2,2))## rank_pval rank_cor reg_pval reg_coef reg_coef_lb reg_coef_ub
## dataset 1 0.209 0.14 0.023 0.90 0.13 1.67
## dataset 2 0.086 0.29 0.055 1.58 -0.04 3.20
## dataset 3 0.015 0.13 0.031 0.63 0.06 1.19
7) Failsafe N or File-Drawer Number
## Fail-safe N Calculation Using the Rosenthal Approach
##
## Observed Significance Level: 0.0072
## Target Significance Level: 0.05
##
## Fail-safe N: 24
## Fail-safe N Calculation Using the Rosenthal Approach
##
## Observed Significance Level: <.0001
## Target Significance Level: 0.05
##
## Fail-safe N: 393
## Fail-safe N Calculation Using the Rosenthal Approach
##
## Observed Significance Level: <.0001
## Target Significance Level: 0.05
##
## Fail-safe N: 50364
Note: The values obtained above differ slightly from those in the book due to two reasons:
- Some values were rounded intermittently in the book.
- The
fsn()function uses (approximate) Wald-type tests for computing the \(p\)-values which are used as input for Stouffer’s method, while t-tests were conducted for datasets 1 and 3 in the book.
## Fail-safe N Calculation Using the Orwin Approach
##
## Average Effect Size: 0.0578
## Target Effect Size: 0.0500
##
## Fail-safe N: 3
# failsafe N for dataset 2 using Orwin's method
fsn(yi, vi, data=dat2, 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
## Fail-safe N Calculation Using the Orwin Approach
##
## Average Effect Size: 0.4713
## Target Effect Size: 0.1500
##
## Fail-safe N: 343
Note: As originally described by Orwin (1983), the calculations are
based on unweighted averages, but to replicate the results in
the book, we need to supply the sampling variances as the second
argument to the fsn() function, so that (inverse-variance)
weighted averages are used. However, note that for dataset 3, the mean
correlation used as input for Orwin’s method in the book (i.e., 0.21) is
based on a random-effects model and hence the result above differs from
what is shown in the book. With a simple trick, we can however reproduce
the results in the book exactly. For this, we just create a vector of
the same length as the original data with the mean effect size used as
input and pass this to the fsn() function.
# failsafe N using Orwin's method calculated exactly as described in the book
fsn1 <- fsn(rep(0.06, nrow(dat1)), type="Orwin", target=0.05)$fsn
fsn2 <- fsn(rep(0.185, nrow(dat2)), type="Orwin", target=0.049)$fsn
fsn3 <- fsn(rep(0.21, nrow(dat3)), type="Orwin", target=0.15)$fsn
c(fsn1 = fsn1, fsn2 = fsn2, fsn3 = fsn3)## fsn1 fsn2 fsn3
## 4 103 64
Although the computations underlying Fisher’s method 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)
}# Fisher's test for dataset 1
pval1 <- c(pnorm(dat1$yi / sqrt(dat1$vi), lower.tail=FALSE))
fisher(pval1)## combined p-values with: Fisher's method
## number of p-values combined: 19
## test statistic: 70.125 ~ chi-square(df = 38)
## adjustment: none
## combined p-value: 0.001159823
# Fisher's test for dataset 2
pval2 <- c(pnorm(dat2$yi / sqrt(dat2$vi), lower.tail=FALSE))
fisher(pval2)## combined p-values with: Fisher's method
## number of p-values combined: 37
## test statistic: 158.246 ~ chi-square(df = 74)
## adjustment: none
## combined p-value: 4.587673e-08
# Fisher's test for dataset 3
pval3 <- c(pnorm(dat3z$yi / sqrt(dat3z$vi), lower.tail=FALSE))
fisher(pval3)## combined p-values with: Fisher's method
## number of p-values combined: 160
## test statistic: 2418.5 ~ chi-square(df = 320)
## adjustment: none
## combined p-value: 3.477446e-318
## Fail-safe N Calculation Using the Rosenthal Approach
##
## Observed Significance Level: <.0001
## Target Significance Level: 0.05
##
## Fail-safe N: 37
## Fail-safe N Calculation Using the Rosenthal Approach
##
## Observed Significance Level: <.0001
## Target Significance Level: 0.05
##
## Fail-safe N: 100
## Fail-safe N Calculation Using the Rosenthal Approach
##
## Observed Significance Level: <.0001
## Target Significance Level: 0.05
##
## Fail-safe N: 3152
The results again differ slightly due to the two reasons outlined earlier.
8) The Trim and Fill Method
# log odds ratios and corresponding standard errors as given in Table 8.1
yi <- c(-0.20, -0.07, 0.04, 0.16, 0.21, 0.27, 0.53, 0.56, 0.80, 1.08, 2.11)
sei <- c( 0.41, 0.18, 0.30, 0.53, 0.51, 0.33, 0.74, 1.08, 0.62, 0.66, 1.55)## Random-Effects Model (k = 11; tau^2 estimator: DL)
##
## tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.0807)
## 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 = 10) = 7.4675, p-val = 0.6807
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.1081 0.1174 0.9204 0.3574 -0.1221 0.3383
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## pred ci.lb ci.ub pi.lb pi.ub
## 1.11 0.89 1.40 0.89 1.40
## Estimated number of missing studies on the left side: 5 (SE = 2.0410)
##
## Random-Effects Model (k = 16; tau^2 estimator: DL)
##
## tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.0844)
## 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 = 15) = 14.8758, p-val = 0.4604
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.0147 0.1115 0.1322 0.8948 -0.2037 0.2332
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## pred ci.lb ci.ub pi.lb pi.ub
## 1.01 0.82 1.26 0.82 1.26
# Figure 8.2
par(mar=c(5,4,2,2))
funnel(taf, yaxis="seinv", xlab="ln(Odds Ratio)", back=NA, las=1,
xlim=c(-3,3), ylim=c(1e-6,6), steps=7)
# trim and fill analysis for dataset 2
res <- rma(yi, vi=vi, data=dat2, method="DL")
taf <- trimfill(res)
pred.res <- predict(res, transf=exp)
pred.taf <- predict(taf, transf=exp)
tab <- data.frame(
k = c(res$k, taf$k),
k0 = c(0, taf$k0),
theta = c(pred.res$pred, pred.taf$pred),
ci.lb = c(pred.res$ci.lb, pred.taf$ci.lb),
ci.ub = c(pred.res$ci.ub, pred.taf$ci.ub))
rownames(tab) <- c("Observed", "Filled")
dfround(tab, c(0,0,2,2,2))## k k0 theta ci.lb ci.ub
## Observed 37 0 1.24 1.13 1.36
## Filled 44 7 1.19 1.08 1.31
# Figure 8.3
par(mar=c(5,4,2,2))
funnel(taf, yaxis="seinv", xlab="ln(Odds Ratio)", back=NA, las=1,
xlim=c(-0.8,0.9), at=c(-0.8, -0.4, 0, 0.4, 0.8), ylim=c(1e-6,12), steps=7)
# trim and fill analysis for dataset 3 - structured interviews
res <- rma(yi, vi=vi, data=subset(dat3z, struct=="s"), method="DL")
taf <- trimfill(res)
taf## Estimated number of missing studies on the left side: 19 (SE = 6.7575)
##
## Random-Effects Model (k = 125; tau^2 estimator: DL)
##
## tau^2 (estimated amount of total heterogeneity): 0.0637 (SE = 0.0162)
## tau (square root of estimated tau^2 value): 0.2523
## I^2 (total heterogeneity / total variability): 87.13%
## H^2 (total variability / sampling variability): 7.77
##
## Test for Heterogeneity:
## Q(df = 124) = 963.6530, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.2053 0.0263 7.7997 <.0001 0.1537 0.2569 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pred.res <- predict(res, transf=transf.ztor)
pred.taf <- predict(taf, transf=transf.ztor)
tab <- data.frame(
k = c(res$k, taf$k),
k0 = c(0, taf$k0),
theta = c(pred.res$pred, pred.taf$pred),
ci.lb = c(pred.res$ci.lb, pred.taf$ci.lb),
ci.ub = c(pred.res$ci.ub, pred.taf$ci.ub))
rownames(tab) <- c("Observed", "Filled")
dfround(tab, c(0,0,3,3,3))## k k0 theta ci.lb ci.ub
## Observed 106 0 0.266 0.223 0.308
## Filled 125 19 0.202 0.153 0.251
# trim and fill analysis for dataset 3 - unstructured interviews
res <- rma(yi, vi=vi, data=subset(dat3z, struct=="u"), method="DL")
taf <- trimfill(res)
taf## Estimated number of missing studies on the left side: 0 (SE = 3.3481)
##
## Random-Effects Model (k = 39; tau^2 estimator: DL)
##
## tau^2 (estimated amount of total heterogeneity): 0.0123 (SE = 0.0078)
## tau (square root of estimated tau^2 value): 0.1107
## I^2 (total heterogeneity / total variability): 69.35%
## H^2 (total variability / sampling variability): 3.26
##
## Test for Heterogeneity:
## Q(df = 38) = 123.9901, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.1891 0.0257 7.3546 <.0001 0.1387 0.2395 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pred.res <- predict(res, transf=transf.ztor)
pred.taf <- predict(taf, transf=transf.ztor)
tab <- data.frame(
k = c(res$k, taf$k),
k0 = c(0, taf$k0),
theta = c(pred.res$pred, pred.taf$pred),
ci.lb = c(pred.res$ci.lb, pred.taf$ci.lb),
ci.ub = c(pred.res$ci.ub, pred.taf$ci.ub))
rownames(tab) <- c("Observed", "Filled")
dfround(tab, c(0,0,3,3,3))## k k0 theta ci.lb ci.ub
## Observed 39 0 0.187 0.138 0.235
## Filled 39 0 0.187 0.138 0.235
Note: One could continue to do the same with other subgroups, as illustrated in Table 8.3.
# trim and fill analysis for dataset 1
res <- rma(yi, vi=vi, data=dat1, method="DL")
taf <- trimfill(res)
taf## Estimated number of missing studies on the left side: 3 (SE = 2.9465)
##
## Random-Effects Model (k = 22; tau^2 estimator: DL)
##
## tau^2 (estimated amount of total heterogeneity): 0.0474 (SE = 0.0266)
## tau (square root of estimated tau^2 value): 0.2176
## I^2 (total heterogeneity / total variability): 61.80%
## H^2 (total variability / sampling variability): 2.62
##
## Test for Heterogeneity:
## Q(df = 21) = 54.9740, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.0273 0.0636 0.4282 0.6685 -0.0975 0.1520
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tab <- data.frame(
k = c(res$k, taf$k),
k0 = c(0, taf$k0),
theta = c(coef(res), coef(taf)),
ci.lb = c(res$ci.lb, taf$ci.lb),
ci.ub = c(res$ci.ub, taf$ci.ub))
rownames(tab) <- c("Observed", "Filled")
dfround(tab, c(0,0,3,3,3))## k k0 theta ci.lb ci.ub
## Observed 19 0 0.084 -0.023 0.191
## Filled 22 3 0.027 -0.097 0.152
Note: One could do the same within the low- and high-contact subgroups, as illustrated in Table 8.4.
9) Selection Method Approaches
# random-effects model for dataset 2 (unadjusted estimates)
res.un <- rma(yi, vi, data=dat2, method="ML")
res.un## Random-Effects Model (k = 37; tau^2 estimator: ML)
##
## tau^2 (estimated amount of total heterogeneity): 0.0204 (SE = 0.0165)
## tau (square root of estimated tau^2 value): 0.1427
## I^2 (total heterogeneity / total variability): 27.62%
## H^2 (total variability / sampling variability): 1.38
##
## Test for Heterogeneity:
## Q(df = 36) = 47.4979, p-val = 0.0952
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.2171 0.0486 4.4712 <.0001 0.1219 0.3123 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Vevea and Hedges (1995) model (weight function estimated)
sel <- selmodel(res.un, type="stepfun", steps=c(0.05, 0.10, 0.50, 1.00))
sel## Random-Effects Model (k = 37; tau^2 estimator: ML)
##
## tau^2 (estimated amount of total heterogeneity): 0.0307 (SE = 0.0271)
## tau (square root of estimated tau^2 value): 0.1753
##
## Test for Heterogeneity:
## LRT(df = 1) = 3.9537, p-val = 0.0468
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.1217 0.1299 0.9373 0.3486 -0.1328 0.3763
##
## Test for Selection Model Parameters:
## LRT(df = 3) = 7.0661, p-val = 0.0698
##
## Selection Model Results:
##
## k estimate se zval pval ci.lb ci.ub
## 0 < p <= 0.05 7 1.0000 --- --- --- --- ---
## 0.05 < p <= 0.1 8 2.4221 1.6609 0.8562 0.3919 0.0000 5.6773
## 0.1 < p <= 0.5 16 0.9775 0.8204 -0.0274 0.9782 0.0000 2.5855
## 0.5 < p <= 1 6 0.3967 0.4692 -1.2857 0.1986 0.0000 1.3164
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Figure 9.1
par(mar=c(5,4,2,2))
plot(sel, scale=TRUE, rug=FALSE, bty="n", xaxt="n")
axis(side=1, at=c(0,sel$steps))
# Table 9.3: Specification of weights for a priori weight functions
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),
weak.1 = c(1, 0.99, 0.95, 0.90, 0.80, 0.75, 0.65, 0.60, 0.55, 0.50, 0.50, 0.50, 0.50, 0.50),
stro.1 = c(1, 0.99, 0.90, 0.75, 0.60, 0.55, 0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.20, 0.10),
weak.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),
stro.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 weak.1 stro.1 weak.2 stro.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.90 0.75 0.90 0.75
## 5 0.250 0.80 0.60 0.80 0.60
## 6 0.350 0.75 0.55 0.75 0.50
## 7 0.500 0.65 0.50 0.60 0.25
## 8 0.650 0.60 0.45 0.60 0.25
## 9 0.750 0.55 0.40 0.75 0.50
## 10 0.900 0.50 0.35 0.80 0.60
## 11 0.950 0.50 0.30 0.90 0.75
## 12 0.990 0.50 0.25 0.95 0.90
## 13 0.995 0.50 0.20 0.99 0.99
## 14 1.000 0.50 0.10 1.00 1.00
# apply a priori weight functions
res.weak.1 <- selmodel(res.un, type="stepfun", steps=ssp$steps, delta=ssp$weak.1)
res.stro.1 <- selmodel(res.un, type="stepfun", steps=ssp$steps, delta=ssp$stro.1)
res.weak.2 <- selmodel(res.un, type="stepfun", steps=ssp$steps, delta=ssp$weak.2)
res.stro.2 <- selmodel(res.un, type="stepfun", steps=ssp$steps, delta=ssp$stro.2)
# Table 9.4 (but transposed)
tab <- data.frame(logOR = c(coef(res.un), coef(res.weak.1)$beta, coef(res.stro.1)$beta,
coef(res.weak.2)$beta, coef(res.stro.2)$beta),
varcomp = c(res.un$tau2, res.weak.1$tau2, res.stro.1$tau2,
res.weak.2$tau2, res.stro.2$tau2))
rownames(tab) <- c("unadjusted", "weak one-tailed selection", "strong one-tailed selection",
"weak two-tailed selection", "strong two-tailed selection")
dfround(tab, c(2,3))## logOR varcomp
## unadjusted 0.22 0.020
## weak one-tailed selection 0.16 0.019
## strong one-tailed selection 0.13 0.018
## weak two-tailed selection 0.19 0.015
## strong two-tailed selection 0.15 0.009
# 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 (selection depending on both T and sigma)
res <- metagen(yi, sqrt(vi), method.tau="ML", data=dat2)
sav <- copas(res)
#plot(sav)
summary(sav)## Copas selection model analysis
##
## p.publ 95%-CI tau^2 tau p.trt p.rsb N
## 1.0000 0.2171 [ 0.1182; 0.3160] 0.0204 0.1427 < 0.0001 0.0371 0
## 0.9296 0.2003 [ 0.1010; 0.2996] 0.0189 0.1374 < 0.0001 0.0641 2
## 0.8612 0.1801 [ 0.0820; 0.2782] 0.0184 0.1358 0.0003 0.1142 4
## 0.8068 0.1602 [ 0.0648; 0.2556] 0.0182 0.1349 0.0010 0.1657 5
## 0.6867 0.1394 [ 0.0226; 0.2562] 0.0163 0.1276 0.0193 0.2431 11
## 0.5989 0.1197 [-0.0036; 0.2429] 0.0153 0.1236 0.0571 0.3482 16
## 0.4901 0.0991 [-0.0356; 0.2337] 0.0142 0.1192 0.1494 0.4843 26
##
## Adjusted estimate 0.1801 [ 0.0820; 0.2782] 0.0184 0.1358 0.0003 0.1142 4
## Unadjusted estimate 0.2171 [ 0.1219; 0.3123] 0.0204 0.1427 < 0.0001
##
## Significance level for test of residual selection bias: 0.1
##
## min max
## range of gamma0: -0.4535 2.0000
## range of gamma1: 0.0000 0.1496
##
## Largest standard error (SE): 0.735
##
## Range of probability publishing trial with largest SE:
## min max
## 0.3251 0.9862
##
## Calculation of orthogonal line:
##
## level nobs adj.r.square slope se.slope
## 0.04 6 0.9916557 -1.1278233 0.04622816
## 0.06 13 0.9868226 -1.6880865 0.05628042
## 0.08 19 0.9873009 -2.1817836 0.05830170
## 0.10 27 0.9900491 -2.6612960 0.05231498
## 0.12 26 0.9975215 -3.3608319 0.03350338
## 0.14 26 0.9542003 -3.4079163 0.14918114
## 0.14 5 -0.2955325 0.1343782 0.45419392
## 0.16 25 0.9855488 -4.7644080 0.11772888
## 0.18 23 0.9927567 -4.6554077 0.08476573
## 0.20 23 0.9921361 -5.1308523 0.09737182
##
## 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
Note that this analysis is a bit different than the one presented in the book (where parameters \(a\) and \(b\) are fixed to specific values).
The analyses of dataset 3 (on the validity of the employment interview) presented in this chapter cannot be replicated because the meta-regression model considered here includes a moderator that is not actually part of the dataset provided in Appendix A (whether the interview was conducted for research or for administrative reasons). Hence, we will move on to the analysis of dataset 1.
# meta-regression model for dataset 1 (unadjusted estimates)
res.un <- rma(yi, vi, mods = ~ weeks, data=dat1, method="FE")
res.un## Fixed-Effects with Moderators Model (k = 19)
##
## I^2 (residual heterogeneity / unaccounted variability): 0.00%
## H^2 (unaccounted variability / sampling variability): 0.99
## R^2 (amount of heterogeneity accounted for): 47.82%
##
## Test for Residual Heterogeneity:
## QE(df = 17) = 16.7735, p-val = 0.4698
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 17.2649, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.3494 0.0792 4.4143 <.0001 0.1943 0.5046 ***
## weeks -0.3709 0.0893 -4.1551 <.0001 -0.5459 -0.1960 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Vevea and Hedges (1995) model (weight function estimated)
sel <- selmodel(res.un, type="stepfun", steps=c(0.05, 0.30, 1.00))
sel## Fixed-Effects with Moderators Model (k = 19)
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 11.0191, p-val = 0.0009
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.3206 0.1115 2.8745 0.0040 0.1020 0.5391 **
## weeks -0.3625 0.1092 -3.3195 0.0009 -0.5766 -0.1485 ***
##
## Test for Selection Model Parameters:
## LRT(df = 2) = 0.4323, p-val = 0.8056
##
## Selection Model Results:
##
## k estimate se zval pval ci.lb ci.ub
## 0 < p <= 0.05 4 1.0000 --- --- --- --- ---
## 0.05 < p <= 0.3 6 0.8924 0.7987 -0.1347 0.8928 0.0000 2.4578
## 0.3 < p <= 1 9 0.5603 0.6422 -0.6846 0.4936 0.0000 1.8191
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Figure 9.3
par(mar=c(5,4,2,2))
plot(sel, scale=TRUE, rug=FALSE, bty="n", xaxt="n")
axis(side=1, at=c(0,sel$steps))
# apply a priori weight functions
res.weak.1 <- selmodel(res.un, type="stepfun", steps=ssp$steps, delta=ssp$weak.1)
res.stro.1 <- selmodel(res.un, type="stepfun", steps=ssp$steps, delta=ssp$stro.1)
res.weak.2 <- selmodel(res.un, type="stepfun", steps=ssp$steps, delta=ssp$weak.2)
res.stro.2 <- selmodel(res.un, type="stepfun", steps=ssp$steps, delta=ssp$stro.2)
# Table 9.10 (but transposed)
tab <- data.frame(logOR = rbind(coef(res.un), coef(res.weak.1)$beta, coef(res.stro.1)$beta,
coef(res.weak.2)$beta, coef(res.stro.2)$beta))
colnames(tab) <- names(coef(res.un))
rownames(tab) <- c("unadjusted", "weak one-tailed selection", "strong one-tailed selection",
"weak two-tailed selection", "strong two-tailed selection")
dfround(tab, 2)## intrcpt weeks
## unadjusted 0.35 -0.37
## weak one-tailed selection 0.32 -0.37
## strong one-tailed selection 0.30 -0.36
## weak two-tailed selection 0.33 -0.35
## strong two-tailed selection 0.29 -0.31
As far as I can tell, the copas() function from the metasens package
cannot handle meta-regression models. Hence, this part cannot be
reproduced (without implementing the methods described on pages
157-158).
11) Software for Publication Bias
# fit an equal-effects model to dataset 2 and obtain the estimated odds ratio (with 95% CI)
res <- rma(yi, vi, data=dat2, method="EE", slab=paste0(author, ", ", year))
predict(res, transf=exp, digits=3)## pred ci.lb ci.ub
## 1.204 1.119 1.295
# Figure 11.1
par(mar=c(4,4,2,2))
wi <- fmtx(weights(res), digits=2)
forest(res, atransf=exp, psize=1, header=TRUE, xlim=c(-4.5,7), ylim=c(-0.6,40),
at=log(c(0.2, 0.5, 1, 2, 5, 10)), efac=c(0,1),
ilab=wi, ilab.xpos=3.3, ilab.pos=2, order="prec")
text(2.8, 39, "Weight", font=2)
wi <- rev(weights(res)[order(dat2$vi)])
invisible(sapply(seq_along(wi), function(i) rect(3.5, i-0.35, 3.5+wi[i]/20, i+0.35, col="black")))
# Figure 11.2
par(mar=c(5,4,2,2))
funnel(res, xlim=c(-2,2), ylim=c(0,0.8), las=1, digits=1,
back=NA, shade=NA, hlines=NA, lty=1, pch=21, bg="white")
## Rank Correlation Test for Funnel Plot Asymmetry
##
## Kendall's tau = 0.144, p = 0.216
## Regression Test for Funnel Plot Asymmetry
##
## Model: weighted regression with multiplicative dispersion
## Predictor: standard error
##
## Test for Funnel Plot Asymmetry: t = 2.382, df = 35, p = 0.023
## Limit Estimate (as sei -> 0): b = 0.005 (CI: -0.169, 0.179)
## Fail-safe N Calculation Using the Rosenthal Approach
##
## Observed Significance Level: <.0001
## Target Significance Level: 0.05
##
## Fail-safe N: 393
## Fail-safe N Calculation Using the Orwin Approach
##
## Average Effect Size: 0.1858
## Target Effect Size: 0.0488
##
## Fail-safe N: 104
## pred ci.lb ci.ub
## 1.168 1.088 1.255
# Figure 11.4
par(mar=c(5,4,2,2))
funnel(taf, xlim=c(-2,2), ylim=c(0,0.8), las=1, digits=1,
back=NA, shade=NA, hlines=NA, lty=1,
pch=21, bg=c("white", "gray10"), legend=list(show=NA))
# Figure 11.5: Cumulative forest plot
sav <- cumul(res, order=vi)
wi <- cumsum(weights(res)[order(dat2$vi)])
wi <- fmtx(wi, digits=2)
par(mar=c(4,4,2,2))
forest(sav, atransf=exp, psize=1, header=TRUE, xlim=c(-2,3), ylim=c(1,40),
at=log(c(0.5, 1, 2)), efac=c(0,1), ilab=wi, ilab.xpos=1.4, ilab.pos=2)
text(1.2, 39, "Weight", font=2)
wi <- rev(cumsum(weights(res)[order(dat2$vi)]))
invisible(sapply(seq_along(wi), function(i) rect(1.5, i-0.35, 1.5+wi[i]/250, i+0.35, col="black")))
License
This documented is licensed under the following license: CC Attribution-Noncommercial-Share Alike 4.0 International.