--- title: "Sexual selection does not drive hindwing tail elaboration in a moon moth, Actias luna - Supplementary Information" author: "Code compiled by: Juliette J. Rubin" date: "1/15/2023" output: pdf_document --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = F, warning=FALSE, message=FALSE,tidy.opts = list(width.cutoff=60),tidy=T) library(lme4) library(tidyr) library(dplyr) library(emmeans) library(car) library(DHARMa) library(knitr) library(MuMIn) library(sjPlot) library(broom) ``` # Defining the parameters Here are definitions for the parameters in all models: Date = the date on which trial was run Cage = the cage identity in which the trial was run. This does not necessarily correlate to the cage size or moth treatments, as the cages were given new identifiers across the data collection seasons based on cage availability, number of moths available to participate in the experiment, etc. Instead, this parameter, in conjunction with the date, can be used to identify which males (rows) were pit against each other in the mating competition. Tent_type = the model of cage in which trial took place. BA = Big Agnes Tiger Wall 3, Eureka = Eureka Suma, BioQuip = BioQuip collapsible cage. Cage_category = The size of the tent in which trials were conducted (large vs small -- large tents were either Eureka Suma + or Big Agnes Tiger Wall 3 (~180 x 230 x 110cm), small tent was BioQuip (~60 x 60 x 90cm)). Temp = Temperature reading at a single time point ~3 hours after trial start (units = degree F). Hum = Humidity reading at a single time point ~3 hours after trial start (units = %). SQM = An integrated measurement of the nocturnal light environment taken with a Sky Quality Meter (larger numbers are indicative of darker environments). Measurement taken at a single time point ~3 hours after trial start (units = magnitudes/arcsecond). Size_Rwings = The sum of the surface area of the right fore- and hindwings of each male moth, measured in imageJ. Used as a proxy for moth size (units = cm^2). Treatment = what experimental treatment the moth was assigned Success = binomial code for whether the moth mated with the female or not (in each) dyad = an identifer for the experimental set to which the moth belonged, as described below. Dyad 1 = intact v. ablated males, dyad 2 = sham control v. ablated control males, dyad 3 = intact v. glue control males. Note: We scaled the numerical variables in the models to allow for comparison across variable type ```{r data} #Reading in the full data set data<-read.csv("Rubin&Kawahara2023_BehavEcol_MatingData.csv", strip.white=T) #making some potentially necessary modifications to variable class data$Cage_category<-as.factor(data$Cage_category) data$Treatment<-as.factor(data$Treatment) data$Tent_type<-as.factor(data$Tent_type) ``` # Defining the experimental sets and building models Because each dyad was run as a separate experiment, we are dividing the data set up into the dyads that were introduced to females and exploring them separately. The dyads are as follows: Intact vs Ablated: male with tails intact vs. male with tails removed Sham control vs Ablated control: male with tails removed and reglued vs. male with tails removed and glue added to the base of the hindwings Intact vs glue control: male with tails intact vs. male with tails intact and glue added to the base of the hindwings. Glue was applied to the underside of the moth wings whenever it was used (see full paper for more detail). ```{r subseting_dyads} ##creating the subsetted data sets #Intact vs Ablated only IvA<-subset(data, Treatment !="Sham_control"& Treatment !="Ablated_control" & Treatment !="Intact_glue" & Treatment !="Intact2") #Sham vs Ablated control only SvAC<-subset(data, Treatment !="Intact"& Treatment !="Ablated" & Treatment !="Intact_glue" & Treatment !="Intact2") #Intact vs Glue control only IvGC<-subset(data, Treatment !="Sham_control"& Treatment !="Ablated_control" & Treatment !="Intact" & Treatment !="Ablated") ``` ## Intact vs. Ablated We start with the primary test between a male with tails and a male with tails removed. ```{r Intact-ablated_models,echo=TRUE} #model set mod1<-glm(Success~Treatment,data=IvA,family="binomial") mod2<-glm(Success~Treatment+scale(Size_Rwings),data=IvA,family="binomial") mod3<-glm(Success~Treatment+scale(SQM),data=IvA,family="binomial") mod4<-glm(Success~Treatment+scale(Temp),data=IvA,family="binomial") mod5<-glm(Success~Treatment+scale(Hum),data=IvA,family="binomial") mod6<-glm(Success~Treatment+Cage_category,data=IvA,family="binomial") #Here is the model selection table model.sel(mod1,mod2,mod3,mod4,mod5,mod6,rank=BIC) ``` Model 2 has the lowest BIC score, so we will proceed with interpreting this model. \newpage We will first perform model checks to ensure this model does not demonstrate over or under-dispersion. We will then extract parameter estimates and p-values. ```{r outputs_intact-ablated} #first checking that treatment and size are not multicollinear; they are not vif1<-vif(mod2) #checking for over or under-dispersion check_mod2 <- simulateResiduals(fittedModel = mod2,n=1000, plot=T) #n.s. kable(tidy(mod2),caption = "Table output for binomial model comparing mating success between A. luna males with tails intact and tails removed, accounting for size of each male", col.names = c("Predictor", "parameter estimate", "SE", "t", "p"), digits = c(0, 2, 3, 2, 3)) #summary statistics and model output tables summ_mod2<-emmeans(mod2, "Treatment", type="response", level=0.95) cat("here are the parameter estimates, back-transformed from the logit scale") print(summ_mod2) ``` \newpage Estimating the intercept -- when males are the same size, do females still show a preference? ```{r size_mod_intact-ablated} #restructing df to wide format newIvA<-IvA %>% pivot_wider( names_from = Treatment, values_from = c(Success, Size_Rwings) ) #Finding the difference between intact and ablated moths newIvA$size_diff<-newIvA$Size_Rwings_Intact - newIvA$Size_Rwings_Ablated #building the intercept model size_mod<-glm(Success_Intact~size_diff,data=newIvA,family="binomial") #Here is the summary for the size intercept model summary(size_mod) ``` So, difference in male sizes does not have a significant effect on trial outcome \newpage ## Sham control vs. Ablated control We then conducted additional trials to test for the effect of removing the tail - comparing the mating success of males whose tails we removed and reglued and males whose tails we removed and glue only added to the base of the hindwing. ```{r damage_control_models, echo=TRUE} #Sham vs ablated control models #model set mod7<-glm(Success~Treatment,data=SvAC,family="binomial") mod8<-glm(Success~Treatment+scale(Size_Rwings),data=SvAC,family="binomial") mod9<-glm(Success~Treatment+scale(SQM),data=SvAC,family="binomial") mod10<-glm(Success~Treatment+scale(Temp),data=SvAC,family="binomial") mod11<-glm(Success~Treatment+scale(Hum),data=SvAC,family="binomial") mod12<-glm(Success~Treatment+Tent_type,data=SvAC,family="binomial") #Here is the model selection table model.sel(mod7,mod8,mod9,mod10,mod11,mod6,mod12,rank=BIC) ``` mod7 and mod8 are within 2 BIC of each other. However, mod7 has a simpler structure (treatment as the sole predictor versus treatment + size) and thus, given this limited sample size, we will stick with this simpler model. We will first present to model checks for mod7 and then proceed to extracting the summary statistics. ```{r model_output_damage-control} #checking for over or under-dispersion check_mod7 <- simulateResiduals(fittedModel = mod7,n=1000, plot=T) #looks good kable(tidy(mod7),caption = "Table output for binomial model comparing mating success between A. luna males with tails removed and reglued (sham control) and A. luna males with tails removed and glue only added (ablated control)", col.names = c("Predictor", "parameter estimate", "SE", "t", "p"), digits = c(0, 2, 3, 2, 3)) #summary statistics and model output tables summ_mod7<-emmeans(mod7, "Treatment", type="response", level=0.95) cat("here are the parameter estimates, back-transformed from the logit scale") print(summ_mod7) ``` \newpage Estimating the intercept -- when males are the same size, do females still show a preference? ```{r size_model-damage_controls} #restructing df to wide format newSvAC<-SvAC %>% pivot_wider( names_from = Treatment, values_from = c(Success, Size_Rwings) ) #Finding the difference between intact and ablated moths newSvAC$size_diff<-newSvAC$Size_Rwings_Sham_control - newSvAC$Size_Rwings_Ablated_control #building the intercept model size_mod2<-glm(Success_Sham_control~size_diff,data=newSvAC,family="binomial") #Here is the summary of the size intercept model summary(size_mod2) ``` So again, difference in male sizes does not have a significant effect on trial outcome ## Intact vs. glue control To control for any sensory masking effect that the glue may have had in the damage control trials, we compared the mating success between males with intact hindwings and males with intact hindwings with glue added. Note: due to tent availability, the smaller tents (BioQuip) were not used for this dyad. Therefore, the Cage_category parameter has been removed. ```{r glue_control_models, echo=TRUE} # model set mod13<-glm(Success~Treatment,data=IvGC,family="binomial") mod14<-glm(Success~Treatment+scale(Size_Rwings),data=IvGC,family="binomial") mod15<-glm(Success~Treatment+scale(Temp),data=IvGC,family="binomial") mod16<-glm(Success~Treatment+scale(Hum),data=IvGC,family="binomial") #Here is the model selection table model.sel(mod13,mod14,mod15,mod16,rank=BIC) ``` Here again mod13 and 14, which include Treatment alone and Treatment + size, perform similarly.Given that adding the size term did not significantly improve model fit, we will provide model checks and interpretation for the simplest model only. ```{r model_outputs-glue_controls} #checking for over or under-dispersion check_mod13 <- simulateResiduals(fittedModel = mod13,n=1000, plot=T) #looks good kable(tidy(mod13),caption = "Table output for binomial model comparing mating success between A. luna males with intact tails and A. luna males with glue added to intact tails", col.names = c("Predictor", "parameter estimate", "SE", "t", "p"), digits = c(0, 2, 3, 2, 3)) #summary statistics and model output tables summ_mod13<-emmeans(mod13, "Treatment", type="response", level=0.95) cat("here are the parameter estimates, back-transformed from the logit scale") print(summ_mod13) ``` \newpage Estimating the intercept -- when males are the same size, do females still show a preference? ```{r size_model-glue_control} #restructing df to wide format newIvGC<-IvGC %>% pivot_wider( names_from = Treatment, values_from = c(Success, Size_Rwings) ) #Finding the difference between intact and ablated moths newIvGC$size_diff<-newIvGC$Size_Rwings_Intact2 - newIvGC$Size_Rwings_Intact_glue #building the intercept model size_mod3<-glm(Success_Intact2~size_diff,data=newIvGC,family="binomial") #Here is a summary of the size intercept model summary(size_mod3) ``` Again, difference in male sizes does not have a significant effect on trial outcome ## Comparing all experimental sets To get a sense of whether the outcome of each dyad statistically differed from each other, we compared the mating success of all treatments within one model. Here, we use dyad ID as a random intercept, to account for the grouping arrangements of the treatments. ```{r all_data,fig.cap="Comparing mating success across all dyads. Geoms are depicted in the true pairs (i.e., Intact - Glue control; Intact - Ablated; Sham control - Ablated control. Central points are marginal estimates and error bars are 95% confidence intervals"} #reordering treatments for better comparisons & visualization data2 <- mutate(data, Treatment = ordered(Treatment, c("Intact2", "Intact_glue","Intact", "Ablated","Sham_control","Ablated_control"))) all_data<-glmer(Success~0+Treatment+(1|dyad),data=data2,family="binomial") cat("here is the model summary, with p-values included. You can see that only the intact & ablated dyad vary from the ~50% mean of the other dyads") summary(all_data) cat("here are the marginal estimates, back-transformed from the logit scale") emmeans(all_data, "Treatment", type = "response") plot_model(all_data,type="pred",pred.style="eff",ci.lvl = 0.95, transform="plogis") ``` So we find that the proportional differences in mating success between the intact and ablated dyad did differ from the outcomes of the other dyads -- that is, only in this group was the mating sucess rate different from ~50% for each of the treatments. \newpage # Analyzing video recorded behavior We video recorded a subset of mating trials and here use an exact binomial test to measure whether the male that bopepd the female more was more successful in mating with her than the male that bopped the female less over the course of an experimental night. ```{r behavior} #reading in the data behav_dat<-read.csv("Rubin&Kawahara_video-data.csv") #modifying the class of the variables behav_dat$M_more_bops_mated<-as.factor(behav_dat$M_more_bops_mated) behav_dat$Num_bops_M1<-as.numeric(behav_dat$Num_bops_M1) #exact binomial test, testing whether bopping led to more successful mating #in 11 of 18 recorded trials the male that bopped the female got the mating binom.test(x=11, n=18, p=.5) ``` p-value > 0.05. We therefore do not find a significant difference in bopping behavior between males who succeeded in mating with the female and males who failed. **Please see the .rmd file for the full code & outputs. For ease of reading, the code chunks have been hidden from this pdf document.