library(metafor)
library(xlsx)#Esse pacote precisa ter o java 64-bits instalado no computador
library(sciplot)
library(doBy)
library(lme4)
library(ggplot2)
library(Hmisc)
library(Rmisc)
library(vegan)
library(scales)
library(tidyverse)

#===========================================================================================================================================================
## NATURAL ENEMIES
rm(list = ls())

data=read.table("enemy.txt", sep = "\t", h=T, dec=",")
summary(data)
View(data)
nrow(data)
(id_es=paste("ID", 1:nrow(data), sep = "_"))

#Calculating the "overall effect" of ant presence on NE.
overall=rma.mv(yi, vi, random= list(~1|authors/id_es), struct="CS", method="REML", digits=5, data=data)
summary(overall)

### I? 
800.62-239
(561.62/800.62)*100

tiff(file="Funnel ne.tif", width=3200, height=3200, res=600, compression = "lzw")
funnel(overall)
dev.off()

#Egger test
eggs= rma.mv(yi,vi, mods=~vi,random= list(~1|authors/id_es), struct="CS", method="REML", digits=5, data=data)
summary(eggs)

#### Pest type as moderators
model_pest=rma.mv(yi, vi, mods= ~pest_type, random=list(~1|authors/id_es), struct="CS", method="REML", digits=5, data=data)
summary(model_pest)

#########Experiment time
md1<-rma.mv(yi, vi, mods= ~experiment_time, random=list(~1|authors/id_es), struct="CS", method="REML", data=data)
summary(md1)

#### Crop size
data$crop_size <- as.numeric(data$crop_size)
str(data)
cs=rma.mv(yi, vi, mods= ~crop_size, random=list(~1|authors/id_es), struct="CS", method="REML", digits=5, data=data)
summary(cs)

#### Enemies movement
md2=rma.mv(yi, vi, mods= ~mov, random=list(~1|authors/id_es), struct="CS", method="REML", digits=5, data=data)
summary(md2)

#### Enemies specialization
md3=rma.mv(yi, vi, mods= ~enemy_group, random=list(~1|authors/id_es), struct="CS", method="REML", digits=5, data=data)
summary(md3)

#Show estimate
md3=rma.mv(yi, vi, mods= ~enemy_group-1, random=list(~1|authors/id_es), struct="CS", method="REML", digits=5, data=data)
summary(md3)

###########Ant Size##############
data$ant_size <- as.numeric(data$ant_size)
testout2=subset(data, data$ant_size!="NA")
(id_es=paste("ID", 1:nrow(testout2), sep = "_"))
as=rma.mv(yi, vi, mods= ~ant_size, random=list(~1|authors/id_es), struct="CS", method="REML", digits=5, data=testout2)
summary(as)

#########################  PLANT DAMAGE #########################################################
rm(list = ls())

data=read.table("damage.txt", sep = "\t", h=T, dec=",")
summary(data)
View(data)
(id_es=paste("ID", 1:nrow(data), sep = "_"))
#Calculating the "overall effect" of ant presence on pest
overall=rma.mv(yi, vi, random= list(~1|authors), struct="UN", method="REML", digits=5, data=data)
summary(overall)

### I? 
1865.06-176
(1689.06/1865.06)*100

tiff(file="Funnel pd.tif", width=3200, height=3200, res=600, compression = "lzw")
funnel(overall)
dev.off()

#Egger test.
eggs= rma.mv(yi,vi, mods=~vi,random= list(~1|authors/id_es), struct="CS", method="REML", digits=5, data=data)
summary(eggs)

data$crop_size <- as.numeric(data$crop_size)
data$experiment_time <- as.numeric(data$experiment_time)

#### Crop size as moderators
mp=rma.mv(yi, vi, mods= ~crop_size, random=list(~1|authors/id_es), struct="CS", method="REML", digits=5, data=data)
summary(mp)

#########Experiment time
md1<-rma.mv(yi, vi, mods= ~experiment_time, random=list(~1|authors/id_es), struct="CS", method="REML", data=data)
summary(md1)

#### Pest type (honeydew Vs non-honeydew)
mpest1=rma.mv(yi, vi, mods= ~pest_type, random=list(~1|authors/id_es), struct="CS", method="REML", digits=5, data=data)
summary(mpest1)

#Show estimate
mpest1=rma.mv(yi, vi, mods= ~pest_type-1, random=list(~1|authors/id_es), struct="CS", method="REML", digits=5, data=testout)
summary(mpest1)

#### Crop system 
md2=rma.mv(yi, vi, mods= ~crop_system, random=list(~1|authors/id_es), struct="CS", method="REML", digits=5, data=data)
summary(md2)

#Show estimate
md2=rma.mv(yi, vi, mods= ~crop_system-1, random=list(~1|authors/id_es), struct="CS", method="REML", digits=5, data=data)
summary(md2)

#### Pest group 
mpest=rma.mv(yi, vi, mods= ~pest, random=list(~1|authors/id_es), struct="CS", method="REML", digits=5, data=data)
summary(mpest)

#Show estimate
mpest=rma.mv(yi, vi, mods= ~pest-1, random=list(~1|authors/id_es), struct="CS", method="REML", digits=5, data=data)
summary(mpest)

###########Ant Size##############
data$ant_size <- as.numeric(data$ant_size)
testout1=subset(data, data$ant_size!="NA")

(id_es=paste("ID", 1:nrow(testout1), sep = "_"))
as_pd=rma.mv(yi, vi, mods= ~ant_size, random=list(~1|authors/id_es), struct="CS", method="REML", digits=5, data=testout1)
summary(as_pd)


######################################### CROP YIELD ######################################################
rm(list = ls())

data=read.table("yield.txt", sep = "\t", h=T, dec=",")
summary(data)
View(data)
(id_es=paste("ID", 1:nrow(data), sep = "_"))

#Calculating the "overall effect" of ant presence on pest
overall=rma.mv(yi, vi, random= list(~1|authors/id_es), struct="CS", method="REML", digits=5, data=data)
summary(overall)

### I2 
226.73-72
(154.73/226.73)*100

tiff(file="Funnel pd.tif", width=3200, height=3200, res=600, compression = "lzw")
funnel(overall)
dev.off()

#Egger test.
eggs= rma.mv(yi,vi, mods=~vi,random= list(~1|authors/id_es), struct="CS", method="REML", digits=5, data=data)
summary(eggs)

#########Experiment time
data$experiment_time <- as.numeric(data$experiment_time)
met<-rma.mv(yi, vi, mods= ~experiment_time, random=list(~1|authors/id_es), struct="CS", method="REML", data=data)
summary(met)

###########Ant Size##############
data$ant_size <- as.numeric(data$ant_size)
as_cy=rma.mv(yi, vi, mods= ~ant_size, random=list(~1|authors/id_es), struct="CS", method="REML", digits=5, data=data)
summary(as_cy)

#### Crop system 
mcs=rma.mv(yi, vi, mods= ~crop_system, random=list(~1|authors/id_es), struct="CS", method="REML", digits=5, data=data)
summary(mcs)

#Show estimate
mcs1=rma.mv(yi, vi, mods= ~crop_system-1, random=list(~1|authors/id_es), struct="CS", method="REML", digits=5, data=data)
summary(mcs1)