Solution Session 4
Dr. Thomas Pollet, Northumbria University (thomas.pollet@northumbria.ac.uk)
2020-03-16 | disclaimer
Questions.
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
yiin this dataset, as a Cohen’s d.
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, trimfillData<-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.0106I 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, unionData <- 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^2Funnel 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.1790311Begg 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.5482Duval & 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 asymmetryContour 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)
Intermezzo: Why no fail-safe N’s?
The effect wasn’t significant to begin with. So, in order to reduce its significance we would have added 0 studies.
p-uniform
This is not optimal, but it looks like for our purpose we’ll have to convert our treatment effect (yi or Te) to a metric which we can feed to p-uniform. If you were to have the M and SD for each condition, you could use those. However, in this case, they are not reported in our dataset. From the materials we read that the effect sizes are standard mean differences (‘SMD’), though it has not been noted as to whether these are Cohen’s d or not. For our purpose we will treat them as Cohen’s d . So we are going to have to do some data-wrangling
Remember: \[r=\frac{d}{\sqrt{d^2+A}}\] and A is our correction factor.
\[A= \frac{(n_1+n_2)^2}{n_1n_2}\]
Correction factor A
So let"s use dplyrto calculate our correction factor A.
require(dplyr)
Data <- Data %>% mutate(correction_factor_A= ((n1i+n2i)^2)/(n1i*n2i))Have a look at the dataframe, you should notice that studies where treatment and control sample sizes were the same have A=4.
Pearson correlation from SMD estimate
Data <- Data %>% mutate(pearson_r= TE/(sqrt((TE^2)+correction_factor_A)))Finally we’ll also need to weigh those correlation coefficients, here we simple sum n1i and n2i to obtain N (here labelled pearson_weight).
Data <- Data %>% mutate(pearson_weight= (n1i+n2i))P-uniform
require(puniform)## Loading required package: puniform## Warning: package 'puniform' was built under R version 3.5.2## Loading 'puniform' package (Version 0.1.1)puni_star(ri = Data$pearson_r, ni=Data$pearson_weight, side='right')##
## Method: ML (k = 19; ksig = 3)
##
## Estimating effect size p-uniform*
##
## est ci.lb ci.ub L.0 pval
## 0.0105 NA NA 0.3262 0.5679
##
## ===
##
## Estimating between-study variance p-uniform*
##
## tau2 tau2.lb tau2.ub L.het pval
## 0 0 0.0068 0 1
##
## ===
##
## Publication bias test p-uniform*
##
## L.pb pval
## 3.5264 0.1715Perhaps unsurprisingly as there are only 3 significant estimates, the p uniform method does not suggest issues.
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 7 on publication bias 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.
Session info.
sessionInfo()## R version 3.5.1 (2018-07-02)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Sierra 10.12.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] puniform_0.1.1 forcats_0.4.0 stringr_1.4.0 purrr_0.3.3
## [5] readr_1.3.1 tidyr_1.0.2 tibble_2.1.3 ggplot2_3.2.1
## [9] tidyverse_1.2.1 dplyr_0.8.4 metafor_2.1-0 Matrix_1.2-14
## [13] meta_4.9-6
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.0.0 xfun_0.12 haven_2.1.1 lattice_0.20-35
## [5] colorspace_1.4-1 vctrs_0.2.3 generics_0.0.2 htmltools_0.4.0
## [9] yaml_2.2.1 rlang_0.4.4 pillar_1.4.3 withr_2.1.2
## [13] glue_1.3.1 modelr_0.1.2 readxl_1.3.1 lifecycle_0.1.0
## [17] munsell_0.5.0 gtable_0.3.0 cellranger_1.1.0 rvest_0.3.2
## [21] evaluate_0.14 knitr_1.28 fansi_0.4.1 broom_0.5.1
## [25] Rcpp_1.0.3 backports_1.1.5 scales_1.1.0 jsonlite_1.6.1
## [29] hms_0.5.1 digest_0.6.25 stringi_1.4.6 grid_3.5.1
## [33] cli_2.0.1 tools_3.5.1 magrittr_1.5 lazyeval_0.2.2
## [37] crayon_1.3.4 pkgconfig_2.0.3 ellipsis_0.3.0 xml2_1.2.2
## [41] lubridate_1.7.4 assertthat_0.2.1 rmarkdown_1.11 httr_1.4.0
## [45] rstudioapi_0.11 R6_2.4.1 nlme_3.1-137 compiler_3.5.1