############################################################### #### A Practical Tutorial on Conducting Meta-Analysis in R #### #### by A. C. Del Re #### ############################################################### #### CITATION: #### Del Re, A. C. (2015). A Practical Tutorial on Conducting #### Meta-Analysis in R. The Quantitative Methods for #### Psychology, 11 (1), 37-50. #### THIS FILE WILL GENERATE THE DATASETS THAT ARE USED AS A #### RUNNING EXAMPLE IN THE TUTORIAL. PART TWO OF THIS FILE #### (tutorial_ma_models.R) PROVIDES CODE TO REPLICATE THE #### MODELS/ANALYSES IN THE TUTORIAL #### LIBRARIES #### # # INSTALL PACKAGES (IF NOT ALREADY IN THE LIBRARY) # install.packages("mvtnorm") # install.packages("compute.es") # install.packages("MAd") # install.packages("metafor") # install.packages("ggplot2") # install.packages("gridExtra") # install.packages("Gmisc") # LOAD THE PACKAGES INTO THE CURRENT R SESSION library(mvtnorm) # TO SIMULATE MULTIVARIATE DATA library(compute.es) # TO COMPUTE EFFECT SIZES library(MAd) # META-ANALYSIS PACKAGE library(metafor) # META-ANALYSIS PACKAGE library(ggplot2) # GRAPHICS PACKAGE library(gridExtra) # ARRANGE GRID FOR GRAPHICS library(Gmisc) # TABLE GENERATION library(scales) # ADDITIONAL GRAPHICS #### SIMULATE FICTIONAL PSYCHOTHERAPY DATASET #### set.seed(323232) # TO REPLICATE DATA # POPULATION MEANS AND SDs: mu1 <- 17 # CONTROL GRP POST-TREATMENT OUTCOME 1 MEAN mu2 <- 8 # CONTROL GRP POST-TREATMENT OUTCOME 2 MEAN mu3 <- 20 # TREATMENT GRP POST-TREATMENT OUTCOME 1 MEAN mu4 <- 10 # TREATMENT GRP POST-TREATMENT OUTCOME 2 MEAN sig1 <- 6 # T AND C GRP POST-TREATMENT OUTCOME 1 SD sig2 <- 2.45 # T AND C GRP POST-TREATMENT OUTCOME 2 SD n <- 30 # AVERAGE SAMPLE SIZE FOR T AND C GRPS t_dose <- 10 # AVERAGE TREATMENT DURATION ("DOSAGE") distress <- 65 # AVERAGE PATIENT DISTRESS (SCALE FROM 1 TO 100; 1=NO DISTRESS) t_dose.sig <- 9 # AVERAGE TREATMENT DURATION SD distress.sig <- 25 # AVERAGE PATIENT DISTRESS SD nsamp <- 1000 # NUMBER OF REPLICATIONS (1000 STUDIES) # TRUE POPULATION CORR BETWEEN OUTCOMES AT r=.5, # TRUE POPULATION CORR BETWEEN MODERATORS AT r=.6 # TRUE POPULATION CORR BETWEEN MODERATORS AND OUTCOMES: # r(O1,M1)=.8; r(O2,M1)=.8; r(O1,M2)=.34; r(O2,M2)=.34 # O1, O2, M1, M2 a <-matrix(c(1, .5, .8, .34, # O1 .5, 1, .8, .34, # O2 .8, .8, 1, .6, # M1 .34, .34, .6, 1), # M2 4, byrow=TRUE) stdevs <- c(sig1,sig2, t_dose.sig, distress.sig) # SD OF VARIABLES b <- stdevs %*% t(stdevs) # N*N MATRIX WITH stdev[i]*stdev[j] sigma <- b * a # COVAR MATRIX # out IS A DATAFRAME TO STORE RESULTS, ONE ROW PER SAMPLE out <- data.frame(id=1:1000, nT = rep(NA,nsamp), m1T = rep(NA,nsamp), sd1T = rep(NA,nsamp), m2T = rep(NA,nsamp), sd2T = rep(NA,nsamp), nC = rep(NA,nsamp), m1C = rep(NA,nsamp), sd1C = rep(NA,nsamp), m2C = rep(NA,nsamp), sd2C = rep(NA,nsamp), dose = rep(NA,nsamp), stress = rep(NA,nsamp) ) for(iii in 1:nsamp) { # RANDOM SAMPLES NC <- round(abs(rnorm(1, mean = n, sd = 7)), 0) # C GRP N NT <- NC # T GRP N (EQUAL SAMPLE SIZES FOR T AND C GRPS) y12 <- rmvnorm(n = n, mean = c(mu1,mu2, t_dose, distress), sigma = sigma) # C GRP OUTCOMES y34 <- rmvnorm(n = n, mean = c(mu3,mu4, t_dose, distress), sigma = sigma) # T GRP OUTCOMES out$nT[iii] <- round(mean(NT), 0) out$nC[iii] <- round(mean(NC), 0) out$m1T[iii] <- mean(y34[,1]) out$m1C[iii] <- mean(y12[,1]) out$m2T[iii] <- mean(y34[,2]) out$m2C[iii] <- mean(y12[,2]) out$sd1T[iii] <- sd(y34[,1]) out$sd1C[iii] <- sd(y12[,1]) out$sd2T[iii] <- sd(y34[,2]) out$sd2C[iii] <- sd(y12[,2]) out$dose[iii] <- round(mean( y34[,3]), 0) out$stress[iii] <- round(mean( y34[,4]), 0) } # ALL STUDIES OUTCOME 1 dat.sim1 <- mes(m.1=m1T, m.2=m1C, sd.1=sd1T, sd.2=sd1C, n.1=nT, n.2=nC, level=95, dig=2, id=id, data=out)[,c(1,3,4,13,14,19)] # ALL STUDIES OUTCOME 2 dat.sim2 <- mes(m.1=m2T, m.2=m2C, sd.1=sd2T, sd.2=sd2C, n.1=nT, n.2=nC, level=95, dig=2, id=id, data=out)[,c(1,3,4,13,14,19)] # DICHOTOMIZE STRESS INTO HIGH AND LOW DISTRESS GROUPS FOR DEMO PURPOSE out$stress <- with(out, ifelse(stress<=(median(stress)-5), "low", 'high')) out$stress <- as.factor(out$stress) # TAKE A RANDOM SAMPLE OF 8 STUDIES FROM THE "UNIVERSE" OF STUDIES (K=1000 HERE) set.seed(932329) # TO REPLICATE DATA dat.sim.raw <- out[sample(out$id, 8),] row.names(dat.sim.raw)<-NULL #### COMPUTE EFFECT SIZES #### # USING THE compute.es PACKAGE # OUTCOME 1 EFFECT SIZES res1 <- mes(m.1=m1T, m.2=m1C, sd.1=sd1T, sd.2=sd1C, n.1=nT, n.2=nC, level=95, dig=2, id=id, data=dat.sim.raw)[,c(1,3,4,13,14,19)] names(res1) <- c("id", "nT", "nC", "g", "var.g", "pval.g") res1$outcome <- "g1" # ADD COLUMN FOR OUTCOME # OUTCOME 2 EFFECT SIZES res2 <- mes(m.1=m2T, m.2=m2C, sd.1=sd2T, sd.2=sd2C, n.1=nT, n.2=nC, level=95, dig=2, id=id, data=dat.sim.raw)[,c(1,3,4,13,14,19)] names(res2) <- c("id", "nT", "nC", "g", "var.g", "pval.g") res2$outcome <- "g2" # ADD COLUMN FOR OUTCOME dat.sim.es <-rbind(res1, res2) # BIND EACH DATA SET TOGETHER #### R SESSION INFORMATION #### sessionInfo() # R version 3.1.1 (2014-07-10) # Platform: x86_64-w64-mingw32/x64 (64-bit) # # locale: # [1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252 # [3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C # [5] LC_TIME=English_United States.1252 # # attached base packages: # [1] splines grid stats graphics grDevices utils datasets methods base # # other attached packages: # [1] quantreg_5.05 SparseM_1.05 scales_0.2.4 Gmisc_0.6.7 Hmisc_3.14-5 # [6] survival_2.37-7 lattice_0.20-29 gridExtra_0.9.1 ggplot2_1.0.0 metafor_1.9-4 # [11] Matrix_1.1-4 Formula_1.1-2 MAd_0.8-2 compute.es_0.2-4 mvtnorm_1.0-0 # # loaded via a namespace (and not attached): # [1] acepack_1.3-3.3 cluster_1.15.2 colorspace_1.2-4 digest_0.6.4 foreign_0.8-61 # [6] gtable_0.1.2 labeling_0.3 latticeExtra_0.6-26 MASS_7.3-33 munsell_0.4.2 # [11] nnet_7.3-8 plyr_1.8.1 proto_0.3-10 RColorBrewer_1.0-5 Rcpp_0.11.3 # [16] reshape2_1.4 rpart_4.1-8 sp_1.0-15 stringr_0.6.2 tools_3.1.1