#metafor
library(metafor)

# set working directory 

metadat2 <- read.table("Dixit_Rinput_RawData_MetaAnalysis_HelpingEggSize.txt",header=T)

# funnel plot
dat<-metadat2
res_all <- rma.mv(yi=Z, V=Variance, data=dat, slab=StudyRef)
funnel(res_all, level=c(90, 95, 99), shade=c("white", "gray", "darkgray"), 
       yaxis="seinv",xlab="effect size (Zr)",ylab="precision (1/SE)",digits=c(1,0),xlim=c(-1,1),ylim=c(3,22))

#publication bias test
dat<-metadat2
res_pb <- rma(yi=Z, vi=Variance, data=dat, random=~1|StudyRef)
regtest(res_pb, model = "rma", predictor = "sei")

#summary of model
res_sum <- rma.mv(yi=Z, V=Variance, data=dat, slab=StudyRef, random=~1|StudyRef)
res_sum

#model with study species as random term
res_sum2 <- rma.mv(yi=Z, V=Variance, data=dat, slab=StudyRef, random=~1|SpeciesNo)
res_sum2 
       
       
       
#model with moderator of type of maternal post-natal response, with study reference as random term
res_sum1 <- rma.mv(yi=Z, V=Variance, data= metadat2, mods=~ Mothers_reduce_postnatal_investment, random=~1| StudyRef)
res_sum1

#model with moderator of type of maternal post-natal response, with study species as random term
res_sum2 <- rma.mv(yi=Z, V=Variance, data= metadat2, mods=~ Mothers_reduce_postnatal_investment, random=~1| SpeciesNo)
res_sum2

#model with moderator of type of maternal post-natal response, with no random term
res_sum3 <- rma.mv(yi=Z, V=Variance, data= metadat2, mods=~ Mothers_reduce_postnatal_investment)
res_sum3

#model of subset of data with mothers reducing post-natal care in helper presence
res_additive <- rma.mv(yi=Z, V=Variance, data=metadat2[metadat2$Mothers_reduce_postnatal_investment =="Yes",])
res_additive

#model of subset of data with mothers reducing post-natal care in helper presence, with study species as random term
res_additive_id <- rma.mv(yi=Z, V=Variance, data=metadat2[metadat2$Mothers_reduce_postnatal_investment =="Yes",],random=~1| SpeciesNo)
res_additive_id

#model of subset of data with mothers not reducing post-natal care in helper presence
res_compensated <- rma.mv(yi=Z, V=Variance, data=metadat2[metadat2$Mothers_reduce_postnatal_investment =="No",])
res_compensated

#model of subset of data with mothers not reducing post-natal care in helper presence
res_compensated_id <- rma.mv(yi=Z, V=Variance, data=metadat2[metadat2$Mothers_reduce_postnatal_investment =="No",],random=~1| SpeciesNo)
res_compensated_id



op <- par(oma=c(1.0,1.5,0.5,1.5), mar=c(4.0,4.0,2.2,1.2), mfrow=c(2,1))

funnel(res_additive, level=c(90, 95, 99), shade=c("white", "gray", "darkgray"), 
       yaxis="seinv",xlab="effect size (Zr)",ylab="precision (1/SE)",digits=c(1,0),xlim=c(-1,1),ylim=c(3,22))
       title(main="a) mothers decrease postnatal investment")

funnel(res_compensated, level=c(90, 95, 99), shade=c("white", "gray", "darkgray"), 
       yaxis="seinv",xlab="effect size (Zr)",ylab="precision (1/SE)",digits=c(1,0),xlim=c(-1,1),ylim=c(3,22))
       title(main="b) mothers do not decrease postnatal investment")








#model with moderator of type of total helping efforts with group size, with study reference as random term
metadat3<-metadat2[is.na(metadat2$Helper_benefits_increase_linearly)==FALSE,]
res_sum1 <- rma.mv(yi=Z, V=Variance, data= metadat3, mods=~ Helper_benefits_increase_linearly, random=~1| StudyRef)
res_sum1

#model with study species as random term
res_sum2 <- rma.mv(yi=Z, V=Variance, data= metadat3, mods=~ Helper_benefits_increase_linearly, random=~1| SpeciesNo)
res_sum2

#model with no random term
res_sum3 <- rma.mv(yi=Z, V=Variance, data= metadat3, mods=~ Helper_benefits_increase_linearly)
res_sum3

#model of subset of data with species in which total helping efforts increase linearly with group size
res_additive <- rma.mv(yi=Z, V=Variance, data= metadat3[metadat3$Helper_benefits_increase_linearly =="Yes",])
res_additive

#model of subset of data with species in which total helping efforts increase linearly with group size, with study species as random term
res_additive_id <- rma.mv(yi=Z, V=Variance, data= metadat3[metadat3$Helper_benefits_increase_linearly =="Yes",],random=~1| SpeciesNo)
res_additive_id

#model of subset of data with species in which total helping efforts does not increase linearly with group size
res_compensated <- rma.mv(yi=Z, V=Variance, data= metadat3[metadat3$Helper_benefits_increase_linearly =="No",])
res_compensated

#model of subset of data with species in which total helping efforts does not increase linearly with group size, with study species as random term
res_compensated_id <- rma.mv(yi=Z, V=Variance, data= metadat3[metadat3$Helper_benefits_increase_linearly =="No",],random=~1| SpeciesNo)
res_compensated_id