Questions.

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

Load and manipulate the data

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

Funnel plot

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

Begg and Mazumdar test

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

Duval & Tweedie’s trim & fill procedure.

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

Contour enhanced funnel plot

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 plot

radial(model_reml2)