Questions.

The data for this question are standardized mean differences from a meta-analysis of 49 experimental effects of teacher expectancy on pupil IQ (Raudenbush, 1984). They are part of the metaforpackage and you can find out more about them in the help file in the metafor package.

The typical study in this dataset administered a pre-test to a sample of students. Subsequently, the teachers were told that a randomly selected subsample were ‘bloomers’ (students with substantial potential intellectual growth). All of the students were then administered a post-test and it was expected that the ones identified as bloomers would show a significantly higher increment in IQ growth than the control group (i.e., students not identified as bloomers). You might have heard of the Pygmalion effect, (see here)

Inspect the data and note that it has yi= standardized mean difference and vi= corresponding sampling variance.

Load and manipulate the data

There are 19 studies.

library(meta)
## Loading 'meta' package (version 4.9-5).
## Type 'help(meta)' for a brief overview.
library(metafor)
## 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 
## 1     1   Rosenthal et al. 1974     2   group  aware  77 339  0.0300 
## 2     2        Conn et al. 1968    21   group  aware  60 198  0.1200 
## 3     3        Jose & Cody 1971    19   group  aware  72  72 -0.1400 
## 4     4 Pellegrini & Hicks 1972     0   group  aware  11  22  1.1800 
## 5     5 Pellegrini & Hicks 1972     0   group  blind  11  22  0.2600 
## 6     6  Evans & Rosenthal 1969     3   group  aware 129 348 -0.0600 
##       vi 
## 1 0.0156 
## 2 0.0216 
## 3 0.0279 
## 4 0.1391 
## 5 0.1362 
## 6 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)
## 
## 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))

Fixed effects model

Fixed<-metagen(TE,
        seTE,
        data=Data,
        studlab=paste(author),
        comb.fixed = TRUE,
        comb.random = FALSE,
        prediction=TRUE,
        sm="SMD")

These results map on to what is reported in Raudenbush et al. (2009): the effect size estimate is .06, with a 95 %CI from -.01 ,z = 1.66 , p = .098. Note that there is significant heterogeneity, e.g., Q(18)= 35.83, p = .007.

Fixed
##                          SMD            95%-CI %W(fixed)
## Rosenthal et al.      0.0300 [-0.2148; 0.2748]       8.5
## Conn et al.           0.1200 [-0.1681; 0.4081]       6.2
## Jose & Cody          -0.1400 [-0.4674; 0.1874]       4.8
## Pellegrini & Hicks    1.1800 [ 0.4490; 1.9110]       1.0
## Pellegrini & Hicks    0.2600 [-0.4633; 0.9833]       1.0
## Evans & Rosenthal    -0.0600 [-0.2618; 0.1418]      12.5
## Fielder et al.       -0.0200 [-0.2218; 0.1818]      12.5
## Claiborn             -0.3200 [-0.7512; 0.1112]       2.7
## Kester                0.2700 [-0.0515; 0.5915]       4.9
## Maxwell               0.8000 [ 0.3081; 1.2919]       2.1
## Carter                0.5400 [-0.0519; 1.1319]       1.5
## Flowers               0.1800 [-0.2569; 0.6169]       2.7
## Keshock              -0.0200 [-0.5864; 0.5464]       1.6
## Henrikson             0.2300 [-0.3384; 0.7984]       1.6
## Fine                 -0.1800 [-0.4918; 0.1318]       5.3
## Grieger              -0.0600 [-0.3874; 0.2674]       4.8
## Rosenthal & Jacobson  0.3000 [ 0.0277; 0.5723]       6.9
## Fleming & Anttonen    0.0700 [-0.1139; 0.2539]      15.1
## Ginsburg             -0.0700 [-0.4112; 0.2712]       4.4
## 
## Number of studies combined: k = 19
## 
##                        SMD            95%-CI    z p-value
## Fixed effect model  0.0604 [-0.0111; 0.1318] 1.66  0.0979
## Prediction interval        [-0.2701; 0.4487]             
## 
## Quantifying heterogeneity:
## tau^2 = 0.0259; 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

Random effects model (DL method)

Perhaps unsurprisingly one would draw same the fundamental conclusion that there is no significant effect (though note this analysis the effect is estimated to be slightly larger .0893 vs. .0604).

Random_dl<-metagen(TE,
        seTE,
        data=Data,
        studlab=paste(author),
        comb.fixed = FALSE,
        comb.random = TRUE,
        hakn = FALSE,
        sm="SMD")
Random_dl
##                          SMD            95%-CI %W(random)
## Rosenthal et al.      0.0300 [-0.2148; 0.2748]        7.5
## Conn et al.           0.1200 [-0.1681; 0.4081]        6.6
## Jose & Cody          -0.1400 [-0.4674; 0.1874]        5.8
## Pellegrini & Hicks    1.1800 [ 0.4490; 1.9110]        1.9
## Pellegrini & Hicks    0.2600 [-0.4633; 0.9833]        1.9
## Evans & Rosenthal    -0.0600 [-0.2618; 0.1418]        8.5
## Fielder et al.       -0.0200 [-0.2218; 0.1818]        8.5
## Claiborn             -0.3200 [-0.7512; 0.1112]        4.2
## Kester                0.2700 [-0.0515; 0.5915]        5.9
## Maxwell               0.8000 [ 0.3081; 1.2919]        3.5
## Carter                0.5400 [-0.0519; 1.1319]        2.7
## Flowers               0.1800 [-0.2569; 0.6169]        4.1
## Keshock              -0.0200 [-0.5864; 0.5464]        2.8
## Henrikson             0.2300 [-0.3384; 0.7984]        2.8
## Fine                 -0.1800 [-0.4918; 0.1318]        6.1
## Grieger              -0.0600 [-0.3874; 0.2674]        5.8
## Rosenthal & Jacobson  0.3000 [ 0.0277; 0.5723]        6.9
## Fleming & Anttonen    0.0700 [-0.1139; 0.2539]        9.0
## Ginsburg             -0.0700 [-0.4112; 0.2712]        5.5
## 
## Number of studies combined: k = 19
## 
##                         SMD            95%-CI    z p-value
## Random effects model 0.0893 [-0.0200; 0.1987] 1.60  0.1094
## 
## Quantifying heterogeneity:
## tau^2 = 0.0259; 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
## - DerSimonian-Laird estimator for tau^2

Random effects model (Hartung-Knapp-Sidik-Jonkman method)

When this correction is applied we see that the confidence interval is somewhat widened, compared to when no correction is applied.

model_hksj<-metagen(TE,
        seTE,
        data=Data,
        studlab=paste(author),
        comb.fixed = FALSE,
        comb.random = TRUE,
        method.tau = "SJ",
        hakn = TRUE,
        prediction=TRUE,
        sm="SMD")
model_hksj
##                          SMD            95%-CI %W(random)
## Rosenthal et al.      0.0300 [-0.2148; 0.2748]        6.6
## Conn et al.           0.1200 [-0.1681; 0.4081]        6.2
## Jose & Cody          -0.1400 [-0.4674; 0.1874]        5.8
## Pellegrini & Hicks    1.1800 [ 0.4490; 1.9110]        2.8
## Pellegrini & Hicks    0.2600 [-0.4633; 0.9833]        2.9
## Evans & Rosenthal    -0.0600 [-0.2618; 0.1418]        6.9
## Fielder et al.       -0.0200 [-0.2218; 0.1818]        6.9
## Claiborn             -0.3200 [-0.7512; 0.1112]        4.9
## Kester                0.2700 [-0.0515; 0.5915]        5.9
## Maxwell               0.8000 [ 0.3081; 1.2919]        4.3
## Carter                0.5400 [-0.0519; 1.1319]        3.6
## Flowers               0.1800 [-0.2569; 0.6169]        4.8
## Keshock              -0.0200 [-0.5864; 0.5464]        3.8
## Henrikson             0.2300 [-0.3384; 0.7984]        3.8
## Fine                 -0.1800 [-0.4918; 0.1318]        5.9
## Grieger              -0.0600 [-0.3874; 0.2674]        5.8
## Rosenthal & Jacobson  0.3000 [ 0.0277; 0.5723]        6.3
## Fleming & Anttonen    0.0700 [-0.1139; 0.2539]        7.1
## Ginsburg             -0.0700 [-0.4112; 0.2712]        5.7
## 
## Number of studies combined: k = 19
## 
##                         SMD            95%-CI    t p-value
## Random effects model 0.1133 [-0.0355; 0.2622] 1.60  0.1270
## Prediction interval         [-0.4914; 0.7181]             
## 
## Quantifying heterogeneity:
## tau^2 = 0.0771; 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
## - Sidik-Jonkman estimator for tau^2
## - Hartung-Knapp adjustment for random effects model

REML estimation.

This is REML estimation which maps on Table 16.2 (Simple Random Effects model) (as an aside, note \(\tau^2\) estimate in this model as opposed to the above).

model_reml<-metagen(TE,
        seTE,
        data=Data,
        studlab=paste(author),
        comb.fixed = FALSE,
        comb.random = TRUE,
        method.tau = "REML",
        hakn = FALSE,
        prediction=TRUE,
        sm="SMD")
model_reml
##                          SMD            95%-CI %W(random)
## Rosenthal et al.      0.0300 [-0.2148; 0.2748]        7.7
## Conn et al.           0.1200 [-0.1681; 0.4081]        6.6
## Jose & Cody          -0.1400 [-0.4674; 0.1874]        5.7
## Pellegrini & Hicks    1.1800 [ 0.4490; 1.9110]        1.7
## Pellegrini & Hicks    0.2600 [-0.4633; 0.9833]        1.7
## Evans & Rosenthal    -0.0600 [-0.2618; 0.1418]        9.1
## Fielder et al.       -0.0200 [-0.2218; 0.1818]        9.1
## Claiborn             -0.3200 [-0.7512; 0.1112]        4.0
## Kester                0.2700 [-0.0515; 0.5915]        5.8
## Maxwell               0.8000 [ 0.3081; 1.2919]        3.3
## Carter                0.5400 [-0.0519; 1.1319]        2.4
## Flowers               0.1800 [-0.2569; 0.6169]        3.9
## Keshock              -0.0200 [-0.5864; 0.5464]        2.6
## Henrikson             0.2300 [-0.3384; 0.7984]        2.6
## Fine                 -0.1800 [-0.4918; 0.1318]        6.0
## Grieger              -0.0600 [-0.3874; 0.2674]        5.7
## Rosenthal & Jacobson  0.3000 [ 0.0277; 0.5723]        7.0
## Fleming & Anttonen    0.0700 [-0.1139; 0.2539]        9.7
## Ginsburg             -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

Forest plot.

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.0     ✔ readr   1.3.1
## ✔ tibble  2.1.3     ✔ purrr   0.3.2
## ✔ tidyr   0.8.3     ✔ stringr 1.4.0
## ✔ ggplot2 3.2.0     ✔ forcats 0.4.0
## ── Conflicts ─────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
Data <-Data %>% unite(author_year, c(author, year_b), sep = " ", remove = FALSE)

Let’s redo our model but now with author_year.

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

Note that this is but one solution.

Unfortunately this doesn’t fit nicely, unless you mess about with the chunk options in knitr. So let’s also export it as .pdf . Have a look at the meta manual manual, you can modify a bunch of other things.

forest(model_reml2,
       sortvar=TE,
       xlim = c(-1,2),
       rightlabs = c("SMD","95% CI","weight"),
       leftlabs = c("Author", "SMD", "SE"),
       pooled.totals = FALSE,
       smlab = "",
       text.random = "Overall effect",
       print.tau2 = FALSE,
       col.diamond = "darkblue",
       col.diamond.lines = "black",
       col.predict = "black",
       print.I2.ci = TRUE,
       digits.sd = 2
)

pdf(file='forestplot_pygmalion.pdf',width=10,height=8) 
forest(model_reml,
       sortvar=TE,
       xlim = c(-1,2),
       rightlabs = c("SMD","95% CI","weight"),
       leftlabs = c("Author", "N","Mean","SD","N","Mean","SD"),
       lab.e = "Intervention",
       pooled.totals = FALSE,
       smlab = "",
       text.random = "Overall effect",
       print.tau2 = FALSE,
       col.diamond = "darkblue",
       col.diamond.lines = "black",
       col.predict = "black",
       print.I2.ci = TRUE,
       digits.sd = 2
)

Note that there is some unnecessary duplication in the plot but it seems that ‘leftlabs’ cannot be modified easily… .

Here is an example of using ‘JAMA’.

pdf(file='JAMA_forestplot_pygmalion.pdf', width=10,height=8) 
forest(model_reml2,
       sortvar=TE,
      layout='JAMA'
)

You should get something like below.