Using the `dat.raudenbush1985`

from week 3. Rerun the random-effects meta-analysis with REML estimation.

- Make a funnel plot with study labels. What do you make of it?
- Make a funnel plot with inverse variance on the y-axis.
- Perform Egger’s test for funnel plot asymmetry.
- Perform Begg & Mazumdar’s test for funnel plot asymmetry. Based on
*both*Egger’s test and this test what do you conclude. - Make a contour-enhanced funnel plot. What does it reveal?
- Make a radial plot.
**Bonus**perform and interpret a*p*uniform analysis. For the purpose of this analysis treat the treatment effect`yi`

in this dataset, as a Cohen’s*d*.

There are 19 studies.

`library(meta)`

`## Warning: package 'meta' was built under R version 3.5.2`

```
## Loading 'meta' package (version 4.9-6).
## Type 'help(meta)' for a brief overview.
```

`library(metafor)`

`## Warning: package 'metafor' was built under R version 3.5.2`

`## Loading required package: Matrix`

```
## Loading 'metafor' package (version 2.1-0). For an overview
## and introduction to the package please type: help(metafor).
```

```
##
## Attaching package: 'metafor'
```

```
## The following objects are masked from 'package:meta':
##
## baujat, forest, funnel, funnel.default, labbe, radial, trimfill
```

```
Data<-dat.raudenbush1985
head(Data)
```

```
## study author year weeks setting tester n1i n2i yi vi
## 1 1 Rosenthal et al. 1974 2 group aware 77 339 0.0300 0.0156
## 2 2 Conn et al. 1968 21 group aware 60 198 0.1200 0.0216
## 3 3 Jose & Cody 1971 19 group aware 72 72 -0.1400 0.0279
## 4 4 Pellegrini & Hicks 1972 0 group aware 11 22 1.1800 0.1391
## 5 5 Pellegrini & Hicks 1972 0 group blind 11 22 0.2600 0.1362
## 6 6 Evans & Rosenthal 1969 3 group aware 129 348 -0.0600 0.0106
```

I have made the names as in our slides. Note that this is unnecessary duplication (as TE = *yi* ) but then it maps on nicely to our code.

`library(dplyr)`

`## Warning: package 'dplyr' was built under R version 3.5.2`

```
##
## Attaching package: 'dplyr'
```

```
## The following objects are masked from 'package:stats':
##
## filter, lag
```

```
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
```

`Data <- Data %>% mutate(TE=yi, seTE=sqrt(vi))`

As in the solution of exercise 3, let’s combine, author and year.

Let’s bracket the year as per convention. Here I rely on base R and this snippet.

`Data$year_b <- paste0("(", format(unlist(Data[,3])),")")`

Then we combine as done in here. Here I have opted for a ‘tidyverse’ solution.

`library(tidyverse)`

`## ── Attaching packages ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──`

```
## ✓ ggplot2 3.2.1 ✓ readr 1.3.1
## ✓ tibble 2.1.3 ✓ purrr 0.3.3
## ✓ tidyr 1.0.2 ✓ stringr 1.4.0
## ✓ ggplot2 3.2.1 ✓ forcats 0.4.0
```

`## Warning: package 'ggplot2' was built under R version 3.5.2`

`## Warning: package 'tibble' was built under R version 3.5.2`

`## Warning: package 'tidyr' was built under R version 3.5.2`

`## Warning: package 'purrr' was built under R version 3.5.2`

`## Warning: package 'stringr' was built under R version 3.5.2`

`## Warning: package 'forcats' was built under R version 3.5.2`

```
## ── Conflicts ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x tidyr::expand() masks Matrix::expand()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x tidyr::pack() masks Matrix::pack()
## x tidyr::unpack() masks Matrix::unpack()
```

`Data <-Data %>% unite(author_year, c(author, year_b), sep = " ", remove = FALSE)`

Let’s redo our model but now with `author_year`

.

```
require(meta)
model_reml2<-metagen(TE,
seTE,
data=Data,
studlab=paste(author_year),
comb.fixed = FALSE,
comb.random = TRUE,
method.tau = "REML",
hakn = FALSE,
prediction=TRUE,
sm="SMD")
model_reml2
```

```
## SMD 95%-CI %W(random)
## Rosenthal et al. (1974) 0.0300 [-0.2148; 0.2748] 7.7
## Conn et al. (1968) 0.1200 [-0.1681; 0.4081] 6.6
## Jose & Cody (1971) -0.1400 [-0.4674; 0.1874] 5.7
## Pellegrini & Hicks (1972) 1.1800 [ 0.4490; 1.9110] 1.7
## Pellegrini & Hicks (1972) 0.2600 [-0.4633; 0.9833] 1.7
## Evans & Rosenthal (1969) -0.0600 [-0.2618; 0.1418] 9.1
## Fielder et al. (1971) -0.0200 [-0.2218; 0.1818] 9.1
## Claiborn (1969) -0.3200 [-0.7512; 0.1112] 4.0
## Kester (1969) 0.2700 [-0.0515; 0.5915] 5.8
## Maxwell (1970) 0.8000 [ 0.3081; 1.2919] 3.3
## Carter (1970) 0.5400 [-0.0519; 1.1319] 2.4
## Flowers (1966) 0.1800 [-0.2569; 0.6169] 3.9
## Keshock (1970) -0.0200 [-0.5864; 0.5464] 2.6
## Henrikson (1970) 0.2300 [-0.3384; 0.7984] 2.6
## Fine (1972) -0.1800 [-0.4918; 0.1318] 6.0
## Grieger (1970) -0.0600 [-0.3874; 0.2674] 5.7
## Rosenthal & Jacobson (1968) 0.3000 [ 0.0277; 0.5723] 7.0
## Fleming & Anttonen (1971) 0.0700 [-0.1139; 0.2539] 9.7
## Ginsburg (1970) -0.0700 [-0.4112; 0.2712] 5.4
##
## Number of studies combined: k = 19
##
## SMD 95%-CI z p-value
## Random effects model 0.0837 [-0.0175; 0.1849] 1.62 0.1051
## Prediction interval [-0.2256; 0.3930]
##
## Quantifying heterogeneity:
## tau^2 = 0.0188; H = 1.41 [1.08; 1.84]; I^2 = 49.8% [14.7%; 70.4%]
##
## Test of heterogeneity:
## Q d.f. p-value
## 35.83 18 0.0074
##
## Details on meta-analytical method:
## - Inverse variance method
## - Restricted maximum-likelihood estimator for tau^2
```

Wen note that there is a ‘gap’ in the funnel plot where studies should be. Note that there is one study with a very large effect size (1.18) by Pellegrini and Hicks (1972) (aware group). The study by Maxwell (1970) also has a very large effect size. Note, as mentioned in the slides, that publication bias is but one potential source for funnel plot asymmetry.

`funnel(model_reml2,xlab = "SMD", studlab = TRUE) # adds label on X`

`funnel(model_reml2, yaxis="invvar", main="Inverse Sampling Variance")`

Egger’s test is not significant but it is very, very close *t*(17)= 2.038, *p*= .057 .

`metabias(model_reml2)`

```
##
## Linear regression test of funnel plot asymmetry
##
## data: model_reml2
## t = 2.038, df = 17, p-value = 0.05743
## alternative hypothesis: asymmetry in funnel plot
## sample estimates:
## bias se.bias slope
## 1.6243012 0.7970252 -0.1790311
```

This rank correlation method would lead to a very similar conclusion as above: *p* = .074 (vs. *p* = .057). Note that this test has low statistical power, i.e. one would need a large number of studies. Based on this I would be inclined that the tests do not allow us to rule out publication bias. (Again, note as before, that other processes can cause funnel plot asymmetry and that our tests are based on such asymmetry.)

`metabias(model_reml2,method.bias = "rank")`

```
##
## Rank correlation test of funnel plot asymmetry
##
## data: model_reml2
## z = 1.7865, p-value = 0.07403
## alternative hypothesis: asymmetry in funnel plot
## sample estimates:
## ks se.ks
## 51.0000 28.5482
```

This procedure would have added 3 studies, which would have reduced the effect dramatically. The estimated effect drops from .084 to .028, so its estimate is around 1/3 of the original.

Note that there is a broad discussion on trim & fill procedure, though in our particular case, one would draw the same conclusion: after one adjusts for publication bias, there is even less evidence for a ‘Pygmalion’ effect in this set.

`trimfill(model_reml2)`

```
## SMD 95%-CI %W(random)
## Rosenthal et al. (1974) 0.0300 [-0.2148; 0.2748] 5.9
## Conn et al. (1968) 0.1200 [-0.1681; 0.4081] 5.6
## Jose & Cody (1971) -0.1400 [-0.4674; 0.1874] 5.2
## Pellegrini & Hicks (1972) 1.1800 [ 0.4490; 1.9110] 2.5
## Pellegrini & Hicks (1972) 0.2600 [-0.4633; 0.9833] 2.6
## Evans & Rosenthal (1969) -0.0600 [-0.2618; 0.1418] 6.3
## Fielder et al. (1971) -0.0200 [-0.2218; 0.1818] 6.3
## Claiborn (1969) -0.3200 [-0.7512; 0.1112] 4.4
## Kester (1969) 0.2700 [-0.0515; 0.5915] 5.3
## Maxwell (1970) 0.8000 [ 0.3081; 1.2919] 3.9
## Carter (1970) 0.5400 [-0.0519; 1.1319] 3.3
## Flowers (1966) 0.1800 [-0.2569; 0.6169] 4.3
## Keshock (1970) -0.0200 [-0.5864; 0.5464] 3.4
## Henrikson (1970) 0.2300 [-0.3384; 0.7984] 3.4
## Fine (1972) -0.1800 [-0.4918; 0.1318] 5.4
## Grieger (1970) -0.0600 [-0.3874; 0.2674] 5.2
## Rosenthal & Jacobson (1968) 0.3000 [ 0.0277; 0.5723] 5.7
## Fleming & Anttonen (1971) 0.0700 [-0.1139; 0.2539] 6.4
## Ginsburg (1970) -0.0700 [-0.4112; 0.2712] 5.1
## Filled: Carter (1970) -0.4891 [-1.0809; 0.1028] 3.3
## Filled: Maxwell (1970) -0.7491 [-1.2410; -0.2571] 3.9
## Filled: Pellegrini & Hicks (1972) -1.1291 [-1.8600; -0.3981] 2.5
##
## Number of studies combined: k = 22 (with 3 added studies)
##
## SMD 95%-CI z p-value
## Random effects model 0.0282 [-0.1170; 0.1734] 0.38 0.7035
## Prediction interval [-0.5707; 0.6271]
##
## Quantifying heterogeneity:
## tau^2 = 0.0769; H = 1.67 [1.33; 2.10]; I^2 = 64.3% [43.7%; 77.3%]
##
## Test of heterogeneity:
## Q d.f. p-value
## 58.75 21 < 0.0001
##
## Details on meta-analytical method:
## - Inverse variance method
## - Restricted maximum-likelihood estimator for tau^2
## - Trim-and-fill method to adjust for funnel plot asymmetry
```

```
meta::funnel(model_reml2, xlab="SMD",level = 0.95, contour = c(0.95, 0.975, 0.99),
col.contour = c("darkgray", "gray", "lightgray"),
lwd = 2, cex = 2, pch = 16)
legend(1.25, 0.10,
c("0.05 > p > 0.025", "0.025 > p > 0.01", "< 0.01"),
fill = c("darkgray", "gray", "lightgray"), bg=c("white"))
```

Only a few studies have statistically significant effect sizes (grey background), others do not (white background). Interestingly, it shows that most studies did not produce significant effects. As is clear, there are 2 studies with (very) large effects and low precision (i.e., high standard error). As in the previous plots, there is some indication of publication bias as studies with comparatively lower effect size estimates and low precision are missing (lower quadrant).

`radial(model_reml2)`