###############################################################
#### 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 PROVIDES CODE TO REPLICATE THE
#### MODELS/ANALYSES IN THE TUTORIAL. PART ONE OF THIS FILE
#### (tutorial_ma_data.R) WILL GENERATE THE DATASETS THAT 
#### ARE USED AS A RUNNING EXAMPLE 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


#### SOURCE PART ONE OF THE TUTORIAL FILES ####

# source('C:/Users/acd/tutorial_ma_data.R')

# NOTE: CHANGE PATH (E.G., 'YOUR_COMPUTER_PATH' BELOW)  TO THE LOCATION OF THE 
# tutorial_ma_data.R FILE (BE SURE TO USE '/' INSTEAD OF '\')
# E.G.,
# source('YOUR_COMPUTER_PATH'/tutorial_ma_data.R', echo=TRUE)
# AND REMOVE '#' BEFORE source()

#### AGGREGATE DEPENDENT EFFECT SIZES ####

dat.sim.agg <- agg(id  =  id, es  =  g, var  =  var.g, n.1  =  nT, n.2  = 
                   nC, cor  = 0.5,  method  =  "BHHR", data  =  dat.sim.es)

# ADD MODERATORS BACK INTO THE DATASET
dat.sim.final <- merge(dat.sim.agg, dat.sim.raw[, c(1,12:13)], by='id')


#### RANDOM EFFECTS MODELS ####

# OMNIBUS
m0 <- mareg(es~1, var=var, data=dat.sim.final)
summary(m0)

# META-REGRESSION
m1 <- mareg(es~dose, var=var, data=dat.sim.final)
summary(m1)  # SINGLE MODERATOR

m2 <- mareg(es~stress, var=var, data=dat.sim.final)
summary(m2)  # SINGLE MODERATOR

m3 <- mareg(es~dose + stress, var=var, data=dat.sim.final)
summary(m3)  # MULTIPLE MODERATOR

## CENTER THE INTERCEPT FOR BETTER INTERPRETATION

# HERE INTERCEPT IS CENTERED AT THE AVERAGE NUMBER OF SESSIONS
summary(mareg(es~I(dose-mean(dose)), var=var, data=dat.sim.final))

# HERE INTERCEPT IS CENTERED FOR DOSE IN THE MULTIPLE MODERATOR ANALYSIS
summary(mareg(es~I(dose-mean(dose))+stress, var=var, data=dat.sim.final))

# CORRELATED MODERATORS--CONFOUNDS? EXAMINE CORRELATION BETWEEN MODERATORS
# (IN THIS CASE THE CATEGORICAL MOD TAKES ON ONLY 2 LEVELS, AND CONVERTING TO 
# A CONTINUOUS VARIABLE IS APPROPRIATE)
with(dat.sim.final, cor(dose, as.numeric(stress))) # r=.58 (LARGE CORRELATION)


confint(m3)

#' MODEL COMPARISONS
# likelihood ratio test 
# alternative to Wald test, this conducts a full vs reduced model
# comparison.

m0a <- mareg(es~1, var=var, data=dat.sim.final, method = "ML") # EMPTY OMNIBUS MODEL
m3a <- mareg(es~dose+stress, var=var, data=dat.sim.final, method = "ML")  # FULL MODEL

anova(m0a, m3a) 
# NOTE: IN THIS SIMULATED DATA EXTRACTING A RANDOM SELECTION 
# OF K=8 (FROM A K=1000 POPULATION), EACH MODERATOR ACCOUNTS FOR
# ALL OF THE BETWEEN-STUDY VARIABILITY (REDUCING TAU^2 TO 0).
# THIS IS UNLIKELY TO HAPPEN WHEN EXAMINING REAL DATA WITH LARGER K.


####' DIAGNOSTICS
funnel(m0) # EXAMINE PUBLICATION BIAS VISUALLY
fsn(yi=m0$yi, vi=m0$vi)  # NO PUB BIAS

#### GRAPHICS  ####

# FIGURES FOR SIMULATED DATA
fig1a <- ggplot(dat.sim1, aes(x=g)) + 
  geom_histogram(binwidth=.1, colour="black", fill="blue") +
  geom_vline(aes(xintercept=median(g, na.rm=T)),   
             color="red", linetype="dashed", size=1)+
  ggtitle("Outcome One in Study Population")+
  xlab("g")+
  theme_bw()  

fig1b <- ggplot(dat.sim2, aes(x=g)) + 
  geom_histogram(binwidth=.1, colour="black", fill="yellow") +
  geom_vline(aes(xintercept=median(g, na.rm=T)),   
             color="red", linetype="dashed", size=1)+
  ggtitle("Outcome Two in Study Population")+
  xlab("g")+
  theme_bw()


grid.arrange(fig1a, fig1b, ncol=2, main="Effect Sizes Generated from Simulation")

# VISUALIZE THE HETEROGENEITY (fig2)

# FOR THE TUTORIAL GRAPHIC, CHANGE ID VARIABLES FOR BETTER DISPLAY 
dat.sim.final$id <- paste("Study", 1:length(dat.sim.final$id))

forest(m0)
text(-2.4 , 9.7, "Study",   pos=4,  cex=.7)
text(3.81,  9.7, "Observed g [95% CI]", pos=2, cex=.7)

# DOSE MOD (CONTINUOUS)

# USING FUNCTION IN MAd PACKAGE (HERE IT DOESNT MATCH UP WELL SO WE
# WILL USE THE CODE BELOW plotcon CODE)
plotcon(g=es, var=es, mod=dose, data=dat.sim.final, modname = "Treatment dose", 
        title = "Scatterplot examining dose moderator")

# LOOKS BETTER:
fig3a <- ggplot(data=dat.sim.final,aes(x=dose, y=es, size=1/var), na.rm=TRUE) +
  geom_point(alpha=.5) +   
  geom_abline(intercept=m1$b[1], slope=m1$b[2]) + 
  #ylim(0,1.25)+
  expand_limits(x = 5, y=0) +
  xlab("Treatment dose") + ylab("Effect Size") +  
  ggtitle("(a) Scatterplot examining dose moderator") +
  theme_bw() +
  theme(legend.position = "none")+
  scale_x_continuous(breaks= pretty_breaks())

# STRESS MOD (CATEGORICAL)
dat.sim.final$stress <- relevel(dat.sim.final$stress, ref="low")

# USING FUNCTION IN MAd PACKAGE
plotcat(g=es, var=es, mod=stress, data=dat.sim.final, modname = "Stress",
        title = "Boxplot examining stress moderator") 

# FOR ADDED FLEXIBILITY USE:
fig3b <- ggplot(dat.sim.final,  aes(factor(stress), es, weight = 1/var), na.rm=TRUE) + 
  geom_boxplot(aes(weight = 1/var), outlier.size=3,na.rm=TRUE) + 
  geom_jitter(aes(shape=factor(stress), color=factor(stress), size=1/var), alpha=.5) + 
  xlab("Stress") + 
  ylab("Effect Size")  +
  ggtitle("(b) Boxplot examining stress moderator") +
  theme_bw() +
  theme(legend.position = "none")


fig3c <- ggplot(dat.sim.final,  aes(factor(stress), dose), na.rm=TRUE) + 
  geom_boxplot(outlier.size=3,na.rm=TRUE) + 
  geom_jitter(aes(shape=factor(stress), color=factor(stress), size=1/var), alpha=.5) + 
  coord_flip()+
  xlab("Stress") + 
  ylab("Treatmen dose")  +
  ggtitle("(c) Confounding between moderators") +
  theme_bw() +
  theme(legend.position = "none")

fig3d <- ggplot(data=dat.sim.final,aes(x=dose, y=es, size=1/var), na.rm=TRUE) +
  geom_point(alpha=.5, aes(shape=stress, color=factor(stress))) +   
  geom_abline(intercept=m1$b[1], slope=m1$b[2]) + 
  #ylim(0,1.25)+
  expand_limits(x = 5, y=0) +
  xlab("Treatment dose") + ylab("Effect Size") +  
  ggtitle("(d) Effect sizes by dose and stress moderators") +
  theme_bw() +
  theme(legend.position = "none")+
  scale_x_continuous(breaks= pretty_breaks())

grid.arrange(fig3a, fig3b, fig3c, fig3d, ncol=2, main="Visualizing continuous and categorical moderator variables")


####  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