#### Introduction ####

#Meta Analysis Example
#Difference between Males and Females 
#in Maladaptive Perfectionism

dat<-read.csv(file.choose())
names(dat)
head(dat)

library(metafor)

#### Fixed Effects Model ####

#Fixed Effects Model (method="FE")
#Raw mean difference (measure="MD")

ma1<-rma(data=dat,
         n1i=N_male,n2i=N_female,
         sd1i=SD_male,sd2i=SD_female,
         m1i=M_male,m2i=M_female,
         method="FE",measure="MD")
summary(ma1)

#Fixed Effects Model (method="FE") 
#Standardized Mean Difference (measure="SMD")

ma2<-rma(data=dat,
         n1i=N_male,n2i=N_female,
         sd1i=SD_male,sd2i=SD_female,
         m1i=M_male,m2i=M_female,
         method="FE",measure="SMD")
summary(ma2)

#### Random Effects Model ####

#Random Effects Model
#Restricted Maximum Likelihood Estimation (method="REML")
#Raw Mean Difference (measure="MD")

ma3<-rma(data=dat,
         n1i=N_male,n2i=N_female,
         sd1i=SD_male,sd2i=SD_female,
         m1i=M_male,m2i=M_female,
         method="REML",measure="MD")
summary(ma3)

#Random Effects Model
#Restricted Maximum Likelihood Estimation (method="REML")
#Raw Mean Difference (measure="SMD")

ma4<-rma(data=dat,
         n1i=N_male,n2i=N_female,
         sd1i=SD_male,sd2i=SD_female,
         m1i=M_male,m2i=M_female,
         method="REML",measure="SMD")
summary(ma4)


#### Some "hand" Calculations ####

#Calculating Q/tau-squared "by hand"
#DerSimonian Laird Method
dat_es<-escalc(data=dat,n1i=N_male,n2i=N_female,sd1i=SD_male,sd2i=SD_female,
       m1i=M_male,m2i=M_female, measure="SMD")

#Weights
w<-(1/dat_es$vi)

#Weighted Mean Effect Size
mes<-sum(w*dat_es$yi)/sum(w)

#Q statistic
(Q<-sum(w*(dat_es$yi-mes)^2))
1-pchisq(Q,length(w)-1)

#Calculate tau/tau-squared
cc<-sum(w)-(sum(w^2)/sum(w))
tau_sq<-(Q-(length(w)-1))/cc
tau_sq
sqrt(tau_sq)

#% of between study variance that is true heterogeneity
(Isq<-((Q-(length(w)-1))/Q)*100)

#Check calculations, note method="DL"
ma5<-rma(data=dat,n1i=N_male,n2i=N_female,sd1i=SD_male,sd2i=SD_female,
         m1i=M_male,m2i=M_female,method="DL",measure="SMD")
summary(ma5)


#### Diagnostics on the Random Effect Model ####

#Forest Plot for the Standardized Mean Difference
forest(ma4,level=95)

#Identify Influential Cases
ic<-influence(ma4)
plot(ic)
(cd_cut<-4/(length(dat$Study)))
any(ic$inf$cook.d > cd_cut)

#Influential Cases
plot(leave1out(ma4)$estimate)

#Publication Bias - Funnel Plot
funnel(ma4)
regtest(ma4) #Test for Funnel Plot Asymmetry

#Explore the effect of Filling in Studies 
#Missing Due to Publication Bias
tf<-trimfill(ma4)
funnel(tf)

# Original and Trim & Fill Confidence Intervals on
# Cohen's D
c(ma4$ci.lb,ma4$ci.ub)
c(tf$ci.lb,tf$ci.ub)

#Four helpful plots with one command
plot(ma4)



#### Moderator Analyses ####

library(ggplot2)

#Quality as a Moderator of Effect Size?
ggplot(data=dat_es, aes(x=Quality,y=yi)) + geom_point(size=5)

#Perfectionism Measure as a Moderator of Effect Size
boxplot(dat_es$yi ~ dat_es$Perf_Measure)

#Random Effects Model
#Restricted Maximum Likelihood Estimation (method="REML")
#Raw Mean Difference (measure="SMD")
#With moderator centred "Perfectionism Measure"

ma6<-rma(data=dat,
         n1i=N_male,n2i=N_female,
         sd1i=SD_male,sd2i=SD_female,
         m1i=M_male,m2i=M_female,
         method="REML",measure="SMD",
         mods=~Perf_Measure)
summary(ma6)

#Forest Plot with Fitted Values
forest(ma6)
forest(ma6,order="obs")

plot(ma6)

# Likelihood Ratio Test for Moderator "Perf_Measure"
# Must use ML, not REML

ma7a<-rma(data=dat,
          n1i=N_male,n2i=N_female,
          sd1i=SD_male,sd2i=SD_female,
          m1i=M_male,m2i=M_female,
          method="ML",measure="SMD",
          mods=NULL)
ma7b<-rma(data=dat,
          n1i=N_male,n2i=N_female,
          sd1i=SD_male,sd2i=SD_female,
          m1i=M_male,m2i=M_female,
          method="ML",measure="SMD",
          mods=~Perf_Measure)
anova(ma7a,ma7b)


#Restricted Maximum Likelihood Estimation (method="REML")
#Raw Mean Difference (measure="SMD")
#With moderator centred "Study Quality"
dat$Quality<-dat$Quality-mean(dat$Quality)

ma8<-rma(data=dat,
         n1i=N_male,n2i=N_female,
         sd1i=SD_male,sd2i=SD_female,
         m1i=M_male,m2i=M_female,
         method="REML",measure="SMD",
         mods=~Quality)
summary(ma8)

# Likelihood Ratio Test for Moderator "Quality"
ma8a<-rma(data=dat,n1i=N_male,n2i=N_female,sd1i=SD_male,sd2i=SD_female,
          m1i=M_male,m2i=M_female,method="ML",measure="SMD",
          mods=NULL)
ma8b<-rma(data=dat,
          n1i=N_male,n2i=N_female,
          sd1i=SD_male,sd2i=SD_female,
          m1i=M_male,m2i=M_female,
          method="ML",measure="SMD",
          mods=~Quality)
anova(ma8a,ma8b)


#Random Effects Model
#Restricted Maximum Likelihood Estimation (method="REML")
#Raw Mean Difference (measure="SMD")
#With moderators "Study Quality" and "Perf_Measure"

ma9<-rma(data=dat,
         n1i=N_male,n2i=N_female,
         sd1i=SD_male,sd2i=SD_female,
         m1i=M_male,m2i=M_female,
         method="REML",measure="SMD",
         mods=~Quality+Perf_Measure)
summary(ma9)

# Likelihood Ratio Test for Moderator "Quality" over and
# above Moderator "Perf_Measure"
ma10<-rma(data=dat,
          n1i=N_male,n2i=N_female,
          sd1i=SD_male,sd2i=SD_female,
          m1i=M_male,m2i=M_female,
          method="ML",measure="SMD",
          mods=~Quality+Perf_Measure)
anova(ma7b,ma10)

# Likelihood Ratio Test for Moderator "Perf_Measure" over and
# above Moderator "Quality"
ma10<-rma(data=dat,
          n1i=N_male,n2i=N_female,
          sd1i=SD_male,sd2i=SD_female,
          m1i=M_male,m2i=M_female,
          method="ML",measure="SMD",
          mods=~Quality+Perf_Measure)
anova(ma8b,ma10)

#### The End! ####