#Load Package
install.packages("metafor")
library(metafor)

install.packages("forestplot")
library(forestplot)

#DataSet
mydata <- read.delim("~/Desktop/anxietydata")

#Describe Data
anxr<-(mydata$AnxietyR)
N<-(mydata$N)
ID<-(mydata$Number)
age<-(mydata$Mean.Age)
Female<-(mydata$PercentageFemale)
SampleType<-(mydata$Percentage.of.Clinical)

#META-ANALYSIS OF IU and ANXIETY

#Calculate Fisher's Z and V for IU and Anxiety
#yi valueas are r to z transformed correlations
#vi values the corresponding sampling variances
ESANX<- escalc(measure = "ZCOR", ri=anxr, ni=N, data = mydata)


#Random Effect Model for IU and Anxiety
REanx<- rma(yi=ESANX$yi, vi=ESANX$vi, data=mydata, digit=3)
predict(REanx, transf = transf.ztor)


forest.rma(REanx, slab = paste(ESANX$Author, ", ", ESANX$Date, sep = ""),showweights = TRUE, transf = transf.ztor)

funnel(REanx)

#Fisher r-to-Z transformation-conversion back to r
transf.ztor(0.692)

#leaveout 
#A single outlying trial could be the source of substantialheterogeneity.
#To identify suspicious cases, a leave-one-out method can be usedwhereby 
#we rerun the meta-analysis, iteratively removing studies.
casesanx<- leave1out(REanx)
which(casesanx$I2 ==min(casesanx$I2))
sum(casesanx$I2<30)



# age as a moderator
MEanxietyage<-rma(yi=ESANX$yi, vi=ESANX$vi, mods= ~(ESANX$Mean.Age), data=mydata) 
predict(MEanxietyage, transf = transf.ztor)

# gender as a moderator
MEanxietygender<-rma(yi=ESANX$yi, vi=ESANX$vi, mods= ~(ESANX$PercentageFemale), data=mydata)

# study location as a moderator
MEanxietylocation<-rma(yi=ESANX$yi, vi=ESANX$vi, mods= ~factor(ESANX$Study.Location), data=mydata)

# IU measure as a moderator
MEanxietyIUmeasure<-rma(yi=ESANX$yi, vi=ESANX$vi, mods= ~factor(ESANX$IU.Measure), data=mydata)  

# anxiety measure as a moderator
MEanxietyAmeasure<-rma(yi=ESANX$yi, vi=ESANX$vi, mods= ~factor(ESANX$Anxiety.Measure), data=mydata)
forest(MEanxietyAmeasure, slab = ESANX$Author,showweights = TRUE)


# sample type as a moderator (percentage of clinical)
MEanxietysampletype<-rma(yi=ESANX$yi, vi=ESANX$vi, mods= ~(ESANX$Percentage.of.Clinical), data=mydata)  

#Rosenthal's Fail-Safe N
fsn(yi=ESANX$yi, vi=ESANX$vi, data=mydata, type = "Rosenthal", alpha = .05)

# Egger Regression test
regtest(REanx, model = "rma", predictor = "sei")
regtest(REanx, model = "rma", predictor = "vi")
regtest(REanx, model = "rma", predictor = "ni")
regtest(REanx, model = "rma", predictor = "ninv")
regtest(REanx, model = "rma", predictor = "sqrtni")
regtest(REanx, model = "rma", predictor = "sqrtninv")

#Begg and Mazumdar rank correlation test
ranktest(REanx)