####Analyses and plots carried out for Lopez-Palafox, T., Macias-Ordonez, R. and Cordero C. The size of signal detection and emission organs in a synchronous firefly: sexual dimorphism, allometry and assortative mating


###Load datasets
data<-read.csv(file="S1 database lopezetal.csv", sep=",", header=T); # Full data set for predictions 1-4
Males<-subset(data, Sex=="Male"); #Dataset for males only
Females<-subset(data, Sex=="Female"); #Dataset for females only
dataassort<-read.csv(file="S2 database lopezetal assortative.csv", sep=",", header=T); #Dataset for mating pairs


#### Predictions 1 and 2 (Table 2): lantern, eye or body size vs probability of capture mating or single

##Males

#absolute lantern and eye size
glmM2<-glm(status10~(Body+Lantern+Eyes2)^2, data=Males, family="binomial")
summary(glmM2)##model results before simplification used in Table 2

glmM3<-update(glmM2,~.-Body:Eyes2)#begining of backwards stepwise simplification
summary(glmM3)
anova(glmM2, glmM3, test="Chisq")
glmM4<-update(glmM3,~.-Body:Lantern)
summary(glmM4)
anova(glmM3, glmM4, test="Chisq")
glmM5<-update(glmM4,~.-Lantern:Eyes2)
summary(glmM5)
anova(glmM4, glmM5, test="Chisq")
glmM6<-update(glmM5,~.-Body)
summary(glmM6)
anova(glmM5, glmM6, test="Chisq")
glmM7<-update(glmM6,~.-Eyes2)
summary(glmM7)
anova(glmM6, glmM7, test="Chisq")
glmM8<-update(glmM7,~.-Lantern)
summary(glmM8)
anova(glmM7, glmM8, test="Chisq")#end of backwards stepwise simplification

#relative lantern and eye size

glmMrel2<-glm(status10~(Body+LanternRel+EyesRel)^2, data=Males, family="binomial")
summary(glmMrel2)##model results before simplification used in Table 2

glmMrel3<-update(glmMrel2,~.-Body:EyesRel)#begining of backwards stepwise simplification
summary(glmMrel3)
anova(glmMrel2, glmMrel3, test="Chisq")
glmMrel4<-update(glmMrel3,~.-Body:LanternRel)
summary(glmMrel4)
anova(glmMrel3, glmMrel4, test="Chisq")
glmMrel5<-update(glmMrel4,~.-LanternRel:EyesRel)
summary(glmMrel5)
anova(glmMrel4, glmMrel5, test="Chisq")
glmMrel6<-update(glmMrel5,~.-Body)
summary(glmMrel6)
anova(glmMrel5, glmMrel6, test="Chisq")
glmMrel7<-update(glmMrel6,~.-EyesRel)
summary(glmMrel7)
anova(glmMrel6, glmMrel7, test="Chisq")
glmMrel8<-update(glmMrel7,~.-LanternRel)
summary(glmMrel8)
anova(glmMrel7, glmMrel8, test="Chisq")#end of backwards stepwise simplification
  
##Females

#absolute lantern and eye size

glmF2<-glm(status10~(Body+Lantern+Eyes2)^2, data=Females, family="binomial")
summary(glmF2)##model results before simplification used in Table 2

glmF3<-update(glmF2,~.-Body:Eyes2)#begining of backwards stepwise simplification
summary(glmF3)
anova(glmF2, glmF3, test="Chisq")
glmF4<-update(glmF3,~.-Body:Lantern)
summary(glmF4)
anova(glmF3, glmF4, test="Chisq")
glmF5<-update(glmF4,~.-Lantern:Eyes2)
summary(glmF5)
anova(glmF4, glmF5, test="Chisq")
glmF6<-update(glmF5,~.-Body)
summary(glmF6)
anova(glmF5, glmF6, test="Chisq")
glmF7<-update(glmF6,~.-Eyes2)
summary(glmF7)
anova(glmF6, glmF7, test="Chisq")
glmF8<-update(glmF7,~.-Lantern)
summary(glmF8)
anova(glmF7, glmF8, test="Chisq")#end of backwards stepwise simplification

#relative lantern and eye size using body size
glmFrel2<-glm(status10~(Body+LanternRel+EyesRel)^2, data=Females, family="binomial")
summary(glmFrel2)##model results before simplification used in Table 2

glmFrel3<-update(glmFrel2,~.-Body:EyesRel)#begining of backwards stepwise simplification
summary(glmFrel3)
anova(glmFrel2, glmFrel3, test="Chisq")
glmFrel4<-update(glmFrel3,~.-Body:LanternRel)
summary(glmFrel4)
anova(glmFrel3, glmFrel4, test="Chisq")
glmFrel5<-update(glmFrel4,~.-LanternRel:EyesRel)
summary(glmFrel5)
anova(glmFrel4, glmFrel5, test="Chisq")
glmFrel6<-update(glmFrel5,~.-Body)
summary(glmFrel6)
anova(glmFrel5, glmFrel6, test="Chisq")
glmFrel7<-update(glmFrel6,~.-EyesRel)
summary(glmFrel7)
anova(glmFrel6, glmFrel7, test="Chisq")
glmFrel8<-update(glmFrel7,~.-LanternRel)
summary(glmFrel8)
anova(glmFrel7, glmFrel8, test="Chisq")#end of backwards stepwise simplification

#relative lantern and eye size using elytra size
LanternRelEl<-Females$Lantern/Females$Elytra;EyesRelEl<-Females$Eyes2/Females$Elytra# this line calculates lantern and eye relative size using elytra size since that variable was not calculated in the database

glmFRelEl2<-glm(status10~(Body+LanternRelEl+EyesRelEl)^2, data=Females, family="binomial")
summary(glmFRelEl2)##model results before simplification used in Table 2

glmFRelEl3<-update(glmFRelEl2,~.-Body:EyesRelEl)#begining of backwards stepwise simplification
summary(glmFRelEl3)
anova(glmFRelEl2, glmFRelEl3, test="Chisq")
glmFRelEl4<-update(glmFRelEl3,~.-Body:LanternRelEl)
summary(glmFRelEl4)
anova(glmFRelEl3, glmFRelEl4, test="Chisq")
glmFRelEl5<-update(glmFRelEl4,~.-LanternRelEl:EyesRelEl)
summary(glmFRelEl5)
anova(glmFRelEl4, glmFRelEl5, test="Chisq")
glmFRelEl6<-update(glmFRelEl5,~.-Body)
summary(glmFRelEl6)
anova(glmFRelEl5, glmFRelEl6, test="Chisq")
glmFRelEl7<-update(glmFRelEl6,~.-EyesRelEl)
summary(glmFRelEl7)
anova(glmFRelEl6, glmFRelEl7, test="Chisq")
glmFRelEl8<-update(glmFRelEl7,~.-LanternRelEl)
summary(glmFRelEl8)
anova(glmFRelEl7, glmFRelEl8, test="Chisq")#end of backwards stepwise simplification


####Prediction 3 (Table 3): Sexual dimorphism

summary(Females$Body)
summary(Males$Body)
summary(lm(Body~sex10, data=data))

summary(Females$Elytra)
summary(Males$Elytra)
summary(lm(Elytra~sex10, data=data))

summary(Females$Lantern)
summary(Males$Lantern)
summary(lm(Lantern~sex10, data=data))

summary(Females$LanternRel)*100
summary(Males$LanternRel)*100
summary(lm(LanternRel~sex10, data=data))

summary(Females$Eyes2)
summary(Males$Eyes2)
summary(lm(Eyes2~sex10, data=data))

summary(Females$EyesRel)*100
summary(Males$EyesRel)*100
summary(lm(EyesRel~sex10, data=data))


####Prediction 4 (Table 4): Allometry
## Note: 95% IC intervals for each slope are estimated through a bootstrap randomization process, therefore the outcome may be slightly different in every estimation, and thus they will differ somewhat from those reported in Table 4. The overall conclusions however should not differ based on these slight differences.

###Females

##Lantern vs body
summary(lm(log(Lantern)~log(Body), data=Females))

countmLanternF<- 0; rndmLanternF<- numeric(); mLanternF<-summary(lm(log(Lantern)~log(Body), data=Females))[[4]][[2]]
for(rn in 1:10000) {
  rndmLanternF[rn]<-summary(lm(log(Lantern)~log(Body), data=Females[sample(nrow(Males),replace = T), ]))[[4]][[2]]
  countmLanternF<-countmLanternF + ifelse(rndmLanternF[rn]>mLanternF,1,0)
}
mLanternF; CImLanternF<-quantile(rndmLanternF,c(0.025,0.975)); CImLanternF

##Eyes2 vs body
summary(lm(log(Eyes2)~log(Body), data=Females))

countmEyes2F<- 0; rndmEyes2F<- numeric(); mEyes2F<-summary(lm(log(Eyes2)~log(Body), data=Females))[[4]][[2]]
for(rn in 1:10000) {
  rndmEyes2F[rn]<-summary(lm(log(Eyes2)~log(Body), data=Females[sample(nrow(Males),replace = T), ]))[[4]][[2]]
  countmEyes2F<-countmEyes2F + ifelse(rndmEyes2F[rn]>mEyes2F,1,0)
}
mEyes2F; CImEyes2F<-quantile(rndmEyes2F,c(0.025,0.975)); CImEyes2F

##Lantern vs Elytra
summary(lm(log(Lantern)~log(Elytra), data=Females))

countmLanternFel<- 0; rndmLanternFel<- numeric(); mLanternFel<-summary(lm(log(Lantern)~log(Elytra), data=Females))[[4]][[2]]
for(rn in 1:10000) {
  rndmLanternFel[rn]<-summary(lm(log(Lantern)~log(Elytra), data=Females[sample(nrow(Males),replace = T), ]))[[4]][[2]]
  countmLanternFel<-countmLanternFel + ifelse(rndmLanternFel[rn]>mLanternFel,1,0)
}
mLanternFel; CImLanternFel<-quantile(rndmLanternFel,c(0.025,0.975)); CImLanternFel

##Eyes2 vs Elytra
summary(lm(log(Eyes2)~log(Elytra), data=Females))

countmEyes2Fel<- 0; rndmEyes2Fel<- numeric(); mEyes2Fel<-summary(lm(log(Eyes2)~log(Elytra), data=Females))[[4]][[2]]
for(rn in 1:10000) {
  rndmEyes2Fel[rn]<-summary(lm(log(Eyes2)~log(Elytra), data=Females[sample(nrow(Males),replace = T), ]))[[4]][[2]]
  countmEyes2Fel<-countmEyes2Fel + ifelse(rndmEyes2Fel[rn]>mEyes2Fel,1,0)
}
mEyes2Fel; CImEyes2Fel<-quantile(rndmEyes2Fel,c(0.025,0.975)); CImEyes2Fel

###Males

##Lantern vs body
summary(lm(log(Lantern)~log(Body), data=Males))

countmLanternM<- 0; rndmLanternM<- numeric(); mLanternM<-summary(lm(log(Lantern)~log(Body), data=Males))[[4]][[2]]
for(rn in 1:10000) {
  rndmLanternM[rn]<-summary(lm(log(Lantern)~log(Body), data=Males[sample(nrow(Males),replace = T), ]))[[4]][[2]]
  countmLanternM<-countmLanternM + ifelse(rndmLanternM[rn]>mLanternM,1,0)
}
mLanternM; CImLanternM<-quantile(rndmLanternM,c(0.025,0.975)); CImLanternM

##Eyes2 vs Body
summary(lm(log(Eyes2)~log(Body), data=Males))

countmEyes2M<- 0; rndmEyes2M<- numeric(); mEyes2M<-summary(lm(log(Eyes2)~log(Body), data=Males))[[4]][[2]]
for(rn in 1:10000) {
  rndmEyes2M[rn]<-summary(lm(log(Eyes2)~log(Body), data=Males[sample(nrow(Males),replace = T), ]))[[4]][[2]]
  countmEyes2M<-countmEyes2M + ifelse(rndmEyes2M[rn]>mEyes2M,1,0)
}
mEyes2M; CImEyes2M<-quantile(rndmEyes2M,c(0.025,0.975));CImEyes2M

##Figure 3
Corner_text <- function(text, location="topleft"){
  legend(location,legend=text, bty ="n", pch=NA) 
  }

ba<-bquote("log body area"~"(mm"^"2"*")")
ea<-bquote("log elytra area"~"(mm"^"2"*")")
la<-bquote("log lantern area"~"(mm"^"2"*")")
es<-bquote("log eyes squared"~"(mm"^"2"*")")

par(mfrow=c(3,2))
par(mar=c(4.5,4.5,0.5,1))
plot(log(Lantern)~log(Body), data=Females,xlab=ba,ylab=la,text(3,0.5, "A"), cex=1.5)
abline(lm(log(Lantern)~log(Body), data=Females),col="red", lwd=2)

plot(log(Eyes2)~log(Body), data=Females,xlab=ba,ylab=es,cex=1.5)
abline(lm(log(Eyes2)~log(Body), data=Females),col="grey", lwd=2.5)

plot(log(Lantern)~log(Elytra), data=Females,xlab=ea,ylab=la,cex=1.5)
abline(lm(log(Lantern)~log(Elytra), data=Females),col="red", lwd=2)

plot(log(Eyes2)~log(Elytra), data=Females,xlab=ea,ylab=es,cex=1.5)
abline(lm(log(Eyes2)~log(Elytra), data=Females),col="grey", lwd=2.5)

plot(log(Lantern)~log(Body), data=Males,xlab=ba,ylab=la,cex=1.5)
abline(lm(log(Lantern)~log(Body), data=Males),col="red", lwd=2)

plot(log(Eyes2)~log(Body), data=Males,xlab=ba,ylab=es,cex=1.5)
abline(lm(log(Eyes2)~log(Body), data=Males),col="grey", lwd=2.5)


####Prediction 5 (Table 5): Assortative mating

###Absolute size

 summary(lm(Bodym~Bodyf, data=dataassort))
 
 summary(lm(Elytram~Elytraf, data=dataassort))
 
 summary(lm(Lanternm~Lanternf, data=dataassort))
 
 summary(lm(Eyes2m~Eyes2f, data=dataassort))
 
###Relative to female body

 summary(lm(LanternRelm~LanternRelf, data=dataassort))
 
 summary(lm(EyesRelm~EyesRelf, data=dataassort))

###Relative to female elytra

 summary(lm(LanternRelm~LanternRelElf, data=dataassort))

 summary(lm(EyesRelm~EyesRelElf, data=dataassort))
 
###Figure 4
 par(mfrow=c(1,1))
 plot(LanternRelm~LanternRelElf, data=dataassort,xlab="Female lantern relative size",ylab="Male lantern relative size",cex=1.5)
 abline(lm(LanternRelm~LanternRelElf, data=dataassort),col="red", lwd=2)