Questions.

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

Load and manipulate the data

There are 19 studies.

library(meta)
library(metafor) # needed later on.
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)
## 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)
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

Subgroup analyses.

Below we perform the models by group.

Interestingly, it seems that the effect is stronger in the blinded group than in the aware group - perhaps they were better studies? However, there is no significant difference between the groups ( Q=.84, p =.367).

tester_subgroup_common<-update.meta(model_reml2, 
                             byvar=tester, 
                             comb.random = TRUE, 
                             comb.fixed = FALSE,
                             tau.common=TRUE)
tester_subgroup_common
##                                 SMD            95%-CI %W(random) tester
## Rosenthal et al. (1974)      0.0300 [-0.2148; 0.2748]        7.7  aware
## Conn et al. (1968)           0.1200 [-0.1681; 0.4081]        6.6  aware
## Jose & Cody (1971)          -0.1400 [-0.4674; 0.1874]        5.7  aware
## Pellegrini & Hicks (1972)    1.1800 [ 0.4490; 1.9110]        1.7  aware
## Pellegrini & Hicks (1972)    0.2600 [-0.4633; 0.9833]        1.7  blind
## Evans & Rosenthal (1969)    -0.0600 [-0.2618; 0.1418]        9.1  aware
## Fielder et al. (1971)       -0.0200 [-0.2218; 0.1818]        9.1  blind
## Claiborn (1969)             -0.3200 [-0.7512; 0.1112]        4.0  aware
## Kester (1969)                0.2700 [-0.0515; 0.5915]        5.8  aware
## Maxwell (1970)               0.8000 [ 0.3081; 1.2919]        3.3  blind
## Carter (1970)                0.5400 [-0.0519; 1.1319]        2.4  blind
## Flowers (1966)               0.1800 [-0.2569; 0.6169]        3.9  blind
## Keshock (1970)              -0.0200 [-0.5864; 0.5464]        2.6  blind
## Henrikson (1970)             0.2300 [-0.3384; 0.7984]        2.6  blind
## Fine (1972)                 -0.1800 [-0.4918; 0.1318]        6.0  aware
## Grieger (1970)              -0.0600 [-0.3874; 0.2674]        5.7  blind
## Rosenthal & Jacobson (1968)  0.3000 [ 0.0277; 0.5723]        7.0  aware
## Fleming & Anttonen (1971)    0.0700 [-0.1139; 0.2539]        9.7  blind
## Ginsburg (1970)             -0.0700 [-0.4112; 0.2712]        5.4  aware
## 
## 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%]
## 
## Quantifying residual heterogeneity:
## tau^2 = 0.0248; H = 1.38; I^2 = 47.5%
## 
## Test of heterogeneity:
##      Q d.f. p-value
##  35.83   18  0.0074
## 
## Results for subgroups (random effects model):
##                  k    SMD            95%-CI     Q  tau^2   I^2
## tester = aware  10 0.0463 [-0.0945; 0.1872] 22.20 0.0248 59.5%
## tester = blind   9 0.1488 [-0.0198; 0.3174] 12.96 0.0248 38.3%
## 
## Test for subgroup differences (random effects model):
##                    Q d.f. p-value
## Between groups  0.84    1  0.3606
## Within groups  35.16   17  0.0059
## 
## Details on meta-analytical method:
## - Inverse variance method
## - Restricted maximum-likelihood estimator for tau^2 (assuming common tau^2 in subgroups)

This model allows an estimate of variance \(\tau^2\) for each group. The conclusion is the same there is no evidence for a significant difference between these groups ( Q=.81, p =.367).

tester_subgroup_sep<-update.meta(model_reml2, 
                             byvar=tester, 
                             comb.random = TRUE, 
                             comb.fixed = FALSE,
                             tau.common=FALSE)
tester_subgroup_sep
##                                 SMD            95%-CI %W(random) tester
## Rosenthal et al. (1974)      0.0300 [-0.2148; 0.2748]        7.7  aware
## Conn et al. (1968)           0.1200 [-0.1681; 0.4081]        6.6  aware
## Jose & Cody (1971)          -0.1400 [-0.4674; 0.1874]        5.7  aware
## Pellegrini & Hicks (1972)    1.1800 [ 0.4490; 1.9110]        1.7  aware
## Pellegrini & Hicks (1972)    0.2600 [-0.4633; 0.9833]        1.7  blind
## Evans & Rosenthal (1969)    -0.0600 [-0.2618; 0.1418]        9.1  aware
## Fielder et al. (1971)       -0.0200 [-0.2218; 0.1818]        9.1  blind
## Claiborn (1969)             -0.3200 [-0.7512; 0.1112]        4.0  aware
## Kester (1969)                0.2700 [-0.0515; 0.5915]        5.8  aware
## Maxwell (1970)               0.8000 [ 0.3081; 1.2919]        3.3  blind
## Carter (1970)                0.5400 [-0.0519; 1.1319]        2.4  blind
## Flowers (1966)               0.1800 [-0.2569; 0.6169]        3.9  blind
## Keshock (1970)              -0.0200 [-0.5864; 0.5464]        2.6  blind
## Henrikson (1970)             0.2300 [-0.3384; 0.7984]        2.6  blind
## Fine (1972)                 -0.1800 [-0.4918; 0.1318]        6.0  aware
## Grieger (1970)              -0.0600 [-0.3874; 0.2674]        5.7  blind
## Rosenthal & Jacobson (1968)  0.3000 [ 0.0277; 0.5723]        7.0  aware
## Fleming & Anttonen (1971)    0.0700 [-0.1139; 0.2539]        9.7  blind
## Ginsburg (1970)             -0.0700 [-0.4112; 0.2712]        5.4  aware
## 
## 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%]
## 
## Quantifying residual heterogeneity:
## H = 1.44 [1.10; 1.88]; I^2 = 51.6% [17.1%; 71.8%]
## 
## Test of heterogeneity:
##      Q d.f. p-value
##  35.83   18  0.0074
## 
## Results for subgroups (random effects model):
##                  k    SMD            95%-CI     Q  tau^2   I^2
## tester = aware  10 0.0466 [-0.0954; 0.1887] 22.20 0.0256 59.5%
## tester = blind   9 0.1474 [-0.0193; 0.3140] 12.96 0.0236 38.3%
## 
## Test for subgroup differences (random effects model):
##                     Q d.f. p-value
## Between groups   0.81    1  0.3672
## 
## Details on meta-analytical method:
## - Inverse variance method
## - Restricted maximum-likelihood estimator for tau^2

Meta-regression with year

Center year

Data <- Data %>% mutate(year_cent = year-mean(year))

Meta-regression

We need to rerun our model as we have added our new moderator variable, note that this overrides our previous model.

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

There is no suggestion that publication year is a viable predictor of the effect size (Q = .41, p= .52).

metareg_pub_year<-metareg(model_reml2,year_cent)
metareg_pub_year
## 
## Mixed-Effects Model (k = 19; tau^2 estimator: REML)
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0201 (SE = 0.0167)
## tau (square root of estimated tau^2 value):             0.1417
## I^2 (residual heterogeneity / unaccounted variability): 42.84%
## H^2 (unaccounted variability / sampling variability):   1.75
## R^2 (amount of heterogeneity accounted for):            0.00%
## 
## Test for Residual Heterogeneity:
## QE(df = 17) = 34.7475, p-val = 0.0067
## 
## Test of Moderators (coefficient 2):
## QM(df = 1) = 0.4081, p-val = 0.5229
## 
## Model Results:
## 
##            estimate      se     zval    pval    ci.lb   ci.ub 
## intrcpt      0.0851  0.0524   1.6239  0.1044  -0.0176  0.1879    
## year_cent   -0.0187  0.0293  -0.6388  0.5229  -0.0761  0.0387    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Plot

bubble.metareg(metareg_pub_year,
              xlab = "Publication Year (centered)",
              ylab = "SMD",
              col.line = "hotpink",
              studlab = TRUE)

Permutation test

We’ll use set.seed() to ensure we get the same results every single time.

The permutation test corroborates our findings from the meta-regression. There is no substantial evidence for an effect of publication year.

set.seed(1981)
permutest(metareg_pub_year, iter=5000, progbar=F)
## 
## Test of Moderators (coefficient 2):
## QM(df = 1) = 0.4081, p-val* = 0.5482
## 
## Model Results:
## 
##            estimate      se     zval   pval*    ci.lb   ci.ub 
## intrcpt      0.0851  0.0524   1.6239  0.4716  -0.0176  0.1879    
## year_cent   -0.0187  0.0293  -0.6388  0.5482  -0.0761  0.0387    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Compare models

Here we make a switch to metafor.

As can be expected this corroborates what we previously found when examining subgroups. Note that I have opted for ML estimation with the Knapp-Hartung estimation. As model comparisons aren’t well defined with REML methods, we switched to ML estimation in order to allow

metareg_tester<-rma(yi=TE, 
              sei=seTE, 
              data=Data, 
              method = "ML", 
              mods = ~ tester, 
              test="knha")
metareg_tester
## 
## Mixed-Effects Model (k = 19; tau^2 estimator: ML)
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0133 (SE = 0.0126)
## tau (square root of estimated tau^2 value):             0.1155
## I^2 (residual heterogeneity / unaccounted variability): 32.78%
## H^2 (unaccounted variability / sampling variability):   1.49
## R^2 (amount of heterogeneity accounted for):            0.00%
## 
## Test for Residual Heterogeneity:
## QE(df = 17) = 35.1572, p-val = 0.0059
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 17) = 0.5293, p-val = 0.4768
## 
## Model Results:
## 
##              estimate      se    tval    pval    ci.lb   ci.ub 
## intrcpt        0.0420  0.0786  0.5348  0.5997  -0.1237  0.2077    
## testerblind    0.0894  0.1229  0.7275  0.4768  -0.1699  0.3487    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
metareg_pub_tester<-rma(yi=TE, 
              sei=seTE, 
              data=Data, 
              method = "ML", 
              mods = ~ tester + year_cent, 
              test="knha")
metareg_pub_tester
## 
## Mixed-Effects Model (k = 19; tau^2 estimator: ML)
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0071 (SE = 0.0099)
## tau (square root of estimated tau^2 value):             0.0842
## I^2 (residual heterogeneity / unaccounted variability): 20.17%
## H^2 (unaccounted variability / sampling variability):   1.25
## R^2 (amount of heterogeneity accounted for):            43.53%
## 
## Test for Residual Heterogeneity:
## QE(df = 16) = 33.8848, p-val = 0.0056
## 
## Test of Moderators (coefficients 2:3):
## F(df1 = 2, df2 = 16) = 0.4516, p-val = 0.6445
## 
## Model Results:
## 
##              estimate      se     tval    pval    ci.lb   ci.ub 
## intrcpt        0.0390  0.0762   0.5122  0.6155  -0.1225  0.2006    
## testerblind    0.0826  0.1190   0.6943  0.4975  -0.1697  0.3349    
## year_cent     -0.0225  0.0329  -0.6844  0.5035  -0.0922  0.0472    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(metareg_tester,metareg_pub_tester)
## 
##         df     AIC     BIC    AICc  logLik    LRT   pval      QE  tau^2 
## Full     4 12.6903 16.4681 15.5475 -2.3452               33.8848 0.0071 
## Reduced  3 11.4095 14.2428 13.0095 -2.7048 0.7192 0.3964 35.1572 0.0133 
##              R^2 
## Full 
## Reduced 46.7824%

Both the AIC/BIC (as well as the Likelihood ratio test (LRT)) do not suggest superiority of one model over another.

Interaction model

We reverted back to the REML model estimation. This model does not support an interaction effect between publication year and tester blindness. Even if it had, we should be very wary as we only have 19 studies in this meta-analysis.

metareg_interaction<-rma(yi=TE, 
              sei=seTE, 
              data=Data, 
              method = "REML", 
              mods = ~ tester * year_cent 
              )
metareg_interaction
## 
## Mixed-Effects Model (k = 19; tau^2 estimator: REML)
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0338 (SE = 0.0242)
## tau (square root of estimated tau^2 value):             0.1840
## I^2 (residual heterogeneity / unaccounted variability): 54.60%
## H^2 (unaccounted variability / sampling variability):   2.20
## R^2 (amount of heterogeneity accounted for):            0.00%
## 
## Test for Residual Heterogeneity:
## QE(df = 15) = 33.7799, p-val = 0.0037
## 
## Test of Moderators (coefficients 2:4):
## QM(df = 3) = 1.0372, p-val = 0.7922
## 
## Model Results:
## 
##                        estimate      se     zval    pval    ci.lb   ci.ub 
## intrcpt                  0.0495  0.0784   0.6310  0.5281  -0.1042  0.2031 
## testerblind              0.1084  0.1217   0.8907  0.3731  -0.1302  0.3470 
## year_cent               -0.0127  0.0399  -0.3186  0.7501  -0.0910  0.0655 
## testerblind:year_cent   -0.0101  0.0737  -0.1364  0.8915  -0.1546  0.1345 
##  
## intrcpt 
## testerblind 
## year_cent 
## testerblind:year_cent 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Acknowledgments and further reading… .

The example is from here.

Note that throughout I have varied the rounding when I reported, you should make your own decisions on how precise you believe your results to be.

Please see the slides for further reading but a good place to start is Chapter 8 on models and approaches to inference in Koricheva, J., Gurevitch, J., & Mengersen, K. (2013). Handbook of Meta-Analysis in Ecology and Evolution. Princeton, NJ: Princeton University Press.

Cited literature

Raudenbush, S. W. (1984). Magnitude of teacher expectancy effects on pupil IQ as a function of the credibility of expectancy induction: A synthesis of findings from 18 experiments. Journal of Educational Psychology, 76(1), 85–97.

Raudenbush, S. W. (2009). Analyzing effect sizes: Random effects models. In H. Cooper, L. V. Hedges, & J. C. Valentine (Eds.), The handbook of research synthesis and meta-analysis (2nd ed., pp. 295–315). New York: Russell Sage Foundation.

The end.