########################################################################################
################################# Microcosm Experiment #################################
########################################################################################

#####################################################################
########################### Survival Data ###########################

#####Loading dataset#####
setwd("C:\\Users\\HP\\Ian\\Active\\Mantid predation study") #Change this line to match the directory where the datafile is stored
d=read.table("Data_Survival.txt",header=T)
str(d)
d$mantisPresence=factor(d$mantisPresence)
##Coding the censorship variable correctly, 0: alive, 1:dead
a=rep(0,length(d$Censored))
a[d$Censored=="N"]=1
d$Censored=a

#####Survival analysis#####
##Load packages
require(survival)
require(ggfortify)
require(survminer)

##Comparing experimentals and controls
#Looking for outliers
boxplot(Survival~mantisPresence, data=d)
d=d[-which(d$Survival==334),] #remove an extreme outlier (334 days; next largest is 92)
#Running the analysis
smod1=survreg(Surv(Survival,Censored)~mantisPresence,dist="exponential",data=d) #model with constant hazard
smod2=survreg(Surv(Survival,Censored)~mantisPresence,data=d) #model with non-constant hazard, default uses Weibull distribution
AIC(smod1,smod2) #non-constant hazard is better (lower AIC value)
smod1a=survreg(Surv(Survival,Censored)~mantisPresence,dist="loglogistic",data=d) #trying other hazard distributions
AIC(smod1a,smod2) #loglogistic is best (also tried gaussian, logistic, loglogistic)
summary(smod1a) 
#Result: Mantid presence is significant (p = 2e-16)
#Plotting survival curves
sObj=Surv(d$Survival,d$Censored)
smod2plot=survfit(sObj~mantisPresence,data=d)
autoplot(smod2plot)
##Final result: butterflies in experimental cages survived much shorter than in control (6.7 vs 28.7 days; p < 0.001)

##Comparing Forms in control trials (without mantises)
#Creating dataframe
d2=d[d$mantisPresence=="0",]
d2=droplevels(d2)
#Looking for outliers
boxplot(Survival~Form, data=d2)
d2=d2[-which(d2$Survival==334),] #remove an extreme outlier (334 days; next largest is 92)
#Running the analysis
smod3=survreg(Surv(Survival,Censored)~Form+cluster(Cage),dist="exponential",data=d2) #model with constant hazard
smod4=survreg(Surv(Survival,Censored)~Form+cluster(Cage),data=d2) #model with non-constant hazard, default uses Weibull distribution
AIC(smod3,smod4) #non-constant hazard is better (lower AIC value)
smod3a=survreg(Surv(Survival,Censored)~Form+cluster(Cage),dist="loglogistic",data=d2) #trying other hazard distributions
AIC(smod3a,smod4) #loglogistic is best (also tried gaussian, logistic, loglogistic)
summary(smod3a) 
#Result: Form is not significant in control trials (p = 0.48)
#Plotting survival curves
sObj2=Surv(d2$Survival,d2$Censored)
smod3aplot_form=survfit(sObj2~Form,data=d2)
autoplot(smod3aplot_form)+scale_colour_manual(values = c('chocolate4','darkolivegreen'))+scale_fill_manual(values = c('tan1','darkolivegreen'))+xlim(c(0,100))
##Final result: there is no difference in longevity between the Spotty (27.2 days) and Wt (29.4 days) in control trials (p = 0.48)

##Comparing Forms in experimental trials (with mantises)
#Creating dataframe
d2=d[d$mantisPresence=="1",]
d2=droplevels(d2)
#Looking for outliers
boxplot(Survival~Form, data=d2) #no extreme outliers
#Running the analysis
smod3=survreg(Surv(Survival,Censored)~Form*mantisGender+cluster(Cage),dist="exponential",data=d2) #model with constant hazard
smod4=survreg(Surv(Survival,Censored)~Form*mantisGender+cluster(Cage),data=d2) #model with non-constant hazard, default uses Weibull distribution
AIC(smod3,smod4) #non-constant hazard is better (lower AIC value)
smod3a=survreg(Surv(Survival,Censored)~Form*mantisGender+cluster(Cage),dist="loglogistic",data=d2)
AIC(smod3a,smod4) #loglogistic is best (also tried gaussian, logistic, loglogistic)
summary(smod3a) #Result: interaction is not significant (p = 0.219), remove from model
smod5=survreg(Surv(Survival,Censored)~Form+mantisGender+cluster(Cage),dist="loglogistic",data=d2)
summary(smod5)
#Results: both Form (p = 0.032) and Mantis Gender (p < 0.01) significant
#Plotting survival curves
sObj2=Surv(d2$Survival,d2$Censored)
smod5plot_form=survfit(sObj2~Form,data=d2) #plot Form only
autoplot(smod5plot_form)+scale_colour_manual(values = c('chocolate4','darkolivegreen'))+scale_fill_manual(values = c('tan1','darkolivegreen'))+xlim(c(0,100))
smod5plot_gender=survfit(sObj2~mantisGender,d2) #plot Mantis Gender only
autoplot(smod5plot_gender)
smod5plot_both=survfit(sObj2~interaction(Form,mantisGender),data=d2) #plot the effects of both
autoplot(smod5plot_both)
##Final result: both Form (Spotty 3.4 days, Wt 9.4 days; p = 0.032) and Mantis Gender (Male 5.7 days, Female 3.3 days; p < 0.01) are significant in experimental trials, but do not interact

######################## Survival Data Ends #########################
#####################################################################

#####################################################################
########################## Fecundity Data ###########################

#####Loading dataset#####
setwd("C:\\Users\\HP\\Ian\\Active\\Mantid predation study") #Change this line to match the directory where the datafile is stored
d3=read.table("Data_Fecundity.txt",header=T)
str(d3)

#####GLMM#####
##Load packages
require(glmmTMB)
require(DHARMa)

##Comparing experimentals and controls
#Looking for outliers
boxplot(EggsButtDay~Mantis, data=d3) #no extreme outliers
#Running the GLMM
mod3=glmmTMB(EggsButtDay~Form*Mantis+(1|Cage),family=gaussian,data=d3) #warning caused by 1|Cage, but model still looks okay and is essentially the same as without the term
mod3sim=simulateResiduals(mod3) #Model checking
plot(mod3sim) #looks OK
summary(mod3) #The interaction is non-significant, to be removed
mod4=glmmTMB(EggsButtDay~Form+Mantis+(1|Cage),family=gaussian,data=d3)
mod4sim=simulateResiduals(mod4)
plot(mod4sim) #looks OK
summary(mod4) #Mantis Presence is significant, Form is marginal (p = 0.0516) and removed
mod5=glmmTMB(EggsButtDay~Mantis+(1|Cage),family=gaussian,data=d3)
mod5sim=simulateResiduals(mod5)
plot(mod5sim) #looks OK
summary(mod5)
##Final result: Mantis presence is highly significant (p = 3.37e-5; 9.3 treatments vs 18.9 controls)

##Comparing Forms in control trials (without mantises)
#Creating dataframe
d3b=d3[d3$Mantis=="N",]
d3b=droplevels(d3b)
#Looking for outliers
boxplot(EggsButtDay~Form, data=d3b) #no extreme outliers
#Running the GLMM
mod7=glmmTMB(EggsButtDay~Form+(1|Cage),family=gaussian,data=d3b) #warning caused by 1|Cage, but model still looks okay and is essentially the same as without the term
mod7sim=simulateResiduals(mod7)
plot(mod7sim) #looks OK
summary(mod7) #nothing significant
#Note: Curves plotted manually from data
##Final result: Form is not sig (Wt 20.3, Spotty 17.5; p = 0.381)

##Comparing Forms in experimental trials (with mantises)
#Creating dataframe
d3a=d3[d3$Mantis=="Y",]
d3a=droplevels(d3a)
#Looking for outliers
boxplot(EggsButtDay~Form, data=d3a) #no extreme outliers
#Running the GLMM
mod5=glmmTMB(EggsButtDay~Form*MantisGender+(1|Cage),family=gaussian,data=d3a) #warning with 1|Cage, but model still looks okay and is essentially the same as without the term
mod5sim=simulateResiduals(mod5) #Model checking
plot(mod5sim) #looks OK
summary(mod5) #nothing significant, remove interaction term
mod6=glmmTMB(EggsButtDay~Form+MantisGender+(1|Cage),family=gaussian,data=d3a)
mod6sim=simulateResiduals(mod6) #Model checking
plot(mod6sim) #looks OK
summary(mod6) #all significant, no further simplification possible
#Note: Curves plotted manually from data
#Final result: Form is sig (Wt 13.4, Spotty 6.8; p = 0.0216); Mantis Gender is sig (Male 13.3, Female 7.0; p = 0.0337)

####################### Fecundity Data Ends #########################
#####################################################################

#####################################################################
######################### Wing Damage Data ##########################

#####Loading dataset#####
setwd("C:\\Users\\HP\\Ian\\Active\\Mantid predation study") #Change this line to match the directory where the datafile is stored
d_fracs=read.table("Data_WingDamage.txt",header=T) #using WingScoringFractions
colnames(d_fracs) #8:21 are Forewing,22:32 are Hindwing
dtemp=d_fracs[,8:32] #Finding the min and max values to have an idea of how to transform data
dtemp[dtemp==0]=NA
min(dtemp,na.rm=T) #0.17
max(dtemp,na.rm=T) #5
d_fracs[,8:32]=floor(d_fracs[,8:32]*100) #Transforming data to do neg.binom

#####Multivariate GLM#####
##Load packages
require(mvabund)

##Comparing experimentals and controls
#Creating dataframe
d=d_fracs
##Forewings
d_fore=d[d$BflyWing=="Fore",c(1:21,33:35)]
#Creating matrix for variables
exV=d_fore[,c(1:7,22)]
deV=mvabund(d_fore[,8:21])
#Doing the GLM
m1=manyglm(deV~MantidPresence*ModalSurvival+Cage,data=d_fore,family="negative.binomial")
plot.manyglm(m1) #looks good
anova(m1,p.uni="none",nBoot=99) #took 6 min 54 sec
#Results: MantidPresence p = 0.01, Cage p=0.01, ModalSurvival = 0.19, interaction 0.78
#Decision: remove interaction
m2=manyglm(deV~MantidPresence+ModalSurvival+Cage,data=d_fore,family="negative.binomial")
plot.manyglm(m2) #looks good
anova(m2,p.uni="none",nBoot=99) #took 4 min 25 sec
#Results: MantidPresence p = 0.01, Cage p=0.01, ModalSurvival = 0.13
#Decision: Remove ModalSurvival
m3=manyglm(deV~MantidPresence+Cage,data=d_fore,family="negative.binomial")
plot.manyglm(m3) #looks good
anova(m3,p.uni="none",nBoot=99) #took 2 min 59 sec
#Results: MantidPresence p = 0.01, Cage p=0.01
#Final result: Significant difference between Controls and Experimentals on forewings, Cage is a "random effect"
##Hindwings
d_hind=d[d$BflyWing=="Hind",c(1:7,22:32,33:35)]
#Creating matrix for variables
exV=d_hind[,c(1:7,19)]
deV=mvabund(d_hind[,8:18])
#Doing the GLM
m1=manyglm(deV~MantidPresence*ModalSurvival+Cage,data=d_hind,family="negative,binomial")
plot.manyglm(m1) #looks OK
anova.manyglm(m1,p.uni="none",nBoot=99) #took 7 min 34 sec
#Results: MantidPresence p = 0.01, ModalSurvival p=0.01, Cage p=0.01, interaction p=0.93
#Decision: Remove interaction
m2=manyglm(deV~MantidPresence+ModalSurvival+Cage,data=d_hind,family="negative,binomial")
plot.manyglm(m2) #looks OK
anova.manyglm(m2,p.uni="none",nBoot=99) #took 4 min 44 sec
#Results: MantidPresence p = 0.01, ModalSurvival p=0.01, Cage p=0.01
#Decision: No further simplification possible
#Final result: Significant difference between Controls and Experimentals on hindwings, Cage and ModalSurvival (both sig) are "random effects"
##Summary: Significant difference between Controls and Experimentals on both wings

##Comparing Forms in control trials (without mantises)
#Creating dataframe
d_fracsC=d_fracs[d_fracs$MantidPresence=="Absent",] 
d=d_fracsC
##Forewings
d_fore=d[d$BflyWing=="Fore",c(1:21,33:35)]
#Creating matrix for variables
exV=d_fore[,c(1:7,22)]
deV=mvabund(d_fore[,8:21])
#Doing the GLM
m1=manyglm(deV~Cage+BflyForm*ModalSurvival,data=d_fore,family="negative.binomial")
plot.manyglm(m1) #looks good
anova.manyglm(m1,p.uni="none",nBoot=99) #took 3 min 53 sec
#Results: Cage p=0.01, BflyForm p = 0.82, Survival p=0.70, Form:Survival p=0.83
#Decision: remove interaction
m2=manyglm(deV~Cage+BflyForm+ModalSurvival,data=d_fore,family="negative.binomial")
plot.manyglm(m2) #looks good
anova.manyglm(m2,p.uni="none",nBoot=99) #took 2 min 46 sec
#Results: Cage p=0.01, BflyForm p = 0.80, Survival p=0.74
#Decision: remove Form
m3=manyglm(deV~Cage+ModalSurvival,data=d_fore,family="negative.binomial")
plot.manyglm(m3) #looks good
anova.manyglm(m3,p.uni="none",nBoot=99) #took 1 min 44 sec
#Results: Cage p=0.01, Survival p = 0.83
#Decision: Remove Survival
m4=manyglm(deV~Cage,data=d_fore,family="negative.binomial")
plot.manyglm(m4) #looks good
anova.manyglm(m4,p.uni="none",nBoot=99) #took 45 sec
#Results: Cage p=0.01
#Decision: No further simplification possible
#Final result: No difference between Forms on forewings in controls, Cage though sig is a "random effect"
##Hindwings
d_hind=d[d$BflyWing=="Hind",c(1:7,22:32,33:35)]
#Creating matrix for variables
exV=d_hind[,c(1:7,19)]
deV=mvabund(d_hind[,8:18])
#Doing the GLM
m1=manyglm(deV~BflyForm*ModalSurvival+Cage,data=d_hind,family="negative,binomial")
plot.manyglm(m1) #looks OK
anova.manyglm(m1,p.uni="none",nBoot=99) #took 3 min 52 sec
#Results: BflyForm p = 0.85, Survival p=0.01, Cage p=0.01, Form:Survival p=0.50
#Remove Form:Survival
m2=manyglm(deV~BflyForm+ModalSurvival+Cage,data=d_hind,family="negative,binomial")
plot.manyglm(m2)
anova.manyglm(m2,p.uni="none",nBoot=99) #took 2 min 44 sec
#Results: BflyForm p = 0.83, Survival p=0.01, Cage p=0.03
#Decision: Remove Form
m3=manyglm(deV~ModalSurvival+Cage,data=d_hind,family="negative,binomial")
plot.manyglm(m3)
anova.manyglm(m3,p.uni="none", nBoot=99) #took 1 min 50 sec
#Results: Survival p=0.01, Cage p=0.04
#Decision: No further simplification possible
#Final result: No difference between Forms on hindwings in controls, Survival and Cage though sig are "random effects"
##Summary: no difference between Spotty and Wildtype in Controls on both fore and hindwings

##Comparing Forms in experimental trials (with mantises)
#Creating dataframe
d_fracsE=d_fracs[d_fracs$MantidPresence=="Present",] 
d=d_fracsE
##Forewings
d_fore=d[d$BflyWing=="Fore",c(1:21,33:35)]
#Creating matrix for variables
exV=d_fore[,c(1:7,22)]
deV=mvabund(floor(d_fore[,8:21]/d_fore[,24]*10))
#Doing the GLM
m3=manyglm(deV~MantidSex*BflyForm,data=d_fore,family="negative.binomial")
plot.manyglm(m3) #looks good
anova(m3,p.uni="none",nBoot=99) #took 17 sec
#Results: MantidSex p=0.4, Form p=0.13, Sex:Form p=0.01
#Decision: Remove interaction
m4=manyglm(deV~MantidSex+BflyForm,data=d_fore,family="negative.binomial")
plot.manyglm(m4) #looks good
anova(m4,p.uni="none",nBoot=99) #took 2 min 46 sec
#Results: MantidSex p=0.333, Form p=0.083
#Decision: Remove MantidSex
m5=manyglm(deV~BflyForm,data=d_fore,family="negative.binomial")
plot.manyglm(m5) #looks good
anova(m5,p.uni="none",nBoot=99) #took 52 sec
#Results: Form p=0.083
#Decision: Compare to a null model
m6=manyglm(deV~1,data=d_fore,family="negative.binomial")
plot.manyglm(m6) #looks good
anova(m5,m6,nBoot=99)
#Results: Form p=0.04
#Decision: Model with Form is better than null model
#Final result: Difference between Wt and Spotty on forewings in Experimental cages
##Hindwings
d_hind=d[d$BflyWing=="Hind",c(1:7,22:32,33:35)]
#Creating matrix for variables
exV=d_hind[,c(1:7,19)]
deV=mvabund(d_hind[,8:18])
#Doing the GLM
m1=manyglm(deV~BflyForm*ModalSurvival+BflyForm*MantidSex,data=d_hind,family="negative,binomial")
plot.manyglm(m1) #looks OK
anova.manyglm(m1,p.uni="none,nBoot=99) #took 1 min 40 sec
#Results: BflyForm p = 0.31, MantidSex p = 0.02, Survival p=0.02, Form:Sex p=0.80, Form:Survival p=0.42
#Decision: remove Form:Sex
m2=manyglm(deV~BflyForm*ModalSurvival+MantidSex,data=d_hind,family="negative,binomial")
plot.manyglm(m2) #looks OK
anova.manyglm(m2,p.uni="none",nBoot=99) #took 1 min 16 sec
#Results: BflyForm p = 0.25, Survival p=0.01, MantidSex p = 0.02, Form:Survival p=0.32
#Decision: remove Form:Survival
m3=manyglm(deV~BflyForm+ModalSurvival+MantidSex,data=d_hind,family="negative,binomial")
plot.manyglm(m3) #looks OK
anova.manyglm(m3,p.uni="none") #took 52 sec
#Results: BflyForm p = 0.19, Survival p=0.01, MantidSex p = 0.06
#Decision: remove Form
m4=manyglm(deV~ModalSurvival+MantidSex,data=d_hind,family="negative,binomial")
plot.manyglm(m4) #looks OK
anova.manyglm(m4,p.uni="none",nBoot=99) #took 33 sec
#Results: Survival p=0.03, MantidSex p = 0.03
#Decision: Survival and Mantid Sex significant, no further simplification possible
#Final result: No difference between Wt and Spotty on hindwings in Experimental cages
##Summary: Signficant difference between Wt and Spotty on forewings but not hindwings

####################### Wing Damage Data Ends #######################
#####################################################################

########################################################################################
############################## Microcosm Experiment Ends ###############################
########################################################################################


########################################################################################
################################### Arena Experiment ###################################
########################################################################################

#####################################################################
######################### First Strike Data #########################

#####Loading dataset#####
setwd("C:\\Users\\HP\\Ian\\Active\\Mantid predation study") #Change this line to match the directory where the datafile is stored
d=read.table("Data_FirstStrike.txt",header=T)
d$AttLoc=as.factor(d$AttLoc)

#####Multinomial model#####
##Load packages
require(nnet)
require(car)

##Multinomial analysis
#Run the analysis
mod1=multinom(AttLoc~Form,data=d)
Anova(mod1) #overall P = 0.06262, marginally sig., look for what is driving the differences
#Extract pairwise p-values by changing reference level
d$AttLoc=relevel(d$AttLoc, ref="BWE")
mod1=multinom(AttLoc~Form,data=d)
z=summary(mod1)$coefficients/summary(mod1)$standard.errors #calculate z-values
p=(1-pnorm(abs(z),0,1))*2 #calculate p-values
p #prints p-values
d$AttLoc=relevel(d$AttLoc, ref="HWE")
mod1=multinom(AttLoc~Form,data=d)
z=summary(mod1)$coefficients/summary(mod1)$standard.errors #calculate z-values
p=(1-pnorm(abs(z),0,1))*2 #calculate p-values
p #prints p-values
d$AttLoc=relevel(d$AttLoc, ref="FWE")
mod1=multinom(AttLoc~Form,data=d)
z=summary(mod1)$coefficients/summary(mod1)$standard.errors #calculate z-values
p=(1-pnorm(abs(z),0,1))*2 #calculate p-values
p #prints p-values
#Summarised p-values
#	Body	FWE	HWE	BWE
#Body	-	0.475	0.106	0.322
#FWE	-	-	0.208	0.085
#HWE	-	-	-	0.031
##Final results: Only significant difference is the proportion between HWE and BWE; FWE vs BWE is marginally significant

#####Binomial GLMM##### (To check whether inclusion of the random effect would affect the multinomial results)
##Load packages
require(lme4)
require(arm) #for binnedplot and confint

##Pairwise comparisons between all categories of first strikes
#Body vs Forewing Eyespots
dBodyvFWE=d[d$AttLoc=="Body"|d$AttLoc=="FWE",]
mod1=glmer(AttLoc~Form+(1|Mantis),family=binomial,data=dBodyvFWE)
binnedplot(predict(mod1,type="response"),resid(mod1,type="response")) #model checking
summary(mod1) #p=0.475
confint(mod1, method="Wald") #95% confidence intervals
#Body vs Hindwing Eyespots
dBodyvHWE=d[d$AttLoc=="Body"|d$AttLoc=="HWE",]
mod2=glmer(AttLoc~Form+(1|Mantis),family=binomial,data=dBodyvHWE)
binnedplot(predict(mod2,type="response"),resid(mod2,type="response")) #model checking
summary(mod2) #p=0.1336
confint(mod2, method="Wald") #95% confidence intervals
#Body vs Both Wings' Eyespots
dBodyvBWE=d[d$AttLoc=="Body"|d$AttLoc=="BWE",]
mod3=glmer(AttLoc~Form+(1|Mantis),family=binomial,data=dBodyvBWE)
binnedplot(predict(mod3,type="response"),resid(mod3,type="response")) #model checking
summary(mod3) #p=0.324
confint(mod3, method="Wald") #95% confidence intervals
#Forewing Eyespots vs Hindwing Eyespots
dFWEvHWE=d[d$AttLoc=="FWE"|d$AttLoc=="HWE",]
mod4=glmer(AttLoc~Form+(1|Mantis),family=binomial,data=dFWEvHWE)
binnedplot(predict(mod4,type="response"),resid(mod4,type="response")) #model checking
summary(mod4) #p=0.208
confint(mod4, method="Wald") #95% confidence intervals
#Forewing Eyespots vs Both Wings' Eyespots
dFWEvBWE=d[d$AttLoc=="FWE"|d$AttLoc=="BWE",]
mod5=glmer(AttLoc~Form+(1|Mantis),family=binomial,data=dFWEvBWE)
binnedplot(predict(mod5,type="response"),resid(mod5,type="response")) #model checking
summary(mod5) #p=0.092
confint(mod5, method="Wald") #95% confidence intervals
#Hindwing Eyespots vs Both Wings' Eyespots
dHWEvBWE=d[d$AttLoc=="HWE"|d$AttLoc=="BWE",]
mod6=glmer(AttLoc~Form+(1|Mantis),family=binomial,data=dHWEvBWE)
binnedplot(predict(mod6,type="response"),resid(mod6,type="response")) #model checking
summary(mod6) #p=0.0308
confint(mod6, method="Wald") #95% confidence intervals
#Summarised p-values
#	Body	FWE	HWE	BWE
#Body	-	0.475	0.134	0.322
#FWE	-	-	0.208	0.092
#HWE	-	-	-	0.031
#Final result: almost identical to p-values from the multinomial regression above, suggesting that the effect of the mantis is not very important and supporting the conclusions

####################### First Strike Data Ends ######################
#####################################################################

########################################################################################
################################ Arena Experiment Ends #################################
########################################################################################