############################################################################ # R script to go along with the talk: # # Viechtbauer, W. (2022). Meta-analysis, open science, and publication bias. # Open Science Workshop: Promoting transparency & replicability in research. # University of Tübingen. ############################################################################ # meta-analysis on the effectiveness of magnesium treatment in acute # myocardial infarction for reducing the risk of mortality # References: # # Teo, K. K., Yusuf, S., Collins, R., Held, P. H., & Peto, R. (1991). Effects # of intravenous magnesium in suspected acute myocardial infarction: Overview # of randomised trials. British Medical Journal, 303(6816), 1499-1503. # https://doi.org/10.1136/bmj.303.6816.1499 # # Horner, S. M. (1992). Efficacy of intravenous magnesium in acute myocardial # infarction in reducing arrhythmias and mortality: Meta-analysis of magnesium # in acute myocardial infarction. Circulation, 86(3), 774-779. # https://doi.org/10.1161/01.cir.86.3.774 # # Yusuf, S., Teo, K. & Woods, K. (1993). Intravenous magnesium in acute # myocardial infarction: An effective, safe, simple, and inexpensive # treatment. Circulation, 87(6), 2043-2046. # https://doi.org/10.1161/01.cir.87.6.2043 # # Egger, M., Davey Smith, G., & Altman, D. G. (Eds.) (2001). Systematic # reviews in health care: Meta-analysis in context (2nd ed.). BMJ Books. ############################################################################ # install the metafor package install.packages("metafor") # load the metafor package library(metafor) ############################################################################ # for a description of the dataset, see: help(dat.egger2001) # copy magnesium treatment dataset to 'dat' and examine the dataset dat <- dat.egger2001 dat # remove studies 8 and 16 dat <- dat[-c(8,16),] # study 16 is the ISIS-4 trial, which wasn't available for the meta-analysis; # and study 8 is a small study with a very low number of events and hence a # very large sampling variance; its inclusion would make this example a bit # less useful didactically, so we'll also exclude it # data are given in terms of 2x2 tables of the form: # # | dead alive | total # ----------+----------------+------- # magnesium | ai n1i-ai | n1i # control | ci n2i-ci | n2i # ----------+----------------+------- # compute the log risk ratios (using +1/2 adjustment for all studies) dat <- escalc(measure="RR", ai=ai, n1i=n1i, ci=ci, n2i=n2i, data=dat, add=1/2, to="all", slab=paste0(study, ", ", year)) dat # variable 'yi' contains the log risk ratios and 'vi' the corresponding # sampling variances; a negative value for yi indicates a lower risk of # mortality in the magnesium group # fit random-effects model res <- rma(yi, vi, data=dat) res # forest plot forest(res, header=TRUE) # estimated average risk ratio (with 95% confidence and prediction interval) predict(res, transf=exp, digits=2) # these results suggest that magnesium treatment leads on average to a 50% # reduction in mortality risk; however, this contradicts the results from the # ISIS-4 study, which found no reduction in risk in the magnesium group # can also obtain a short 'write up' of the analysis with reporter() reporter(res) ############################################################################ # funnel plot funnel(res, ylim=c(0,1.2)) # failsafe-N based on Rosenthal method (examines how many studies it would # take to change the observed p-value to .05; the calculation is based on a # method for combining p-values which is not used in modern meta-analyses) fsn(yi, vi, data=dat) # failsafe-N based on the 'Orwin method' (examines how many studies with null # effects it would take to reduce the average effect size to some target # value; by default, to half its size) fsn(yi, vi, data=dat, type="Orwin") # failsafe-N based on the 'Rosenberg method' (examines how many studies with # null effects it would take to reduce the estimate from an equal-effects # model to non-significance) fsn(yi, vi, data=dat, type="Rosenberg") # regression test for funnel plot asymmetry # meta-regression model with sqrt(vi) as predictor dat$sei <- sqrt(dat$vi) res <- rma(yi, vi, mods = ~ sei, data=dat) res # or use the regtest() function res <- rma(yi, vi, data=dat) regtest(res) # test of excess significance tes(yi, vi, data=dat) # trim and fill method res <- rma(yi, vi, data=dat) taf <- trimfill(res) taf # funnel plot with filled-in studies funnel(taf) # fit selection model (based on an equal-effects model and assuming selection # in favor of significant negative effects) res <- rma(yi, vi, method="EE", data=dat) sel <- selmodel(res, type="logistic", alternative="less") sel plot(sel) # fitting a random-effects selection model also works, but the estimate of # tau^2 drifts towards 0, causing some problems with getting the standard # errors of the estimates; the results are essentially the same though as # those for the equal-effects model res <- rma(yi, vi, data=dat) sel <- selmodel(res, type="logistic", alternative="less") sel # can get the standard errors also for the random-effects selection model when # switching to another way of computing/inverting the 'Hessian' matrix sel <- selmodel(res, type="logistic", alternative="less", control=list(hesspack="pracma")) sel # PET/PEESE methods rma(yi, vi, mods = ~ sei, data=dat) rma(yi, vi, mods = ~ vi, data=dat) # or use regtest() (see the 'limit estimate') res <- rma(yi, vi, data=dat) regtest(res) regtest(res, predictor="vi") # show funnel plot with PET/PEESE results sav <- funnel(res, ylim=c(0,1.2)) reg <- regtest(res) se <- seq(0, 1.3, length=100) lines(coef(reg$fit)[1] + coef(reg$fit)[2]*se, se, lwd=5, col="dodgerblue") reg <- regtest(res, predictor="vi") lines(coef(reg$fit)[1] + coef(reg$fit)[2]*se^2, se, lwd=5, col="firebrick") legend("topright", inset=.02, legend=c("PET","PEESE"), col=c("dodgerblue","firebrick"), lwd=5, bg="white") points(sav$x, sav$y, pch=19) ############################################################################