######################################################################################## ################################# 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 ################################# ########################################################################################