rothstein2005.rmd

---
title: "R Code Corresponding to the Book *Publication Bias in Meta-Analysis* by Rothstein et al. (2005)"
author: |
  | Wolfgang Viechtbauer
  | Maastricht University
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
  html_document:
    # code_download: true
    df_print: default
    toc: true
    number_sections: false
    toc_depth: 3
    toc_float:
      collapsed: true
    theme: default
    # lots of nice themes can be used: https://bootswatch.com/
    highlight: haddockadj.theme
  # rmarkdown::github_document
  # pdf_document:
  #   toc: true
  #   number_sections: false
  #   toc_depth: 3
  # word_document
fig_caption: no
# bibliography: references.bib
---

```{r klippy, echo=FALSE, include=TRUE}
# remotes::install_github("rlesur/klippy")
klippy::klippy(position = c("top", "right"), color="gray20")
```

```{r crayon, echo=FALSE, message=FALSE, include=TRUE}
library(crayon)
options(crayon.enabled = TRUE)
knitr::knit_hooks$set(output = function(x, options){
   paste0(
      '<pre class="r-output"><code>',
      fansi::sgr_to_html(x = htmltools::htmlEscape(x), warn = FALSE),
      '</code></pre>'
   )
})
```

## 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](https://www.metafor-project.org/) and the [package documentation](https://wviechtb.github.io/metafor/).

The package can be installed with:

```{r, eval=FALSE}
install.packages("metafor")
```

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

```{r, message=FALSE}
library(metafor)
```

A few additional notes:

1. Results are only reproduced for chapters containing worked examples.
2. Occasionally, there are some minor discrepancies between the results shown in the book and those obtained below. 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).
3. 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.

```{r, echo=FALSE, message=FALSE, warning=FALSE}
setmfopt(space=FALSE)
setmfopt(style=list(legend=make_style("gray90"), warning=strip_style))
pointsize <- 14

options(width=94)
```

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

```{r}
# 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

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

```{r}
# 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
```

```{r}
# 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
```

```{r, figureA_1, fig.width=8.5, fig.height=10, dev.args=list(pointsize=pointsize), fig.align='center'}
# 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

```{r}
# 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

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

```{r}
# random-effects model

res <- rma(yi, vi, data=dat2, method="DL", digits=2, slab=paste0(author, ", ", year))
predict(res, transf=exp)
```

```{r, figureA_2, fig.width=8.5, fig.height=14, dev.args=list(pointsize=pointsize), fig.align='center'}
# 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

```{r}
# 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

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

```{r}
# random-effects model

res <- rma(yi, vi, data=dat3z, digits=2, method="DL")
predict(res, transf=transf.ztor)
```

```{r}
# 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")
```

```{r, figureA_3, fig.width=8.5, fig.height=5, dev.args=list(pointsize=pointsize), fig.align='center'}
# 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

```{r}
# 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
```

```{r}
# fit an equal-effects model

res <- rma(yi, vi, data=dat, method="EE", digits=3)
res
```

```{r, figure5_6, fig.width=8.5, fig.height=10, dev.args=list(pointsize=pointsize), fig.align='center'}
# 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)
```

```{r}
# fit an equal-effects model to dataset 2

res <- rma(yi, vi, data=dat2, method="EE")
res
```

```{r, figure5_7, fig.width=8.5, fig.height=6.5, dev.args=list(pointsize=pointsize), fig.align='center'}
# 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")
```

```{r}
# fit an equal-effects model to dataset 1

res <- rma(yi, vi, data=dat1, method="EE")
res
```

```{r, figure5_8a, fig.width=8.5, fig.height=6.5, dev.args=list(pointsize=pointsize), fig.align='center'}
# 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")
```

```{r, figure5_8b, fig.width=8.5, fig.height=6.5, dev.args=list(pointsize=pointsize), fig.align='center'}
# 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"))
```

```{r}
# fit an equal-effects model to dataset 3r

res <- rma(yi, vi, data=dat3r, method="EE")
res
```

```{r, figure5_9, fig.width=8.5, fig.height=6.5, dev.args=list(pointsize=pointsize), fig.align='center'}
# 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")
```

```{r}
# fit an equal-effects model to dataset 3z

res <- rma(yi, vi, data=dat3z, method="EE")
res
```

```{r, figure5_10a, fig.width=8.5, fig.height=6.5, dev.args=list(pointsize=pointsize), fig.align='center'}
# 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")
```

```{r, figure5_10b, fig.width=8.5, fig.height=6.5, dev.args=list(pointsize=pointsize), fig.align='center'}
# 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

```{r}
# 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)
```

```{r}
# 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 (excluding the study by Bertschat et al., 1989)

ranktest(yi, vi, data=dat, subset=-8, digits=3, exact=FALSE)

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

```{r}
# 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

# bias coefficient (with 95% CI)

round(c(coef=coef(reg$fit)[[2]], confint(reg$fit)[2,]), digits=2)

# rank correlation test (excluding the study by Bertschat et al., 1989)

reg <- regtest(yi, vi, data=dat, subset=-8, model="lm", digits=3)
reg

# bias coefficient (with 95% CI)

round(c(coef=coef(reg$fit)[[2]], confint(reg$fit)[2,]), digits=2)
```

```{r, figure6_1a, fig.width=8.5, fig.height=7, dev.args=list(pointsize=pointsize), fig.align='center'}
# 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)
```

```{r, figure6_1b, fig.width=8.5, fig.height=7, dev.args=list(pointsize=pointsize), fig.align='center'}
# 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.

```{r}
# fit (mixed-effects) meta-regression model with study year as predictor

res <- rma(yi, vi, mods = ~ year, data=dat, subset=-c(8,16))
res
```

```{r}
# regression test (and show results from the model)

regtest(res, ret.fit=TRUE)
```

```{r}
# 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))
```

***

## 7) Failsafe N or File-Drawer Number

```{r}
# failsafe N for dataset 1

fsn(yi, vi, data=dat1)
```

```{r}
# failsafe N for dataset 2

fsn(yi, vi, data=dat2)
```

```{r}
# failsafe N for dataset 3

fsn(yi, vi, data=dat3z)
```

Note: The values obtained above differ slightly from those in the book due to two reasons:

1. Some values were rounded intermittently in the book.
2. 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.

```{r}
# failsafe N for dataset 1 using Orwin's method

fsn(yi, vi, data=dat1, type="Orwin", target=0.05)
```

```{r}
# failsafe N for dataset 2 using Orwin's method

fsn(yi, vi, data=dat2, type="Orwin", target=log(1.05))
```

```{r}
# failsafe N for dataset 3 using Orwin's method

fsn(yi, vi, data=dat3r, type="Orwin", target=0.15)
```

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.

```{r}
# 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)
```

Although the computations underlying Fisher's method are not difficult, we can make use of the [poolr](https://ozancinar.github.io/poolr/) package to automate things.

```{r, message=FALSE, warning=FALSE}
# install the 'poolr' package (if it is not already installed) and load it

if (!require(poolr)) {
   install.packages("poolr")
   library(poolr)
}
```

```{r}
# Fisher's test for dataset 1

pval1 <- c(pnorm(dat1$yi / sqrt(dat1$vi), lower.tail=FALSE))
fisher(pval1)
```

```{r}
# Fisher's test for dataset 2

pval2 <- c(pnorm(dat2$yi / sqrt(dat2$vi), lower.tail=FALSE))
fisher(pval2)
```

```{r}
# Fisher's test for dataset 3

pval3 <- c(pnorm(dat3z$yi / sqrt(dat3z$vi), lower.tail=FALSE))
fisher(pval3)
```

```{r}
# Fisher failsafe N for datasets 1

fsn(yi, vi, data=dat1, pool="Fisher")
```

```{r}
# Fisher failsafe N for datasets 2

fsn(yi, vi, data=dat2, pool="Fisher")
```

```{r}
# Fisher failsafe N for datasets 3

fsn(yi, vi, data=dat3z, pool="Fisher")
```

The results again differ slightly due to the two reasons outlined earlier.

***

## 8) The Trim and Fill Method

```{r}
# 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)
```

```{r}
# random-effects model

res <- rma(yi, sei=sei, method="DL")
res
```

```{r}
# estimated odds ratio (and 95% CI)

predict(res, transf=exp, digits=2)
```

```{r}
# trim and fill analysis

taf <- trimfill(res)
taf
```

```{r}
# estimated odds ratio (and 95% CI)

predict(taf, transf=exp, digits=2)
```

```{r, figure8_2, fig.width=8, fig.height=7, dev.args=list(pointsize=pointsize), fig.align='center'}
# 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)
```

```{r}
# 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))
```

```{r, figure8_3, fig.width=8, fig.height=7, dev.args=list(pointsize=pointsize), fig.align='center'}
# 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)
```

```{r}
# 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
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))
```

```{r}
# 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
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))
```

Note: One could continue to do the same with other subgroups, as illustrated in Table 8.3.

```{r}
# trim and fill analysis for dataset 1

res <- rma(yi, vi=vi, data=dat1, method="DL")
taf <- trimfill(res)
taf

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

Note: One could do the same within the low- and high-contact subgroups, as illustrated in Table 8.4.

***

## 9) Selection Method Approaches

```{r}
# random-effects model for dataset 2 (unadjusted estimates)

res.un <- rma(yi, vi, data=dat2, method="ML")
res.un
```

```{r}
# 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
```

```{r, figure9_1, fig.width=8, fig.height=6.5, dev.args=list(pointsize=pointsize), fig.align='center'}
# 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))
```

```{r}
# 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
```

```{r}
# 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))
```

```{r, message=FALSE, warning=FALSE}
# install the 'metasens' package (if it is not already installed) and load it

if (!require(metasens)) {
   install.packages("metasens")
   library(metasens)
}
```

```{r}
# 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)
```

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.

```{r}
# meta-regression model for dataset 1 (unadjusted estimates)

res.un <- rma(yi, vi, mods = ~ weeks, data=dat1, method="FE")
res.un
```

```{r}
# Vevea and Hedges (1995) model (weight function estimated)

sel <- selmodel(res.un, type="stepfun", steps=c(0.05, 0.30, 1.00))
sel

```

```{r, figure9_3, fig.width=8, fig.height=6.5, dev.args=list(pointsize=pointsize), fig.align='center'}
# 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))
```

```{r}
# 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)
```

As far as I can tell, the `copas()` function from the [metasens](https://cran.r-project.org/package=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

```{r}
# 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)
```

```{r, figure11_1, fig.width=8.5, fig.height=14, dev.args=list(pointsize=pointsize), fig.align='center'}
# 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")))
```

```{r, figure11_2, fig.width=8.5, fig.height=6.5, dev.args=list(pointsize=pointsize), fig.align='center'}
# 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")
```

```{r}
# rank correlation test (note: a one-sided p-value is given in the book)

ranktest(res, digits=3)
```

```{r}
# Egger's test

regtest(res, model="lm", digits=3)
```

```{r}
# failsafe N

fsn(yi, vi, data=dat2)
```

```{r}
# failsafe N (Orwin's method)

fsn(yi, vi, data=dat2, type="Orwin", target=log(1.05))
```

```{r}
# trim and fill analysis

taf <- trimfill(res)
predict(taf, transf=exp, digits=3)
```

```{r, figure11_4, fig.width=8.5, fig.height=6.5, dev.args=list(pointsize=pointsize), fig.align='center'}
# 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))
```

```{r, figure11_5, fig.width=8.5, fig.height=14, dev.args=list(pointsize=pointsize), fig.align='center'}
# 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](http://creativecommons.org/licenses/by-nc-sa/4.0/).