#R code - summary
#Evolutionary loss of a signalling colour is linked to increased response to conspecific chemicals
#C. Romero-Diaz*, J. A. Pruett, S. M. Campos, A. G. Ossip-Drahos, J. J. Zúñiga-Vega, C. Vital-García, D. K. Hews & E. P. Martins. *cromerod@asu.edu
#Proceedings of the Royal Society B. 2021

 ##Overall models
#Chemical acts:
fit<-glm.nb(Total_chemo~Clade * Belly * Trt + offset(log(Trial_Dur)), data=sp11)
drop1(fit, test="Chi")
summary(fit); confint(fit)
Anova(fit,type = "II")
testInteractions(fit, fixed="Clade", "Belly")
testInteractions(fit, fixed="Belly", "Clade")
testInteractions(fit, fixed="Trt", "Belly")
testInteractions(fit, fixed="Belly", "Trt")
testInteractions(fit, fixed="Clade", "Trt")
testInteractions(fit, fixed="Trt", "Clade")
#Visual displays:
fit<-glm.nb(Total_Disp~Clade * Belly * Trt + offset(log(Trial_Dur)), data=sp11)
summary(fit); confint(fit)
drop1(fit, test="Chi")
Anova(fit,type = "II")
testInteractions(fit, fixed="Clade", "Trt")
testInteractions(fit, fixed="Trt", "Clade")
testInteractions(fit, fixed="Clade", "Belly")
testInteractions(fit, fixed="Belly", "Clade")

with(fit, cbind(res.deviance = deviance, df = df.residual, p = pchisq(deviance, df.residual, lower.tail=FALSE)))
exp(coef(fit))  
fit$theta

 ##Species models
split.Sp<-split(sp11, sp11$Species)
#Chemical acts:
fit<-glm.nb(Total_chemo~Trt+ offset(log(Trial_Dur)), data=split.Sp$cozumelae); summary(fit)	#NB GLM
dunn.test(split.Sp$cozumelae$chem_rate5,split.Sp$cozumelae$Trt, method="holm")			#Non-parametric
fit<-glm.nb(Total_chemo~Trt+ offset(log(Trial_Dur)), data=split.Sp$graciosus); summary(fit)
dunn.test(split.Sp$graciosus$chem_rate5,split.Sp$graciosus$Trt, method="holm")
fit<-glm.nb(Total_chemo~Trt+ offset(log(Trial_Dur)), data=split.Sp$grammicus); summary(fit)
dunn.test(split.Sp$grammicus$chem_rate5,split.Sp$grammicus$Trt, method="holm")
fit<-glm.nb(Total_chemo~Trt+ offset(log(Trial_Dur)), data=split.Sp$jarrovii); summary(fit)
dunn.test(split.Sp$jarrovii$chem_rate5,split.Sp$jarrovii$Trt, method="holm")
fit<-glm.nb(Total_chemo~Trt+ offset(log(Trial_Dur)), data=split.Sp$megalepidurus); summary(fit)
dunn.test(split.Sp$megalepidurus$chem_rate5,split.Sp$megalepidurus$Trt, method="holm")
fit<-glm.nb(Total_chemo~Trt+ offset(log(Trial_Dur)), data=split.Sp$merriami); summary(fit)
dunn.test(split.Sp$merriami$chem_rate5,split.Sp$merriami$Trt, method="holm")
fit<-glm.nb(Total_chemo~Trt+ offset(log(Trial_Dur)), data=split.Sp$occidentalis); summary(fit)
dunn.test(split.Sp$occidentalis$chem_rate5,split.Sp$occidentalis$Trt, method="holm")
fit<-glm.nb(Total_chemo~Trt+ offset(log(Trial_Dur)), data=split.Sp$parvus); summary(fit)
dunn.test(split.Sp$parvus$chem_rate5,split.Sp$parvus$Trt, method="holm")
fit<-glm.nb(Total_chemo~Trt+ offset(log(Trial_Dur)), data=split.Sp$siniferus); summary(fit)
dunn.test(split.Sp$siniferus$chem_rate5,split.Sp$siniferus$Trt, method="holm")
fit<-glm.nb(Total_chemo~Trt+ offset(log(Trial_Dur)), data=split.Sp$variabilis); summary(fit)
dunn.test(split.Sp$variabilis$chem_rate5,split.Sp$variabilis$Trt, method="holm")
fit<-glm.nb(Total_chemo~Trt+ offset(log(Trial_Dur)), data=split.Sp$virgatus); summary(fit)
dunn.test(split.Sp$virgatus$chem_rate5,split.Sp$virgatus$Trt, method="holm")

#Visual displays:
fit<-glm.nb(Total_Disp~Trt+ offset(log(Trial_Dur)), data=split.Sp$cozumelae); summary(fit)
dunn.test(split.Sp$cozumelae$Disp_rate5,split.Sp$cozumelae$Trt, method="holm")
fit<-glm.nb(Total_Disp~Trt+ offset(log(Trial_Dur)), data=split.Sp$graciosus); summary(fit)
dunn.test(split.Sp$graciosus$Disp_rate5,split.Sp$graciosus$Trt, method="holm")
fit<-glm.nb(Total_Disp~Trt+ offset(log(Trial_Dur)), data=split.Sp$grammicus); summary(fit)
dunn.test(split.Sp$grammicus$Disp_rate5,split.Sp$grammicus$Trt, method="holm")
fit<-glm.nb(Total_Disp~Trt+ offset(log(Trial_Dur)), data=split.Sp$jarrovii); summary(fit)
dunn.test(split.Sp$jarrovii$Disp_rate5,split.Sp$jarrovii$Trt, method="holm")
fit<-glm.nb(Total_Disp~Trt+ offset(log(Trial_Dur)), data=split.Sp$megalepidurus); summary(fit)
dunn.test(split.Sp$megalepidurus$Disp_rate5,split.Sp$megalepidurus$Trt, method="holm")
fit<-glm.nb(Total_Disp~Trt+ offset(log(Trial_Dur)), data=split.Sp$merriami); summary(fit)
dunn.test(split.Sp$merriami$Disp_rate5,split.Sp$merriami$Trt, method="holm")
fit<-glm.nb(Total_Disp~Trt+ offset(log(Trial_Dur)), data=split.Sp$occidentalis); summary(fit)
dunn.test(split.Sp$occidentalis$Disp_rate5,split.Sp$occidentalis$Trt, method="holm")
fit<-glm.nb(Total_Disp~Trt+ offset(log(Trial_Dur)), data=split.Sp$parvus); summary(fit)
dunn.test(split.Sp$parvus$Disp_rate5,split.Sp$parvus$Trt, method="holm")
fit<-glm.nb(Total_Disp~Trt+ offset(log(Trial_Dur)), data=split.Sp$siniferus); summary(fit)
dunn.test(split.Sp$siniferus$Disp_rate5,split.Sp$siniferus$Trt, method="holm")
fit<-glm.nb(Total_Disp~Trt+ offset(log(Trial_Dur)), data=split.Sp$variabilis); summary(fit)
dunn.test(split.Sp$variabilis$Disp_rate5,split.Sp$variabilis$Trt, method="holm")
fit<-glm.nb(Total_Disp~Trt+ offset(log(Trial_Dur)), data=split.Sp$virgatus); summary(fit)
dunn.test(split.Sp$virgatus$Disp_rate5,split.Sp$virgatus$Trt, method="holm")

#Post-hocs
summary(glht(fit,linfct=mcp(Trt="Tukey")), test=adjusted(type="holm"))
#(E.g. planned contrasts)#contrast.matrix <- rbind("con vs. STI" =  c(0, 0, 0, 1),"cue vs. STI" =  c(0, 1, 0, -1), "baseline vs. STI" = c(0,0,1,-1))
#summary(glht(fit, contrast.matrix),test=adjusted(type="holm"))

#[Vis Displays-Chem acts] difference - (Figure S3)
sp11$disp_chem<- sp11$Disp_rate5-sp11$chem_rate5
 #Vis:
subset<-sp11[sp11$Trt %in% c("STI", "Baseline"),]; subset<-droplevels(subset)
fit<-glm(disp_chem ~ Belly * Trt, data=subset)
drop1(fit, test="Chi")
 #Chem:
subset<-sp11[sp11$Trt %in% c("con", "cue"),]; subset<-droplevels(subset)
fit2<-glm(disp_chem ~ Belly * Trt, data=subset)
drop1(fit2, test="Chi")