Wolfgang Viechtbauer

Maastricht University

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 several other packages will also be used. To read more about the `metafor`

package, see the package website and the package documentation.

Note that the 'devel' version of `metafor`

needs to be installed as some of the datasets used below are not currently in the official release of the package on CRAN. The 'devel' version of the package can be installed with:

```
install.packages("remotes")
remotes::install_github("wviechtb/metafor")
```

This step will become obsolete once a new release of the `metafor`

package is published on CRAN.

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

`library(metafor)`

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

Finally, let's create a little helper function for formatting some results later on (essentially like `round()`

, but this one does not drop trailing zeros).

```
fc <- function(x, digits=4, ...)
formatC(x, format="f", digits=digits, ...)
```

We will actually start with Appendix A, which provides the three datasets used throughout the book for illustrative purposes.

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

```
# fixed-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="FE", digits=2)
res.h
```

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

```
# fixed-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="FE", digits=2)
res.l
```

```
## Fixed-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")
forest(tmp$yi, tmp$vi, slab=paste0(tmp$author, ", ", tmp$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=tmp$weeks, ilab.xpos=-1.2, ilab.pos=2)
text(-1.6, 24, "Prior Contact", font=2)
addpoly(res.h, row=11, mlab="High", cex=1, width=c(4,5,3), col="white", font=c(sans=2))
addpoly(res.l, row= 1, mlab="Low", cex=1, width=c(4,5,3), col="white", font=c(sans=2))
res <- rma(yi, vi, data=tmp, method="FE")
addpoly(res, row=-1, mlab="Overall", cex=1, width=c(4,5,3), font=c(sans=2))
abline(h=0)
```