cooper2019.rmd

---
title: "R Code Corresponding to the Book *The Handbook of Research Synthesis and Meta-Analysis* by Cooper et al. (2019)"
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 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](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 (e.g., when using MCMC sampling in the chapter on Bayesian methods). 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)
```

***

## 1) Research Synthesis as a Scientific Process

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:

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

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

***

## 10) Evaluating Coding Decisions

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.

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

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

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

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

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

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

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

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

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

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

```{r}
# cell counts and marginal totals for coders 1 and 2

addmargins(table(c2, c1))

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

```{r}
# agreement rate for coders 1 and 2

mean(c1 == c2)
```

```{r}
# agreement rate for all three coders

mean(c1 == c2 & c1 == c3)
```

```{r}
# agreement rate (in %) between coders 1 and 2

irr::agree(c1c2)

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

```{r}
# agreement rate (in %) between all three coders

irr::agree(all3)
```

```{r}
# unweighted Cohen's kappa for coders 1 and 2

irr::kappa2(c1c2)
```

```{r}
# unweighted Cohen's kappa for all three coders

irr::kappam.fleiss(all3)
```

```{r}
# weighted Cohen's kappa for coders 1 and 2

irr::kappa2(c1c2, weight=0:2)
```

We can also use the `psych` package to compute Cohen's kappa, which also provides corresponding confidence intervals.

```{r}
# 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
res <- psych::cohen.kappa(c1c2, w=W)
print(res, digits=3)
```

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.

```{r}
print(vcd::Kappa(table(c1,c2)), digits=3, CI=TRUE)

# note: the (default) weighting scheme used for computing weighted kappa by
# the function is the one described in the chapter
```

```{r}
# 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))
irr::kripp.alpha(t(c1c2), method="ordinal")
irr::kripp.alpha(t(c1c2), method="ratio")
```

```{r}
# correlation between coders 1 and 2

cor(c1, c2)

# note: the cor() function is part of the 'stats' package, which comes with R
```

```{r}
# mean correlation between all pairs of coders

irr::meancor(all3)
```

```{r}
# intraclass correlation coefficient for coders 1 and 2

psych::ICC(c1c2)

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

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

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

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

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

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

***

## 11) Effect Sizes for Meta-Analysis

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

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

### Effect Sizes for a Comparison of Means

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

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

# note: since n1i=n2i in this example, the results are exactly the same
```

```{r}
# mean change

dat <- escalc("MC", m1i=105, m2i=100, sd1i=10, sd2i=10, ni=50, ri=0.5)
summary(dat)
```

```{r}
# standardized mean difference (Hedges' g)

dat <- escalc("SMD", m1i=103, m2i=100, sd1i=5.5, sd2i=4.5, n1i=50, n2i=50)
summary(dat)
```

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

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

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

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

```{r}
# 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)
round(Vg, digits=4)
round(sqrt(Vg), digits=4)

# note: the results given in the book are not quite correct
```

### Correlations

```{r}
# r-to-z transformed correlation coefficient

dat <- escalc("ZCOR", ri=0.50, ni=100)
summary(dat)

# back-transformation

c(transf.ztor(dat$yi))
```

### Effect Sizes for Comparing Risks

```{r}
# risk difference

dat <- escalc("RD", ai=5, n1i=100, ci=10, n2i=100)
summary(dat)
```

```{r}
# risk ratio (log transformed)

dat <- escalc("RR", ai=5, n1i=100, ci=10, n2i=100)
summary(dat)
```

```{r}
# odds ratio (log transformed)

dat <- escalc("OR", ai=5, n1i=100, ci=10, n2i=100)
summary(dat)
```

```{r}
# odds ratio (log transformed) for a case-control study

dat <- escalc("OR", ai=25, bi=20, ci=75, di=80)
summary(dat)
```

***

## 12) Statistically Analyzing Effect Sizes: Equal- and Random-Effects Models

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

```{r}
# equal-effects model analysis

res <- rma(d, v, data=dat, method="EE")
print(res, digits=3)
```

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](https://wviechtb.github.io/metafor/reference/misc-models.html).

```{r}
# random-effects model analysis

res <- rma(d, v, data=dat, method="DL")
print(res, digits=3)

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

```{r}
# fixed-effects ANOVA-type analysis

res <- rma(d, v, mods = ~ factor(group) - 1, data=dat, method="FE")
print(res, digits=3)

# note: by removing the intercept, the three coefficients directly provide the
# estimated average effect for the three groups
```

```{r}
# weighted grand mean effect size

rma(coef(res), diag(vcov(res)), method="FE", digits=3)
```

```{r}
# partitioning of the Q-statistics

res <- rma(d, v, mods = ~ factor(group), data=dat, method="FE")
res

# not removing the intercept, so the QM-statistic is equal to Q-between

round(res$QM, digits=3) # Q-between
round(res$QE, digits=3) # Q-within

# 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)
round(res2$QE, digits=3) # 0.0004 in the book, but must be exactly 0 since k=1
round(res3$QE, digits=3)

# these add up to Q-within above

round(res1$QE + res2$QE + res3$QE, digits=3)
```

```{r}
# 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)
predict(res, newmods=c(-1,0,1), digits=3)

# note: the results given in the book are slightly off
```

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

# maximum likelihood estimate of tau^2

res <- rma(d, v, mods = ~ factor(group), data=dat, method="ML")
round(res$tau2, digits=3)

# restricted maximum likelihood estimate of tau^2

res <- rma(d, v, mods = ~ factor(group), data=dat, method="REML")
round(res$tau2, digits=3)

# note: the REML estimate is incorrectly claimed to be 0 in the book
```

```{r}
res <- rma(d, v, mods = ~ factor(group) - 1, data=dat, method="DL")
print(res, digits=3)

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

# contrast between group 1 and 3

anova(res, X=c(-1,0,1), digits=3)
predict(res, newmods=c(-1,0,1), digits=3)
```

```{r}
# meta-regression model

res <- rma(d, v, mods = ~ I(log(nitems)), data=dat, method="FE")
print(res, digits=3)

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

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

# maximum likelihood estimate of tau^2

res <- rma(d, v, mods = ~ I(log(nitems)), data=dat, method="ML")
round(res$tau2, digits=3)

# restricted maximum likelihood estimate of tau^2

res <- rma(d, v, mods = ~ I(log(nitems)), data=dat, method="REML")
round(res$tau2, digits=3)

# note: the REML estimate is incorrectly claimed to be 0 in the book
```

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

```{r}
# robust variance estimation

robust(res, cluster=study, digits=3)

# note: the test statistic for log(nitems) is somewhat off in the book
```

***

## 13) Stochastically Dependent Effect Sizes

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

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

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

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

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

```{r}
# multivariate fixed-effects model

res <- rma.mv(d, V, mods = ~ outcome - 1, data=dat)
print(res, digits=3)
```

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

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

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

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

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

if (!require(clubSandwich))
   install.packages("clubSandwich")
```

```{r}
robust(res, cluster=study, clubSandwich=TRUE)

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

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

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

```{r}
robu(d ~ outcome - 1, data=dat, studynum=study, var.eff.size=v, rho=0.7, small=TRUE)

# note: the value of rho actually has no influence on these results
```

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

```{r, include=FALSE}
# an example where tau^2 != 0

tmp <- dat
tmp$v <- tmp$v / 10
res <- robu(d ~ outcome - 1, data=tmp, studynum=study, var.eff.size=v, rho=0.7, small=TRUE)
res
V <- lapply(split(tmp$v, tmp$study), vcalcfun, tau2=c(res$mod_info$tau.sq))
res <- rma.mv(d, V, mods = ~ outcome - 1, data=tmp)
robust(res, cluster=study, clubSandwich=TRUE)
```

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

```{r, include=FALSE, eval=FALSE}
# can also do this

dat <- escalc(yi=d, vi=v, data=dat, var.names=c("d","v"))
agg <- aggregate(dat, cluster=study, struct="CS", rho=1, weighted=FALSE, checkpd=FALSE,
                 select=c("study", "d", "v"))
agg
```

```{r}
# fit standard random-effects model to the dataset with synthetic effects

res <- rma(d, v, data=agg, method="DL")
print(res, digits=3)
```

***

## 14) Bayesian Meta-Analysis

For these analyses, we could do most of what is described in the chapter using the excellent [brms](https://cran.r-project.org/package=brms) package. However, to fully reproduce all of the analyses conducted, we either need to use [WinBUGS](https://www.mrc-bsu.cam.ac.uk/software/bugs/the-bugs-project-winbugs/) (as was done by the chapter authors) or we can use [JAGS](http://mcmc-jags.sourceforge.net/) (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](https://cran.r-project.org/package=rjags) package. Note that JAGS needs to be installed separately (follow the link above for installation instructions).

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

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

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

# note: including 'conceal' variable (coded 0 = allocation concealment
# adequate, 1 = allocation concealment inadequate)
```

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

```{r}
# compare with results from a non-Bayesian analysis

res.ee <- rma(yi, vi, data=dat, method="EE", digits=2)
predict(res.ee, transf=exp)
```

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

```{r}
# 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
conf.re <- confint(res.re)
round(conf.re$random["tau^2",], digits=2)
```

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

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

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

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

```{r}
# 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)
round(confint(res.me)$random["tau^2",], digits=2)
```

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

```{r}
# shrunken estimates for Aerdts (1991) and Blair (1991)

round(bma.re[4:5,], digits=2)
```

```{r}
# expected results for a new study

round(bma.re[6:7,], digits=2)
```

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

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

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

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

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

```{r, echo=FALSE, message=FALSE, warning=FALSE}
res.re <- rma(yi, vi, data=dat, method="DL", digits=2)
pred.re <- predict(res.re, transf=exp)
conf.re <- confint(res.re)

k <- length(dat$yi)
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, 0.0001)
   tau ~ dunif(0, 2)
   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.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,])

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,])

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,])

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.01, 1/1.64^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,])

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,]), c(bma.ip[2,], bma.ip[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",
                   "Bayesian RE model, log-normal(-2.01,1.64^2) prior for tau^2")
```

```{r}
# Table 14.5

round(tab, digits=2)

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

***

## 16) Model-Based Meta-Analysis and Related Approaches

```{r}
# Table 16.1: Relationships Between CSAI Subscales and Sport Performance

dat <- dat.craft2003

# examine data for study 3

dat[dat$study == 3,]

# note: the data are in long format; we can recreate Table 16.1 as follows

tab <- do.call(rbind, lapply(split(dat, dat$study), function(x)
   data.frame(x$study[1], x$ni[1], x$sport[1], data.frame(rbind(x$ri)))))
names(tab) <- c("ID", "ni", "type", paste0(dat$var1[1:6], "/", dat$var2[1:6]))
tab
```

The `rcalc()` function from the `metafor` package can be used to compute the variance-covariance (var-cov) matrix of the correlation coefficients. The typical usage would be as follows.

```{r}
tmp <- rcalc(ri ~ var1 + var2 | study, ni=ni, data=dat)
V <- tmp$V
tmp <- tmp$dat

# examine data for study 3

tmp[tmp$study == 3,]

# var-cov matrix for study 3

round(V[tmp$study == 3, tmp$study == 3], digits=3)
```

However, as described in the chapter, it may be advantageous to compute the var-cov matrix based on the (sample-size weighted) mean correlations instead. We can obtain these mean correlations by running a fixed-effects meta-regression model on the data with a moderator that indicates the 6 possible pairs of correlations and setting the weights equal to the sample sizes.

```{r, warning=FALSE}
dat$v. <- paste0(dat$var1, ".", dat$var2)
dat$v. <- factor(dat$v.,
   levels=c("acog.perf", "asom.perf", "conf.perf", "acog.asom", "acog.conf", "asom.conf"))

res <- rma(ri, 0, weights=ni, mods = ~ v. - 1, data=dat, method="FE")
round(coef(summary(res)), digits=3)

# note: the 'estimate' values are the sample-size weighted mean correlations;
# all other results are irrelevant / meaningless (by having set the sampling
# variances to 0 for all outcomes)
```

Next, we need to add these mean correlations back to the dataset. This is easily done by computing the fitted values. However, care must be taken since some rows were excluded from the model fitting (as some of the correlations are missing in studies 6 and 17). To make sure that we get a vector with fitted values that is of the same length as the original data, we need to adjust the `na.action` option.

```{r}
options(na.action = "na.pass")
dat$fit <- fitted(res)
options(na.action = "na.omit")

# examine data for study 3

dat[dat$study == 3,]
```

Now we can compute the var-cov matrix with these mean correlations as outcome. The equations for the variances and covariances given in the chapter use $n_i$ in the denominator, but the `rcalc()` function uses $n_i - 1$, so we just set the sample size to `ni+1` to obtain the same results as in the chapter.

```{r}
tmp <- rcalc(fit ~ var1 + var2 | study, ni=ni+1, data=dat)
V <- tmp$V

# can remove the 'fit' variable from 'dat'

dat$fit <- NULL

# var-cov matrix for study 3

round(V[dat$study == 3, dat$study == 3], digits=3)
```

```{r, warning=FALSE}
# fixed-effects multivariate model

res <- rma.mv(ri, V, mods = ~ v. - 1, data=dat)
res
```

```{r}
# restructure estimated mean correlations into a 4x4 lower triangular matrix

R <- vec2mat(coef(res), dimnames=c("perf", "acog", "asom", "conf"))
round(R, digits=3)
```

```{r}
# var-cov matrix of the estimated mean correlations

round(vcov(res), digits=5)
```

```{r, warning=FALSE}
# BA1 method to test for heterogeneity

res <- lapply(levels(dat$v.), function(l) rma.mv(ri, V, data=dat, subset=v.==l))
ba1 <- data.frame(Q    = sapply(res, function(x) x$QE),
                  df   = sapply(res, function(x) x$k-1),
                  pval = sapply(res, function(x) x$QEp))
ba1$sig <- ifelse(ba1$pval <= 0.05/6, "*", "")
rownames(ba1) <- levels(dat$var1.var2)
dfround(ba1, digits=c(1,0,4,NA))
```

```{r, warning=FALSE}
# random-effects multivariate model

res <- rma.mv(ri, V, mods = ~ v. - 1, random = ~ v. | study, struct="UN", data=dat)
res
```

```{r}
# var-cov matrix of the random effects (values below diagonal are correlations)
# note: the T-hat matrix given in the chapter (p. 354) is in a different order

G <- res$G
G[lower.tri(G)] <- cov2cor(G)[lower.tri(G)]
round(G, digits=2)
```

```{r}
# restructure estimated mean correlations into a 4x4 matrix

R <- vec2mat(coef(res), dimnames=c("perf", "acog", "asom", "conf"))
round(R, digits=3)
```

```{r}
# var-cov matrix of the estimated mean correlations

round(vcov(res), digits=4)
```

```{r}
# fit regression model with 'conf' as outcome and 'acog' and 'asom' as predictors

fit1 <- matreg(4, 2:3, R=R, V=vcov(res))
dfround(fit1$tab[1:4], digits=c(2,3,2,3))

# fit regression model with 'perf' as outcome and 'acog', 'asom', and 'conf' as predictors

fit2 <- matreg(1, 2:4, R=R, V=vcov(res))
dfround(fit2$tab[1:4], digits=c(2,3,2,3))
```

```{r}
# variance explained in sport performance by the three subscales

round(fit2$R2 * 100, digits=2)
```

```{r, include=FALSE, eval=FALSE}
library(metaSEM)
model <- '
conf ~ acog + asom
perf ~ acog + asom + conf
acog ~~ 1*acog
asom ~~ 1*asom
acog ~~ asom'
fit <- wls(R, vcov(res), n=1, RAM=lavaan2RAM(model, obs.variables=rownames(R)))
summary(fit)
```

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

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

```{r, figure16_4, fig.width=8, fig.height=7, dev.args=list(pointsize=pointsize), fig.align='center'}
# Figure 16.4: Random-Effects Path Coefficients for Prediction of Sport Performance

B <- R
B[,] <- 0
B[2:3,4] <- fit1$tab$beta
B[2:4,1] <- fit2$tab$beta
P <- B
P[,] <- 1
P[2:3,4] <- fit1$tab$pval
P[2:4,1] <- fit2$tab$pval

g <- graph_from_adjacency_matrix(B, mode="directed", weighted=TRUE)
vertex_attr(g, "name") <- c("Sport\nperformance", "Cognitive",
                            "Somatic", "Self\nconfidence")
L <- matrix(c(0,0, -1,2, 1,2, 0,1.2), nrow=4, ncol=2, byrow=TRUE)
par(mar=c(0,0,0,0))
plot(g, layout=L,
     edge.color = "black", edge.label.color = "black",
     vertex.color = "black", vertex.label.color = "black",
     vertex.shape = "none", vertex.size = 20,
     edge.label = paste0("              ", fmtx(t(B)[t(B) != 0], 2),
                         ifelse(t(P)[t(B) != 0] <= 0.05, "*", "")),
     edge.lty   = ifelse(t(P)[t(B) != 0] <= 0.05, "solid", "dashed"))
```

```{r, warning=FALSE}
# moderator analysis (subgroup by sport type)

dat.T <- dat[dat$sport == "T",]

res <- rma(ri, 0, weights=ni, mods = ~ v. - 1, data=dat.T, method="FE")
options(na.action = "na.pass")
dat.T$fit <- fitted(res)
options(na.action = "na.omit")

tmp <- rcalc(fit ~ var1 + var2 | study, ni=ni+1, data=dat.T)
V.T <- tmp$V

res.T <- rma.mv(ri, V.T, mods = ~ v. - 1, random = ~ v. | study, struct="UN", data=dat.T)

R.T <- R
R.T[lower.tri(R.T)] <- coef(res.T)

fit.T <- matreg(1, 2:4, R=R.T, V=vcov(res.T))

dat.I <- dat[dat$sport == "I",]

res <- rma(ri, 0, weights=ni, mods = ~ v. - 1, data=dat.I, method="FE")
options(na.action = "na.pass")
dat.I$fit <- fitted(res)
options(na.action = "na.omit")

tmp <- rcalc(fit ~ var1 + var2 | study, ni=ni+1, data=dat.I)
V.I <- tmp$V

res.I <- rma.mv(ri, V.I, mods = ~ v. - 1, random = ~ v. | study, struct="UN", data=dat.I)

R.I <- R
R.I[lower.tri(R.I)] <- coef(res.I)

fit.I <- matreg(1, 2:4, R=R.I, V=vcov(res.I))
```

```{r}
# show results for team sports

dfround(fit.T$tab[1:4], digits=c(2,3,2,3))
```

```{r}
# show results for individual sports

dfround(fit.I$tab[1:4], digits=c(2,3,2,3))
```

```{r}
# var-cov matrix of coefficients for team sports

round(fit.T$vb, digits=4)
```

```{r}
# var-cov matrix of coefficients for individual sports

round(fit.I$vb, digits=4)
```

```{r}
# Wald-type test of the difference in slopes of conf -> perf

rma(yi = c(fit.I$tab["conf", "beta"], fit.T$tab["conf", "beta"]),
    vi = c(fit.I$tab["conf", "se"],   fit.T$tab["conf", "se"])^2,
    mods = c(0,1), method = "FE", digits = 3)

# note: z=-0.094 differs from value in chapter (-0.079) due to rounding
```

```{r, warning=FALSE}
# compute matrix with partial correlations based on R

invR <- solve(R)
Rp <- -diag(diag(1/sqrt(invR))) %*% invR %*% diag(diag(1/sqrt(invR)))
round(Rp, digits=3)
```

```{r, include=FALSE}
# TODO: compute var-cov matrix of the partial correlations based on vcov(res)
```

```{r}
# compute partial correlations for each study individually

Ri <- lapply(split(dat, dat$study), function(x) {
   if (anyNA(x$ri)) {
      NULL
   } else {
      R <- vec2mat(x$ri)
      attr(R, "ni") <- x$ni[1]
      R
   }
})
Ri <- Ri[!sapply(Ri, is.null)]

dat.rp <- lapply(Ri, function(x) {
   invRi <- solve(x)
   rp <- -invRi[1,4] / sqrt(invRi[1,1] * invRi[4,4])
   vrp <- (1-rp^2)^2 / (attributes(x)$ni-3-1)
   data.frame(n=attributes(x)$ni, rp=rp, vrp=vrp)
})
dat.rp <- do.call(rbind, dat.rp)

dfround(dat.rp, digits=c(0,3,4))
```

```{r}
# meta-analysis of the partial correlations

res <- rma(rp, vrp, data=dat.rp, digits=3)
res
```

```{r}
# 95% prediction interval

predict(res, pi.type="simple", digits=2)

# note: can also use predict(res, digits=2), but this will give a slightly
# wider prediction interval because this also considers the uncertainty in
# the estimated mean partial correlation
```

A final note about this chapter: If you are planning on conducting analyses of the types described in this chapter, I would recommend also looking into the excellent [metaSEM](https://cran.r-project.org/package=metaSEM) package, which provides additional functionality for 'meta-analytic structural equation modeling'.

***

## 17) Missing Data in Meta-Analysis

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

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

```{r, warning=FALSE}
# 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)
```

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

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

# note: results in the table are different; not sure why
```

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

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

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

impMethod <- make.method(dat)
impMethod["age"] <- "norm"
impMethod

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]

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

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

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

```{r, warning=FALSE}
# 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)

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

***

## 18) Publication Bias

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.

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

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

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

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

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

```{r}
# Egger regression using a weighted regression model with the standard error as predictor

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

```{r}
# Egger regression using a meta-regression model with the standard error as predictor

regtest(res)
```

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

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

```{r}
# rank correlation test

ranktest(res)
```

For a p-curve analysis, one needs to use the online app at http://p-curve.com/.

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

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

```{r, figure18_15, fig.width=8, fig.height=7, dev.args=list(pointsize=pointsize), fig.align='center'}
# p-uniform

par(mar=c(5,4,2,2))
puniform(yi=dat$yi, vi=dat$vi, side="right", plot=TRUE)
```

```{r}
# excess significance test

tes(yi, vi, data=dat)
```

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

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

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

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

```{r}
# Vevea and Hedges (1995) model (with a single cutpoint at p=0.025)

selmodel(res, type="stepfun", steps=0.025)
```

```{r, include=FALSE}
# Vevea and Hedges (1995) model (with more cutpoints)

selmodel(res, type="stepfun", steps=c(0.01, 0.025, 0.05, 0.30, 0.50, 1.00),
         delta=c(1,NA,NA,NA,0.01,NA))

# note: the value of the selection parameter for the interval that contains no
# p-values needs to be constrained to a fixed value to be able to fit this model
```

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

```{r, warning=FALSE}
# Vevea and Wood (2005) model using moderate one-tailed selection

selmodel(res, type="stepfun", steps=ssp$steps, delta=ssp$weights.mod.1)
```

```{r, warning=FALSE}
# Vevea and Wood (2005) model using severe one-tailed selection

selmodel(res, type="stepfun", steps=ssp$steps, delta=ssp$weights.sev.1)
```

```{r, warning=FALSE}
# Vevea and Wood (2005) model using moderate two-tailed selection

selmodel(res, type="stepfun", steps=ssp$steps, delta=ssp$weights.mod.2)
```

```{r, warning=FALSE}
# Vevea and Wood (2005) model using severe two-tailed selection

selmodel(res, type="stepfun", steps=ssp$steps, delta=ssp$weights.sev.2)
```

```{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, figure18_18, fig.width=9, fig.height=7, dev.args=list(pointsize=pointsize), fig.align='center'}
# Copas and Shi (2001) model

res <- metagen(yi, sqrt(vi), method.tau="ML", data=dat)
sav <- copas(res)
plot(sav)
summary(sav)
```

```{r}
# Rücker limit meta-analysis

sav <- limitmeta(res)
sav
```

***

## 19) Interpreting Effect Sizes

```{r}
# converting d-values to U1, U2, and U3

transf.dtou1(1.0)
transf.dtou2(1.0)
transf.dtou3(1.0)
```

```{r}
# converting d-values to common language effect size values

transf.dtocles(0.3)
```

```{r}
# converting d-values to binomial effect-size display values

transf.dtobesd(0.3)

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

```{r}
# converting a 2x2 table to a correlation (phi coefficient)

escalc(measure="PHI", ai=100, bi=50, ci=60, di=90)
```

```{r}
# converting a 2x2 table to an odds ratio

summary(escalc(measure="OR", ai=100, bi=50, ci=60, di=90), transf=exp)

# note: escalc(measure="OR", ...) computes log odds ratios, so need to exponentiate
```

```{r}
# converting a 2x2 table to a risk ratio

summary(escalc(measure="RR", ai=100, bi=50, ci=60, di=90), transf=exp)

# note: escalc(measure="RR", ...) computes log risk ratios, so need to exponentiate
```

```{r}
# converting a 2x2 table to a risk difference

escalc(measure="RD", ai=100, bi=50, ci=60, di=90)
```

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

```{r}
# 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
predict(res, transf=transf.dtou3)

# 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
predict(res, transf=exp, digits=2)
predict(res, transf=function(x) transf.lnortord(x, pc=0.12), digits=2)
```

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](https://www.russellsage.org/publications/handbook-research-synthesis-and-meta-analysis)).

***

## 20) Heterogeneity in Meta-Analysis

```{r}
# data for the illustrative example

dat <- data.frame(yi = seq(10, 90, by=10), vi = 562.5)
dat
```

```{r}
# fit random-effects model

res <- rma(yi, vi, data=dat, method="DL")
res
```

```{r, figure20_2, fig.width=8, fig.height=7, dev.args=list(pointsize=pointsize), fig.align='center'}
# forest plot

par(mar=c(5,4,2,2))
forest(res, header=TRUE, xlim=c(-170, 310))
```

```{r}
# approximate 95% prediction interval (based on equation 20.6)

predict(res, pi.type="simple", digits=0)
```

```{r}
# approximate 95% prediction interval (based on equation 20.7)

predict(res, pi.type="riley", digits=0)
```

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](https://www.metafor-project.org/doku.php/faq#for_random-effects_models_fitt). However, by using argument `pi.type="riley"`, equation 20.7 is used.

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

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

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

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