---
output:
  word_document: default
  html_document: default
---
F---
title: "All analysis_Impact of irradiation and transportation on G. morsitans morsitans"
author: "Mirieri et al.,"
date: "12-06-2024"
output: word_document
---


Load library


```{r}  

library(datasets)
library(gcookbook)
library(ggplot2)
library(plyr)
library(dplyr)
library(lattice)
library(MASS)
library(rcompanion)
library(survival)
library(ranger)
library(ggfortify)
library(rmarkdown)
library(knitr)
library(coxme)
## Graphics
library(lattice)
library(tidyverse)
library(gapminder)
library(FSA)
library(stats)
library(RCA)
#library(broom)
library(sp)
library(ggpubr)
library(AICcmodavg)
library(car)
library(ggthemes)
## Mixed generalized linear models
library(lme4)
library(MuMIn) 
library(nlme)
library(survminer)
```

Working directory
```{r setup, include=FALSE, echo=FALSE}

require("knitr")
opts_knit$set(root.dir ="D:/2024/Final folder_irrtrans_06022024/Raw data_14022022_CMK/June 2024")
```


ANALYSIS OF THE IMPACT OF IRRADIATION AND TRASPORATION(22-day old only/29-day old only)


MAIN FIGURES (MANUSCRIPT FIGURES)



Figure 1:Emergence rate(22 and 29 days)-differences between ages within each treatment


```{r}

tab1= read.csv("Figure 1.csv")
head(tab1)
Emerg_rate<- tab1$emerged / (tab1$unemerged+ tab1$emerged)
tab1$Pupal_age<- as.factor(tab1$Pupal_age)

            
Figure_1<-ggplot(tab1, aes(x=factor(Treatments),y=Emerg_rate, colour=Pupal_age)) +
  geom_boxplot(position=position_dodge(0.8))+
  geom_jitter(position=position_dodge(0.8))
              
              
Figure_1

tiff("Figure 1.tiff", width = 7, height = 4, units = 'in',compression = 'lzw',res = 300)
plot(Figure_1+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + theme(axis.text.x = element_text(colour = "black")) + theme(axis.text.y = element_text(colour = "black")) + xlab(expression(bold("Treatments"))) + ylab(expression(bold("Emergence rate")))
dev.off()


#### Significance of treatments  


##22 days 

Age22_fig1 <- subset(tab1, Pupal_age== "22")


head(Age22_fig1)
Age22_fig1$rate <- Age22_fig1$emerged/ (Age22_fig1$unemerged+ Age22_fig1$emerged)
fm1_22 <- glmer(cbind(emerged, unemerged) ~ Treatment +(1|Replicate), family = binomial, data =Age22_fig1)
summary(fm1_22)



##29 days 

Age29_fig1 <- subset(tab1, Pupal_age== "29")

head(Age29_fig1 )
Age29_fig1 $rate <- Age29_fig1 $emerged/ (Age29_fig1 $unemerged+ Age29_fig1 $emerged)                                                                                                                                                           
fm1_29 <- glmer(cbind(emerged, unemerged) ~ Treatment +(1|Replicate), family = binomial, data = Age29_fig1 )
summary(fm1_29)


```


Figure 2:Flight propensity(22 and 29 days)-differences between ages within each treatment


```{r}
##"flight propensity combined data 
tab2= read.csv("Figure 2.csv")
head(tab2)
Flight_rate<- tab2$out / (tab2$in.+ tab2$out)
tab2$Pupal_age<- as.factor(tab2$Pupal_age)

Figure_2<-ggplot(tab2, aes(x=Treatments,y=Flight_rate, colour=Pupal_age)) +
  geom_boxplot(position=position_dodge(0.8))+
  geom_jitter(position=position_dodge(0.8))
Figure_2

tiff("Figure 2.tiff", width = 7, height = 4, units = 'in',compression = 'lzw',res = 300)
plot(Figure_2+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + theme(axis.text.x = element_text(colour = "black")) + theme(axis.text.y = element_text(colour = "black")) + xlab(expression(bold("Treatments"))) + ylab(expression(bold("Flight propensity")))
dev.off()

###### Significance of treatments 
##22 days 


Age22_fig2 <- subset(tab2, Pupal_age== "22")

head(Age22_fig2 )
Age22_fig2 $rate <- Age22_fig2 $out / (Age22_fig2 $in. + Age22_fig2 $out)
fm2_22 <- glmer(cbind(out,in.) ~ Treatment +(1|Replicate), family = binomial, data = Age22_fig2 )
summary(fm2_22)



##29 days 

Age29_fig2 <- subset(tab2, Pupal_age== "29")

head(Age29_fig2)
Age29_fig2$rate <- Age29_fig2$out / (Age29_fig2$in. + Age29_fig2$out)
fm2_29 <- glmer(cbind(out,in.) ~ Treatment +(1|Replicate), family = binomial, data = Age29_fig2)
summary(fm2_29)


```


Figure 3: Mating ability(22 and 29 days)-differences between ages within each treatment



```{r }

#####################-combined data###Figure 3

tab3= read.csv("Figure 3.csv")
head(tab3)
mating_prop<- tab3$pairs_formed  / (tab3$unformed_pairs+ tab3$pairs_formed )
tab3$Pupal_age<- as.factor(tab3$Pupal_age)

Figure_3<-ggplot(tab3, aes(x=Treatments,y=mating_prop, colour=Pupal_age)) +
  geom_boxplot(position=position_dodge(0.8))+
  geom_jitter(position=position_dodge(0.8))
Figure_3

tiff("Figure 3.tiff", width = 7, height = 4, units = 'in',compression = 'lzw', res = 300)
plot(Figure_3+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + theme(axis.text.x = element_text(colour = "black")) + theme(axis.text.y = element_text(colour = "black")) + xlab(expression(bold("Treatments"))) + ylab(expression(bold("Mating ability")))
dev.off()



#### Significance of treatments  

##22 days 


Age22_fig3 <- subset(tab3, Pupal_age== "22")

head(Age22_fig3 )
mating_prop <- Age22_fig3 $pairs_formed/ (Age22_fig3 $unformed_pairs+ Age22_fig3 $pairs_formed)
fm3_22<- glmer(cbind(pairs_formed, unformed_pairs) ~ Treatment +(1|Replicate), family = binomial, data = Age22_fig3 )
summary(fm3_22)



#29 days

Age29_fig3 <- subset(tab3, Pupal_age== "29")


head(Age29_fig3 )
mating_prop <- Age29_fig3 $pairs_formed/ (Age29_fig3 $unformed_pairs+ Age29_fig3 $pairs_formed)
fm3_29 <- glmer(cbind(pairs_formed, unformed_pairs) ~ Treatment +(1|Replicate), family = binomial, data = Age29_fig3 )
summary(fm3_29)


###### Significance of treatments 



```


Figure 4: Insemination rate(22 and 29 days)-differences between ages within each treatment


```{r }
tab4= read.csv("Figure 4.csv")
head(tab4)
Insemination_rate<- tab4$Inseminated / (tab4$Empty+ tab4$Inseminated)
tab4$Pupal_age<- as.factor(tab4$Pupal_age)

Figure_4<-ggplot(tab4, aes(x=Treatments,y=Insemination_rate, colour=Pupal_age)) +
  geom_boxplot(position=position_dodge(0.8))+
  geom_jitter(position=position_dodge(0.8))
Figure_4

tiff("Figure 4.tiff", width = 7, height = 4, units = 'in',compression = 'lzw', res = 300)
plot(Figure_4+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + theme(axis.text.x = element_text(colour = "black")) + theme(axis.text.y = element_text(colour = "black")) + xlab(expression(bold("Treatments"))) + ylab(expression(bold("Insemination rate")))
dev.off()



#### Significance of treatments  

##22 days 

Age22_fig4 <- subset(tab4, Pupal_age== "22")

head(Age22_fig4)
Insemination_rate<- Age22_fig4$Inseminated / (Age22_fig4$Empty+ Age22_fig4$Inseminated)
fm4_22 <- glmer(cbind(Inseminated,Empty) ~ Treatment +(1|Replicate), family = binomial, data = Age22_fig4)
summary(fm4_22 )


##29 days 

Age29_fig4  <- subset(tab4, Pupal_age== "29")

head(Age29_fig4 )
Insemination_rate<- Age29_fig4 $Inseminated / (Age29_fig4 $Empty+ Age29_fig4 $Inseminated)
fm4_29 <- glmer(cbind(Inseminated,Empty) ~ Treatment +(1|Replicate), family = binomial, data = Age29_fig4 )
summary(fm4_29)



###### Significance of treatments 

```




Figure 5:Mean spermathecal value(22 and 29 days) differences between ages within each treatment


```{r }

tab5= read.csv("Figure 5.csv")
head(tab5)
tab5$Pupal_age<- as.factor(tab5$Pupal_age)

Figure_5<-ggplot(tab5, aes(x=Treatments,y=MSV, colour=Pupal_age)) +
  geom_boxplot(position=position_dodge(0.8))+
  geom_jitter(position=position_dodge(0.8))
Figure_5

tiff("Figure 5.tiff", width = 7, height = 4, units = 'in',compression = 'lzw', res = 300)
plot(Figure_5 +theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + theme(axis.text.x = element_text(colour = "black")) + theme(axis.text.y = element_text(colour = "black")) + xlab(expression(bold("Treatments"))) + ylab(expression(bold("Mean Spermathecal Value (MSV)")))
dev.off()

###### Significance of treatments 

##22 days 


Age22_fig5 <- subset(tab5, Pupal_age== "22")

head(Age22_fig5)
fm5_22<- lme(MSV ~ Treatment, random=~1|Replicate,data = Age22_fig5)
summary(fm5_22)


##29days 

Age29_fig5  <- subset(tab5, Pupal_age== "29")

head(Age29_fig5)
fmS26<- lme(MSV ~ Treatments, random=~1|Replicate,data = Age29_fig5)
summary(fmS26)

```

Figure 6a:Spermathecal fill distribution- 22 days-Chi square values 


```{r }


#first,full dataset
tab <- read.csv("spermfill_22.csv", sep=",")
summary(tab)
head(tab)
chisq.test(tab[,5:6])


tab <- read.csv("spermfill_22.csv", sep=",")
summary(tab)
head(tab)

chisq.test(tab[,2:6])


tab<- read.csv("spermfill_22_ship_110Gy.csv", sep=",")
summary(tab)
head(tab)
chisq.test(tab[,2:6])

tab <- read.csv("spermfill_22_Shipped-0Gy .csv", sep=",")
summary(tab)
head(tab)
chisq.test(tab[,2:6])

tab <- read.csv("spermfill_22_Unshipped-0Gy.csv", sep=",")
summary(tab)
head(tab)
chisq.test(tab[,2:6])


tab <- read.csv("spermfill_22_Unshipped-110Gy.csv", sep=",")
summary(tab)
head(tab)
chisq.test(tab[,2:6])

```

Figure 6b: Spermathecal fill distribution- Chi square values 


```{r }
#first,full dataset
tab <- read.csv("Spermfill_29.csv", sep=",")
summary(tab)
head(tab)
chisq.test(tab[,5:6])

tab <- read.csv("Spermfill_29_Ship-0Gy.csv", sep=",")
summary(tab)
head(tab)
chisq.test(tab[,2:6])

tab <- read.csv("Spermfill_29_Ship-110Gy.csv", sep=",")
summary(tab)
head(tab)
chisq.test(tab[,2:6])

tab <- read.csv("Spermfill_29_Unship-110Gy.csv", sep=",")
summary(tab)
head(tab)
chisq.test(tab[,2:6])

tab <- read.csv("Spermfill_29_Unshipped-0Gy.csv", sep=",")
summary(tab)
head(tab)
chisq.test(tab[,2:6])

```



Kaplan Meier Analysis (km)on the impact of treatments on the duration of survival of blood fed flies(trends)

Figure 7a and 7b


```{r }

##Figure 7a:Comparing difference of trends between treatments(22 days old )

data_2=read.csv("Figure 7a.csv")

#data_2
str(data_2)
head(data_2)
data_2=na.omit(data_2)
names(data_2)


km <- with(data_2, Surv(day, status))
head(km,115)
km_trt_fit <- survfit(Surv(day, status) ~ Treatment, data=data_2)
plot(survfit(Surv(day,status)~Treatments,data=data_2), xlab = "Time(days)", ylab = "Survival rate (22 day-old)", col=c('red','blue','green','orange'), lwd=2, xlim =c(0, 115))
legend('topright', text.font =, cex=1, c("Control","Shipped", "Irradiated","shipped_irradiated"), col=c('red','blue', 'green','orange'), lty = 1, lwd=2, box.lty = 1)


## survival: median Kaplan-Meier estimator 
survfit(Surv(day,status)~Treatment,data=data_2)
#data_2$Treatment <- relevel(data_2$Treatment, ref = "Unshipped-0Gy")
cox.surv <- coxph(Surv(day, status) ~ Treatment, data = data_2, method ="exact")
summary(cox.surv)
```


Figure 7b:Comparing difference of trends between treatments(29 day old )

```{r }
data_3=read.csv("Figure 7b.csv")

#data_3
str(data_3)
head(data_3)
data_3=na.omit(data_3)
names(data_3)


km <- with(data_3, Surv(day, status))
head(km,115)
km_trt_fit <- survfit(Surv(day, status) ~ Treatment, data=data_3)
plot(survfit(Surv(day,status)~Treatments,data=data_3), xlab = "Time(days)", ylab = "Survival rate (29 day-old)", col=c('red','blue','green','orange'), lwd=2, xlim =c(0, 115))
legend('topright', text.font = , cex=1, c("Control","Shipped", "Irradiated","shipped and irradiated"), col=c('red','blue', 'green','orange'), lty = 1, lwd=2, box.lty = 1)



### survival: median Kaplan-Meier estimator 
survfit(Surv(day,status)~Treatment,data=data_3)
cox.surv <- coxph(Surv(day, status) ~ Treatment, data = data_3, method ="exact")
summary(cox.surv)

```



Figure 7c and 7d: Survival trends- Kaplan Meier Analysis (km)



```{r }




##Figure 7c: Comparing difference of trends between treatments(combined data)

data_1=read.csv("Figure 7c.csv")

#data_1
str(data_1)
head(data_1)
data_1=na.omit(data_1)
names(data_1)

#Kaplan Meier Analysis (km)
km <- with(data_1, Surv(day, status))
head(km,115)
km_trt_fit <- survfit(Surv(day, status) ~ Treatments, data=data_1)
plot(survfit(Surv(day,status)~Treatments,data=data_1), xlab = "Time(days)", ylab = "Survival rate(22 and 29 day-old)", col=c('red','blue','green','orange'), lwd=2, xlim =c(0, 115))
legend('topright', text.font = , cex=1, c("Control","Shipped", "Irradiated","shipped and irradiated"), col=c('red','blue', 'green','orange'), lty = 1, lwd=2, box.lty = 1)


### survival: median Kaplan-Meier estimator 
survfit(Surv(day,status)~Treatment,data=data_1)
cox.surv <- coxph(Surv(day, status) ~ Treatment, data = data_1, method ="exact")
summary(cox.surv)


####Figure 7d: Comparing difference of trends between ages(combined data)


data_1=read.csv("Figure 7d.csv")
km <- with(data_1, Surv(day, status))
head(km,115)
km_trt_fit <- survfit(Surv(day, status) ~ Pupal_age, data=data_1)
plot(survfit(Surv(day,status)~Pupal_age,data=data_1), xlab = "Time (days)", ylab = "Survival rate", col=c('red','blue'), lwd=2, xlim =c(0, 115))
legend('topright', text.font = , cex=1, c("22 day-old pupae","29 day-old pupae"), col=c('red','blue'), lty = 1, lwd=2, box.lty = 1)


### survival: median Kaplan-Meier estimator 
survfit(Surv(day,status)~Pupal_age,data=data_1)
cox.surv <- coxph(Surv(day, status) ~ Pupal_age, data = data_1, method ="exact")
summary(cox.surv)
```



Fig 8:Survival at point in time-figures (Binomial method)
Figure 8a:survival at different time points (15,30,60 days) for each treatment (22 days)

```{r }




data_8a<- read.csv("Figure 8a.csv")
head(data_8a)
data_8a$Time_point<- as.factor(data_8a$Time_point)

Figure_8a<-ggplot(data_8a,aes(x=Time_point,y=surv_flies,fill=Time_point
))+geom_boxplot(position=position_dodge(0.8))+
  geom_jitter(position=position_dodge(0.8))+
  geom_boxplot(alpha=0.3) +
  labs(fill = "Time points(days)") + 
  facet_wrap(~Treatments,ncol = 4) +
  theme_bw(base_size = 16)


Figure_8a

tiff("Figure 8a.tiff", width = 7, height = 4, units = 'in', compression = 'lzw',res = 300)
plot(Figure_8a+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + theme(axis.text.x = element_blank()) + xlab(expression(bold("Time points (days)"))) + ylab(expression(bold("No.of surviving flies (Age 22 days)")))
dev.off()
```


Figure 8b:survival at different time points (15,30,60 days) for each treatment(29 days)
```{r }
data_8b<- read.csv("Figure 8b.csv")
head(data_8b)
data_8b$Time_point<- as.factor(data_8b$Time_point)

Figure_8b<-ggplot(data_8b,aes(x=Time_point,y=surv_flies,fill=Time_point
))+geom_boxplot(position=position_dodge(0.8))+
  geom_jitter(position=position_dodge(0.8))+
  geom_boxplot(alpha=0.3) +
  labs(fill = "Time points(days)") + 
  facet_wrap(~Treatments,ncol = 4) +
  theme_bw(base_size = 16)

Figure_8b

tiff("Figure 8b.tiff", width = 7, height = 4, units = 'in',compression = 'lzw', res = 300)
plot(Figure_8b+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + theme(axis.text.x = element_blank()) + xlab(expression(bold("Time points (days)"))) + ylab(expression(bold("No.of surviving flies (Age 29 days) ")))
dev.off()
```


Figure 8c:survival at different time points (15,30,60 days) for each treatment (Combined data(22 and 29 days))
```{r }
data_8c<- read.csv("Figure 8c.csv")
head(data_8c)
data_8c$Time_point<- as.factor(data_8c$Time_point)


Figure_8c<-ggplot(data_8c,aes(x=Time_point,y=surv_flies,fill=Time_point
))+geom_boxplot(position=position_dodge(0.8))+
  geom_jitter(position=position_dodge(0.8))+
  geom_boxplot(alpha=0.3) +
  labs(fill = "Time points(days)") + 
  facet_wrap(~Treatments,ncol = 4) +
  theme_bw(base_size = 16)

Figure_8c

tiff("Figure 8c.tiff", width = 7, height = 4, units = 'in', res = 300)
plot(Figure_8c+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + theme(axis.text.x = element_blank()) + xlab(expression(bold("Time points (days)"))) + ylab(expression(bold("No. of surviving flies (Age 22 and 29 days)")))
dev.off()
```



Figure 8d: survival at different time points- differences between Pupal_age -Regardless of treatment
```{r }
data_8d<- read.csv("Figure 8d.csv")
head(data_8d)
data_8d$Pupal_age<- as.factor(data_8d$Pupal_age)
data_8d$Time_point<- as.factor(data_8d$Time_point)

Figure_8d<-ggplot(data_8d,aes(x=Time_point,y=surv_flies,fill=Time_point))+
  geom_boxplot(alpha=0.3) + 
  labs(fill = "Time points(days)") + 
  facet_wrap(~Pupal_age) +theme_bw(base_size = 16)

Figure_8d

tiff("Figure 8d.tiff", width = 7, height = 4, units = 'in', compression = 'lzw',res = 300)
plot(Figure_8d+theme_tufte() + theme(axis.line = element_line(size = 1, colour = "black")) + theme(legend.title = element_text(face = "bold")) + theme(legend.text = element_text(face = "italic"))) + theme(axis.text.x = element_blank()) + xlab(expression(bold("Time points (days)"))) + ylab(expression(bold("No. of surviving flies (Age 22 and 29 days) ")))
dev.off()



```



Figure  8: Survival at point in time(####Significant Values)-Binomial method


Figure 8a:survival  point in time -22 days
```{r }

data_4<- read.csv("Figure 8a_22.csv")
head(data_4)
boxplot(data_4$surv15 ~ data_4$Treatments, ylab = "Survival rate after 15 days_22")
boxplot(data_4$surv30 ~ data_4$Treatments, ylab = "Survival rate after 30days_22")
boxplot(data_4$surv60 ~ data_4$Treatments, ylab = "Survival rate after 60 days_22")

## survival: binomial model survival rate 15 days(difference in treatments)
fm4_1 <- glmer(cbind(surv15,n - surv15) ~ Treatment +(1|rep), family = binomial, data = data_4)
summary(fm4_1)

## survival: binomial model survival rate 30 days(difference in treatments)
fm4_2 <- glmer(cbind(surv30,n - surv30) ~ Treatment +(1|rep), family = binomial, data = data_4)
summary(fm4_2)

## survival: binomial model survival rate 60 days
fm4_3<- glmer(cbind(surv60,n - surv60) ~ Treatment +(1|rep), family = binomial, data = data_4)
summary(fm4_3)
```



Figure 8a:survival  point in time -29 days


```{r }
data_5<- read.csv("Figure 8a_29.csv")
head(data_5)
boxplot(data_5$surv15/data_5$n ~ data_5$Treatments, ylab = "Survival rate after 15 days_29")
boxplot(data_5$surv30/data_5$n ~ data_5$Treatments, ylab = "Survival rate after 30 days_29")
boxplot(data_5$surv60/data_5$n ~ data_5$Treatments, ylab = "Survival rate after 60 days_29")

### survival: binomial model survival rate 15 days (difference in treatments)
fm5_1 <- glmer(cbind(surv15,n - surv15) ~ Treatment +(1|rep), family = binomial, data = data_5)
summary(fm5_1)

### survival: binomial model survival rate 30 days(difference in treatments)
fm5_2 <- glmer(cbind(surv30,n - surv30) ~ Treatment +(1|rep), family = binomial, data = data_5)
summary(fm5_2)

### survival: binomial model survival rate 60 days(difference in treatments)
fm5_3 <- glmer(cbind(surv60,n - surv60) ~ Treatment +(1|rep), family = binomial, data = data_5)
summary(fm5_3)
```


Survival at point in time-Combined data(22 and 29 days) -Figure 8c and 8d ###significance

```{r }
data_4<- read.csv("Figure 8c_8d_significance.csv")
str(data_4)
summary(data_4)
head(data_4)

#### Individual boxplots for Figure 8c and 8d

Survival_15days<-data_4$surv15

boxplot(Survival_15days ~ data_4$Treatment, ylab = "Survival rate after 15 days (combined data)")

boxplot(Survival_15days ~ data_4$Pupal_age, ylab = "Survival rate after 15 days")

#30
boxplot(data_4$surv30/data_4$n ~ data_4$Treatment, ylab = "Survival rate after 30 days (combined data)")

boxplot(data_4$surv30/data_4$n ~ data_4$Pupal_age, ylab = "Survival rate after 30 days")

#60
boxplot(data_4$surv60/data_4$n ~ data_4$Treatment, ylab = "Survival rate after 60 days(combined data)")

boxplot(data_4$surv60/data_4$n ~ data_4$Pupal_age, ylab = "Survival rate after 60 days")





#### significance (boxplots above) for Figure 8c and 8d

########################## survival: binomial model survival rate 15 days


###Treatments 
fm4_1 <- glmer(cbind(surv15,n - surv15) ~ Treatment +(1|rep), family = binomial, data = data_4)
summary(fm4_1)

###Pupal age

fm4_1 <- glmer(cbind(surv15,n - surv15) ~ Pupal_age +(1|rep), family = binomial, data = data_4)
summary(fm4_1)

########################### survival: binomial model survival rate 30 days

###Treatments 
fm4_2 <- glmer(cbind(surv30,n - surv30) ~ Treatment +(1|rep), family = binomial, data = data_4)
summary(fm4_2)

###Pupal age

fm4_2 <- glmer(cbind(surv30,n - surv30) ~ Pupal_age +(1|rep), family = binomial, data = data_4)
summary(fm4_1)

########################## survival: binomial model survival rate 60 days

###Treatments 
fm4_3<- glmer(cbind(surv60,n - surv60) ~ Treatment +(1|rep), family = binomial, data = data_4)
summary(fm4_3)


###Pupal age

fm4_3<- glmer(cbind(surv60,n - surv60) ~ Pupal_age +(1|rep), family = binomial, data = data_4)
summary(fm4_3)

```











MODELS ON THE IMPACT OF SHOCK ON THE EMERGENCE RATE OF 22 AND 29 DAY OLD PUPAE (COMBINED DATA)


Figure 9a



```{r}
 

tab <- read.csv("Figure 9a.csv")
head(tab)
str(tab)

tab$pcemerg <- tab$emerged / (tab$emerged+tab$unemerged)

boxplot(tab$pcemerg ~ tab$Treatments, xlab ="Treatments", ylab = "Emergence rate")




plot(tab$pcemerg~tab$RH._Max)#*
plot(tab$pcemerg~tab$RH._mean)
plot(tab$pcemerg~tab$Temp_mean)#*
plot(tab$pcemerg~tab$Temp_Max)

#Scatter plots at threshold 5-Max
par(mfrow = c(2,2))
plot(tab$pcemerg~tab$Events._5max)
plot(tab$pcemerg~tab$Duration.ms._5max)
plot(tab$pcemerg~tab$scalar_max_5max)
plot(tab$pcemerg~tab$scalar_mean_5max)
plot(tab$pcemerg~tab$changescalar_max_5max)
plot(tab$pcemerg~tab$changescalar_mean_5max)#*
plot(tab$pcemerg~tab$Changevector_max_5max)    
plot(tab$pcemerg~tab$Changevector_mean_5max)
plot(tab$pcemerg~tab$Angle_max_5max)        
plot(tab$pcemerg~tab$Angle_mean_5max)
plot(tab$pcemerg~tab$Hz_max_5max)    
plot(tab$pcemerg~tab$Hz_mean_5max)
```

#Scatter plots at threshold 5-Mean

```{r}
par(mfrow = c(2,2))
plot(tab$pcemerg~tab$Duration.ms._5mean)#*
plot(tab$pcemerg~tab$scalar_max_5mean)#*
plot(tab$pcemerg~tab$scalar_mean_5mean)
plot(tab$pcemerg~tab$changescalar_max_5mean)
plot(tab$pcemerg~tab$changescalar_mean_5mean)
plot(tab$pcemerg~tab$Changevector_max_5mean)    
plot(tab$pcemerg~tab$Changevector_mean_5mean)
plot(tab$pcemerg~tab$Angle_max_5mean)        
plot(tab$pcemerg~tab$Angle_mean_5mean)
plot(tab$pcemerg~tab$Hz_max_5mean)    
plot(tab$pcemerg~tab$Hz_mean_5mean)
```


#Scatter plots at threshold 10-Max

```{r}
par(mfrow = c(2,2))
plot(tab$pcemerg~tab$Events._mean10)
plot(tab$pcemerg~tab$Duration.ms._10max)#*
plot(tab$pcemerg~tab$scalar_max_10max)
plot(tab$pcemerg~tab$scalar_mean_10max)
plot(tab$pcemerg~tab$changescalar_max_10max)
plot(tab$pcemerg~tab$changescalar_mean_10max)#*
plot(tab$pcemerg~tab$Changevector_max_10max)    
plot(tab$pcemerg~tab$Changevector_mean_10max)
plot(tab$pcemerg~tab$Angle_max_10max)        
plot(tab$pcemerg~tab$Angle_mean_10max)#*
plot(tab$pcemerg~tab$Hz_max_10max)    
plot(tab$pcemerg~tab$Hz_mean_10max)#*
```


#Scatter plots at threshold 10-Mean
```{r}
par(mfrow = c(2,2))
plot(tab$pcemerg~tab$Duration.ms._10mean)
plot(tab$pcemerg~tab$scalar_max_10mean)
plot(tab$pcemerg~tab$scalar_mean_10mean)#*
plot(tab$pcemerg~tab$changescalar_max_10mean)
plot(tab$pcemerg~tab$changescalar_mean_10mean)
plot(tab$pcemerg~tab$Changevector_max_10mean)    
plot(tab$pcemerg~tab$Changevector_mean_10mean)
plot(tab$pcemerg~tab$Angle_max_10mean)        
plot(tab$pcemerg~tab$Angle_mean_10mean)
plot(tab$pcemerg~tab$Hz_max_10mean)    
plot(tab$pcemerg~tab$Hz_mean_10mean)
```


#Scatter plots at threshold 15-Max
```{r}
par(mfrow = c(2,2))
plot(tab$pcemerg~tab$Events_15mean)
plot(tab$pcemerg~tab$Duration.ms._15max)
plot(tab$pcemerg~tab$scalar_max_15max)
plot(tab$pcemerg~tab$scalar_mean_15max)
plot(tab$pcemerg~tab$changescalar_max_15max)
plot(tab$pcemerg~tab$changescalar_mean_15max)#*
cor.test(tab$pcemerg,tab$changescalar_mean_15max)
plot(tab$pcemerg~tab$Changevector_max_15max)    
plot(tab$pcemerg~tab$Changevector_mean_15max)#
cor.test(tab$pcemerg,tab$Changevector_mean_15max)
plot(tab$pcemerg~tab$Angle_max_15max)        
plot(tab$pcemerg~tab$Angle_mean_15max)#*
plot(tab$pcemerg~tab$Hz_max_15max)    
plot(tab$pcemerg~tab$Hz_mean_15max)#*
```



#Scatter plots at threshold 15-Mean

```{r}
par(mfrow = c(2,2))
plot(tab$pcemerg~tab$Duration.ms._15mean)
plot(tab$pcemerg~tab$scalar_max_15mean)
plot(tab$pcemerg~tab$scalar_mean_15mean)#*
plot(tab$pcemerg~tab$changescalar_max_15mean)
plot(tab$pcemerg~tab$changescalar_mean_15mean)#*
plot(tab$pcemerg~tab$Changevector_max_15mean) #*   
plot(tab$pcemerg~tab$Changevector_mean_15mean)
plot(tab$pcemerg~tab$Angle_max_15mean)        
plot(tab$pcemerg~tab$Angle_mean_15mean)
plot(tab$pcemerg~tab$Hz_max_15mean)    
plot(tab$pcemerg~tab$Hz_mean_15mean)
```


#Scatter plots at threshold 20-Max

```{r}
par(mfrow = c(2,2))
plot(tab$pcemerg~tab$Events_20max)
plot(tab$pcemerg~tab$Duration.ms._20max)
plot(tab$pcemerg~tab$scalar_max_20max)
plot(tab$pcemerg~tab$scalar_mean_20max)#*
plot(tab$pcemerg~tab$changescalar_max_20max)
plot(tab$pcemerg~tab$changescalar_mean_20max)
plot(tab$pcemerg~tab$Changevector_max_20max)    
plot(tab$pcemerg~tab$Changevector_mean_20max)
plot(tab$pcemerg~tab$Angle_max_20max)        
plot(tab$pcemerg~tab$Angle_mean_20max)
plot(tab$pcemerg~tab$Hz_max_20max)    
plot(tab$pcemerg~tab$Hz_mean_20max)
```


#Scatter plots at threshold 20-Mean
```{r}
par(mfrow = c(2,2))
plot(tab$pcemerg~tab$Duration.ms._20mean)
plot(tab$pcemerg~tab$scalar_max_20mean)
plot(tab$pcemerg~tab$scalar_mean_20mean)
plot(tab$pcemerg~tab$changescalar_max_20mean)
plot(tab$pcemerg~tab$changescalar_mean_20mean)#*
plot(tab$pcemerg~tab$Changevector_max_20mean) 
plot(tab$pcemerg~tab$Changevector_mean_20mean)
plot(tab$pcemerg~tab$Angle_max_20mean)        
plot(tab$pcemerg~tab$Angle_mean_20mean)
plot(tab$pcemerg~tab$Hz_max_20mean)    
plot(tab$pcemerg~tab$Hz_mean_20mean)
```

##Models of impact of shock on emergence rate of transported pupae 

```{r}
#models threshold 5
fm1 <- glmer(cbind(emerged,unemerged) ~ Irradiation + RH._Max + Temp_mean + (1|Replicate), family = binomial, data = tab)
fm2 <- glmer(cbind(emerged,unemerged) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + scalar_max_5mean + (1|Replicate), family = binomial, data = tab)
fm3 <- glmer(cbind(emerged,unemerged) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Duration.ms._5mean +(1|Replicate), family = binomial, data = tab)
fm4 <- glmer(cbind(emerged,unemerged) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + (1|Replicate), family = binomial, data = tab)
fm5 <- glmer(cbind(emerged,unemerged) ~ Irradiation + RH._Max + Temp_mean + scalar_max_5mean + (1|Replicate), family = binomial, data = tab)
fm6 <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean + changescalar_mean_5max + scalar_max_5mean + (1|Replicate), family = binomial, data = tab)
fm7 <- glmer(cbind(emerged,unemerged) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + scalar_max_5mean + (1|Replicate), family = binomial, data = tab)
fm9 <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean + changescalar_mean_5max +  (1|Replicate), family = binomial, data = tab)
fm10 <- glmer(cbind(emerged,unemerged) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Duration.ms._5mean +(1|Replicate), family = binomial, data = tab)
fm11 <- glmer(cbind(emerged,unemerged) ~ Irradiation + RH._Max + Temp_mean + (1|Replicate), family = binomial, data = tab)
fm12 <- glmer(cbind(emerged,unemerged) ~ Irradiation +Temp_mean + Age +(1|Replicate), family = binomial, data = tab)
fm13 <- glmer(cbind(emerged,unemerged) ~ Irradiation +RH._Max + Age +(1|Replicate), family = binomial, data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm9,fm10,fm11,fm12,fm13)
summary(fm6)

fm11a <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean + Duration.ms._10max + changescalar_mean_10max+Angle_mean_10max+ (1|Replicate), family = binomial, data = tab)
fm12b <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean +changescalar_mean_10max+ scalar_mean_10mean +(1|Replicate), family = binomial, data = tab)
fm13c <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean + Duration.ms._10max + changescalar_mean_10max+ (1|Replicate), family = binomial, data = tab)
fm14 <- glmer(cbind(emerged,unemerged) ~ Irradiation + Duration.ms._10max + scalar_mean_10mean +(1|Replicate), family = binomial, data = tab)
fm15 <- glmer(cbind(emerged,unemerged) ~ Irradiation + changescalar_mean_10max+ scalar_mean_10mean +(1|Replicate), family = binomial, data = tab)
fm16 <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean + Duration.ms._10max + Angle_mean_10max+ (1|Replicate), family = binomial, data = tab)
fm17 <- glmer(cbind(emerged,unemerged) ~ Irradiation +  scalar_mean_10mean +Angle_mean_10max+  (1|Replicate), family = binomial, data = tab)
fm18 <- glmer(cbind(emerged,unemerged) ~ Irradiation + Age + changescalar_mean_10max+ scalar_mean_10mean +(1|Replicate), family = binomial, data = tab)
fm19 <- glmer(cbind(emerged,unemerged) ~ Irradiation + Duration.ms._10max + Age + (1|Replicate), family = binomial, data = tab)
fm14a <- glmer(cbind(emerged,unemerged) ~ Irradiation + Duration.ms._10max +  (1|Replicate), family = binomial, data = tab)
fm14b <- glmer(cbind(emerged,unemerged) ~ Irradiation + Duration.ms._10max + changescalar_mean_10max+ Age +(1|Replicate), family = binomial, data = tab)
fm14c <- glmer(cbind(emerged,unemerged) ~ Irradiation +  changescalar_mean_10max+ (1|Replicate), family = binomial, data = tab)


AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm9,fm10,fm11,fm12,fm13,fm11a,fm12b,fm13c,fm14,fm14a,fm14b,fm14c,fm15,fm16,fm17,fm18,fm19)
summary(fm6)

fm21 <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean + changescalar_mean_15max + Angle_mean_15max+ changescalar_mean_15mean+Changevector_max_15mean+ (1|Replicate), family = binomial, data = tab)
fm22 <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean + changescalar_mean_15max + Changevector_mean_15max+ Angle_mean_15max+ (1|Replicate), family = binomial, data = tab)
fm23 <- glmer(cbind(emerged,unemerged) ~ Irradiation + changescalar_mean_15max + Changevector_mean_15max+ Angle_mean_15max+ (1|Replicate), family = binomial, data = tab)
fm24 <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean +changescalar_mean_15max + changescalar_mean_15mean+  (1|Replicate), family = binomial, data = tab)
fm25 <- glmer(cbind(emerged,unemerged) ~ Irradiation + changescalar_mean_15max + Age+ Angle_mean_15max+ (1|Replicate), family = binomial, data = tab)
fm26 <- glmer(cbind(emerged,unemerged) ~ Irradiation + changescalar_mean_15max + scalar_mean_15mean + (1|Replicate), family = binomial, data = tab)
fm27 <- glmer(cbind(emerged,unemerged) ~ Irradiation + Changevector_mean_15max+Angle_mean_15max+scalar_mean_15mean +Changevector_max_15mean+ (1|Replicate), family = binomial, data = tab)
fm28 <- glmer(cbind(emerged,unemerged) ~ Irradiation + Changevector_mean_15max+(1|Replicate), family = binomial, data = tab)
fm29 <- glmer(cbind(emerged,unemerged) ~ Irradiation + changescalar_mean_15max + Changevector_mean_15max+ (1|Replicate), family = binomial, data = tab)
fm30 <- glmer(cbind(emerged,unemerged) ~ Irradiation + scalar_mean_15mean + Age +(1|Replicate), family = binomial, data = tab)
fm21a <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean +  Changevector_mean_15max+ Angle_mean_15max+changescalar_mean_15mean+ (1|Replicate), family = binomial, data = tab)
fm21b <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean + changescalar_mean_15max + scalar_mean_15mean +Changevector_max_15mean+ (1|Replicate), family = binomial, data = tab)
fm21c <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean + changescalar_mean_15max + Angle_mean_15max+ (1|Replicate), family = binomial, data = tab)
fm21d <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean +changescalar_mean_15mean+Changevector_max_15mean+ (1|Replicate), family = binomial, data = tab)
fm21e <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean +changescalar_mean_15mean+Changevector_max_15mean+ Age +(1|Replicate), family = binomial, data = tab)
fm27a <- glmer(cbind(emerged,unemerged) ~ Irradiation + scalar_mean_15mean +Changevector_max_15mean+ (1|Replicate), family = binomial, data = tab)
fm27b <- glmer(cbind(emerged,unemerged) ~ Irradiation + Changevector_mean_15max+Angle_mean_15max+scalar_mean_15mean + (1|Replicate), family = binomial, data = tab)
fm29a <- glmer(cbind(emerged,unemerged) ~ Irradiation + Changevector_mean_15max+ (1|Replicate), family = binomial, data = tab)



AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm9,fm10,fm11,fm12,fm13,fm11a,fm12b,fm13c,fm14,fm14a,fm14b,fm14c,fm15,fm16,fm17,fm18,fm19,fm21,fm22,fm23,fm24,fm25,fm26,fm27,fm28,fm29,fm29a,fm30,fm21a,fm21b,fm21c,fm21d,fm21e,fm27a,fm27b)
summary(fm29)

fm31 <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean + scalar_mean_20max + changescalar_mean_20mean+ (1|Replicate), family = binomial, data = tab)
fm32 <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean +  changescalar_mean_20mean+ (1|Replicate), family = binomial, data = tab)
fm31a <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean + scalar_mean_20max + (1|Replicate), family = binomial, data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm9,fm10,fm11,fm12,fm13,fm11a,fm12b,fm13c,fm14,fm14a,fm14b,fm14c,fm15,fm16,fm17,fm18,fm19,fm21,fm22,fm23,fm24,fm25,fm26,fm27,fm28,fm29,fm29a,fm30,fm21a,fm21b,fm21c,fm21d,fm21e,fm27a,fm27b,fm31,fm31a,fm32)
summary(fm29)

#combination of thresholds
fm41 <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean + Duration.ms._10max + scalar_max_5mean + Changevector_mean_15max+(1|Replicate), family = binomial, data = tab)
fm42 <- glmer(cbind(emerged,unemerged) ~ Irradiation + scalar_max_5mean + Changevector_mean_15max+ (1|Replicate), family = binomial, data = tab)
fm41a <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean + Duration.ms._10max + changescalar_mean_5max + (1|Replicate), family = binomial, data = tab)
fm41b <- glmer(cbind(emerged,unemerged) ~ Irradiation + Temp_mean + Duration.ms._10max + scalar_max_5mean + (1|Replicate), family = binomial, data = tab)


AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm9,fm10,fm11,fm12,fm13,fm11a,fm12b,fm13c,fm14,fm14a,fm14b,fm14c,fm15,fm16,fm17,fm18,fm19,fm21,fm22,fm23,fm24,fm25,fm26,fm27,fm28,fm29,fm29a,fm30,fm21a,fm21b,fm21c,fm21d,fm21e,fm27a,fm27b,fm31,fm31a,fm32,fm41,fm42,fm41a,fm41b)
summary(fm29)

#the best model remains as fm29 (s12 fig)

##Fig9a
head(tab)
str(tab)
summary(fm29)
plot((emerged / unemerged) ~ fitted(fm29), data = tab)
abline(lm((emerged / unemerged) ~ fitted(fm29), data = tab), col = "red")
cor.test(tab$emerged/tab$unemerged,fitted(fm29))
summary(lm(fitted(fm29)~(tab$emerged/tab$unemerged)))$r.squared


#####Correlation of the factors in the best model(fm29) 
cor.test(tab$changescalar_mean_15max,tab$Changevector_mean_15max)


####Conclusion: significantly correlated. Settled on Changevector_mean_15max
```







MODELS ON THE IMPACT OF SHOCK ON THE flight propensity OF 22 AND 29 DAY OLD PUPAE (COMBINED DATA)

Figure 9b:Flight propensity 




```{r}



tab <- read.csv("Figure 9b.csv")
head(tab)
str(tab)
tab$pcflight <- tab$out / (tab$out+tab$in.)
boxplot(tab$pcflight~ tab$Treatments, xlab ="Treatments", ylab = "Flight propensity")

plot(tab$pcflight~tab$RH._Max)#*
plot(tab$pcflight~tab$RH._mean)
plot(tab$pcflight~tab$Temp_mean)#*
cor.test(tab$pcflight,tab$Temp_mean)
plot(tab$pcflight~tab$Temp_Max)

par(mfrow = c(2,2))
plot(tab$pcflight~tab$Events._5max)


```
#Scatter plots at threshold 5-Max
```{r }
par(mfrow = c(2,2))
plot(tab$pcflight~tab$Duration.ms._5max)#*
plot(tab$pcflight~tab$scalar_max_5max)
plot(tab$pcflight~tab$scalar_mean_5max)#*
plot(tab$pcflight~tab$changescalar_max_5max)
plot(tab$pcflight~tab$changescalar_mean_5max)*#
plot(tab$pcflight~tab$Changevector_max_5max)    
plot(tab$pcflight~tab$Changevector_mean_5max)
plot(tab$pcflight~tab$Angle_max_5max)#*        
plot(tab$pcflight~tab$Angle_mean_5max)
plot(tab$pcflight~tab$Hz_max_5max)#*
plot(tab$pcflight~tab$Hz_mean_5max)


```
#Scatter plots at threshold 5-Mean

```{r }
par(mfrow = c(2,2))
plot(tab$pcflight~tab$Duration.ms._5mean)#*
plot(tab$pcflight~tab$scalar_max_5mean)
plot(tab$pcflight~tab$scalar_mean_5mean)
plot(tab$pcflight~tab$changescalar_max_5mean)
plot(tab$pcflight~tab$changescalar_mean_5mean)
plot(tab$pcflight~tab$Changevector_max_5mean)    
plot(tab$pcflight~tab$Changevector_mean_5mean)
plot(tab$pcflight~tab$Angle_max_5mean)        
plot(tab$pcflight~tab$Angle_mean_5mean)
plot(tab$pcflight~tab$Hz_max_5mean)    
plot(tab$pcflight~tab$Hz_mean_5mean)


```
#Scatter plots at threshold 10-Max

```{r }
par(mfrow = c(2,2))
plot(tab$pcflight~tab$Events_10max)
plot(tab$pcflight~tab$Duration.ms._10max)
plot(tab$pcflight~tab$scalar_max_10max)
plot(tab$pcflight~tab$scalar_mean_10max)#*
plot(tab$pcflight~tab$changescalar_max_10max)
plot(tab$pcflight~tab$changescalar_mean_10max)
plot(tab$pcflight~tab$Changevector_max_10max)    
plot(tab$pcflight~tab$Changevector_mean_10max)#*
plot(tab$pcflight~tab$Angle_max_10max)        
plot(tab$pcflight~tab$Angle_mean_10max)
plot(tab$pcflight~tab$Hz_max_10max)    
plot(tab$pcflight~tab$Hz_mean_10max)#*

```
#Scatter plots at threshold 10-Mean

```{r }

par(mfrow = c(2,2))
plot(tab$pcflight~tab$Duration.ms._10mean)#*
plot(tab$pcflight~tab$scalar_max_10mean)#*
plot(tab$pcflight~tab$scalar_mean_10mean)
plot(tab$pcflight~tab$changescalar_max_10mean)
plot(tab$pcflight~tab$changescalar_mean_10mean)#*
plot(tab$pcflight~tab$Changevector_max_10mean)    
plot(tab$pcflight~tab$Changevector_mean_10mean)#*
plot(tab$pcflight~tab$Angle_max_10mean)        
plot(tab$pcflight~tab$Angle_mean_10mean)
plot(tab$pcflight~tab$Hz_max_10mean)    
plot(tab$pcflight~tab$Hz_mean_10mean)
```


#Scatter plots at threshold 15-Max

```{r }

par(mfrow = c(2,2))
plot(tab$pcflight~tab$Events_15max)
plot(tab$pcflight~tab$Duration.ms._15max)
plot(tab$pcflight~tab$scalar_max_15max)
plot(tab$pcflight~tab$scalar_mean_15max)
plot(tab$pcflight~tab$changescalar_max_15max)
plot(tab$pcflight~tab$changescalar_mean_15max)#*
plot(tab$pcflight~tab$Changevector_max_15max)    
plot(tab$pcflight~tab$Changevector_mean_15max)#*
plot(tab$pcflight~tab$Angle_max_15max)        
plot(tab$pcflight~tab$Angle_mean_15max)#*
cor.test(tab$pcflight,tab$Angle_mean_15max)
plot(tab$pcflight~tab$Hz_max_15max)    
plot(tab$pcflight~tab$Hz_mean_15max)#*

```

#Scatter plots at threshold 15-Max-mean

```{r }
par(mfrow = c(2,2))

plot(tab$pcflight~tab$Duration.ms._15mean)#*
plot(tab$pcflight~tab$scalar_max_15mean)
plot(tab$pcflight~tab$scalar_mean_15mean)#*
plot(tab$pcflight~tab$changescalar_max_15mean)
plot(tab$pcflight~tab$changescalar_mean_15mean)#*
plot(tab$pcflight~tab$Changevector_max_15mean)   
plot(tab$pcflight~tab$Changevector_mean_15mean)#*
plot(tab$pcflight~tab$Angle_max_15mean)        
plot(tab$pcflight~tab$Angle_mean_15mean)
plot(tab$pcflight~tab$Hz_max_15mean)    
plot(tab$pcflight~tab$Hz_mean_15mean)

```

#Scatter plots at threshold 20-Max


```{r }
par(mfrow = c(2,2))

plot(tab$pcflight~tab$Events_20max)
plot(tab$pcflight~tab$Duration.ms._20max)#*
plot(tab$pcflight~tab$scalar_max_20max)
plot(tab$pcflight~tab$scalar_mean_20max)
plot(tab$pcflight~tab$changescalar_max_20max)#*
plot(tab$pcflight~tab$changescalar_mean_20max)
plot(tab$pcflight~tab$Changevector_max_20max)    
plot(tab$pcflight~tab$Changevector_mean_20max)
plot(tab$pcflight~tab$Angle_max_20max)        
plot(tab$pcflight~tab$Angle_mean_20max)#*
plot(tab$pcflight~tab$Hz_max_20max)    
plot(tab$pcflight~tab$Hz_mean_20max)
```


#Scatter plots at threshold 20-Mean
```{r }
par(mfrow = c(2,2))

plot(tab$pcflight~tab$Duration.ms._20mean)
plot(tab$pcflight~tab$scalar_max_20mean)
plot(tab$pcflight~tab$scalar_mean_20mean)
plot(tab$pcflight~tab$changescalar_max_20mean)#*
plot(tab$pcflight~tab$changescalar_mean_20mean)
plot(tab$pcflight~tab$Changevector_max_20mean) #*
plot(tab$pcflight~tab$Changevector_mean_20mean)
plot(tab$pcflight~tab$Angle_max_20mean)        
plot(tab$pcflight~tab$Angle_mean_20mean)#*
plot(tab$pcflight~tab$Hz_max_20mean)    
plot(tab$pcflight~tab$Hz_mean_20mean)
```

#Models of impact of transportation on the flight propensity 

```{r }

factor(tab$Age)

#####Caroline 
#models threshold 5
fm1 <- glmer(cbind(out,in.) ~ Irradiation + RH._Max + Temp_mean + (1|Replicate), family = binomial, data = tab)
fm2 <- glmer(cbind(out,in.) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Changevector_mean_5max + (1|Replicate), family = binomial, data = tab)
fm3 <- glmer(cbind(out,in.) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Duration.ms._5mean +(1|Replicate), family = binomial, data = tab)
fm4 <- glmer(cbind(out,in.) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max +scalar_mean_5max + (1|Replicate), family = binomial, data = tab)
fm5 <- glmer(cbind(out,in.) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Duration.ms._5mean + Changevector_mean_5max + (1|Replicate), family = binomial, data = tab)
fm6 <- glmer(cbind(out,in.) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Angle_mean_5max + Duration.ms._5mean + (1|Replicate), family = binomial, data = tab)
fm7 <- glmer(cbind(out,in.) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Changevector_mean_5max +Duration.ms._5mean +(1|Replicate), family = binomial, data = tab)
fm8 <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean + changescalar_mean_5max + Duration.ms._5max +Angle_mean_5max + (1|Replicate), family = binomial, data = tab)
fm1a<- glmer(cbind(out,in.) ~ Irradiation + RH._Max + Temp_mean + Duration.ms._5max +scalar_mean_5max + (1|Replicate), family = binomial, data = tab)
fm1b <- glmer(cbind(out,in.) ~ Irradiation + RH._Max + Temp_mean + scalar_mean_5max +changescalar_mean_5mean +Changevector_mean_5max+(1|Replicate), family = binomial, data = tab)
fm1c <- glmer(cbind(out,in.) ~ Irradiation + RH._Max + Temp_mean + scalar_mean_5max +changescalar_mean_5max +Changevector_mean_5mean + (1|Replicate), family = binomial, data = tab)
fm1d <- glmer(cbind(out,in.) ~ Irradiation + RH._Max + Temp_mean + scalar_mean_5max +changescalar_mean_5mean +Changevector_mean_5mean +(1|Replicate), family = binomial, data = tab)
fm1e <- glmer(cbind(out,in.) ~ Irradiation + RH._Max + Temp_mean + scalar_mean_5max +(1|Replicate), family = binomial, data = tab)
fm1f <- glmer(cbind(out,in.) ~ Irradiation + RH._Max + Temp_mean + Angle_mean_5max +(1|Replicate), family = binomial, data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm1a,fm1b,fm1c,fm1d,fm1e,fm1f)
summary(fm1)


##models for threshold 10

fm11 <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean + Duration.ms._10max + changescalar_mean_10max+Angle_mean_10max+ (1|Replicate), family = binomial, data = tab)
fm12 <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean + Duration.ms._10max + changescalar_mean_10max+ Changevector_mean_10mean + (1|Replicate), family = binomial, data = tab)
fm13 <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean + Duration.ms._10max + changescalar_mean_10max+ Changevector_mean_10max + (1|Replicate), family = binomial, data = tab)
fm14 <- glmer(cbind(out,in.) ~ Irradiation + Duration.ms._10max + changescalar_mean_10max+Changevector_mean_10max + (1|Replicate), family = binomial, data = tab)
fm15 <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean +Duration.ms._10mean + changescalar_mean_10max+ (1|Replicate), family = binomial, data = tab)
fm11a <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean + Duration.ms._10max + Changevector_mean_10mean +(1|Replicate), family = binomial, data = tab)
fm11b <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean + Duration.ms._10mean + changescalar_mean_10max+scalar_max_10mean +(1|Replicate), family = binomial, data = tab)
fm11c <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean + Duration.ms._10mean + Changevector_mean_10mean +scalar_max_10mean +(1|Replicate), family = binomial, data = tab)
fm11d <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean +Duration.ms._10mean +Changevector_mean_10mean + (1|Replicate), family = binomial, data = tab)
fm12a <- glmer(cbind(in.,out) ~ Temp_mean + Angle_mean_10max+ (1|Replicate), family = binomial, data = tab)
fm12b <- glmer(cbind(in.,out) ~ Temp_mean + Duration.ms._10max+ (1|Replicate), family = binomial, data = tab)
fm12c <- glmer(cbind(in.,out) ~ Temp_mean + Angle_mean_10max+ Duration.ms._10max+(1|Replicate), family = binomial, data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm1a,fm1b,fm1c,fm1d,fm1e,fm1f,fm12,fm13,fm14,fm15,fm11a,fm11b,fm11c,fm11d,fm12a,fm12b,fm12c)
summary(fm12b)
#

##threshold 15
fm21 <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean + changescalar_mean_15max + scalar_mean_15mean + Duration.ms._15mean +(1|Replicate), family = binomial, data = tab)
fm22 <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean + changescalar_mean_15max + Changevector_mean_15mean +(1|Replicate), family = binomial, data = tab)
fm23 <- glmer(cbind(out,in.) ~ Irradiation + changescalar_mean_15max + Changevector_mean_15max+scalar_mean_15mean + (1|Replicate), family = binomial, data = tab)
fm24 <- glmer(cbind(out,in.) ~ Irradiation + changescalar_mean_15max + Changevector_mean_15max+ Angle_mean_15max+Duration.ms._15mean +(1|Replicate), family = binomial, data = tab)
fm25 <- glmer(cbind(out,in.) ~ Irradiation + changescalar_mean_15max + Changevector_mean_15max+ Angle_mean_15max+ (1|Replicate), family = binomial, data = tab)
fm26 <- glmer(cbind(out,in.) ~ Irradiation + changescalar_mean_15max + Changevector_mean_15max+scalar_mean_15mean + Duration.ms._15mean +(1|Replicate), family = binomial, data = tab)
fm27 <- glmer(cbind(out,in.) ~ Irradiation + Changevector_mean_15max+scalar_mean_15mean + (1|Replicate), family = binomial, data = tab)
fm28 <- glmer(cbind(out,in.) ~ Irradiation + Changevector_mean_15max+ changescalar_mean_15mean +(1|Replicate), family = binomial, data = tab)
fm29 <- glmer(cbind(out,in.) ~ Irradiation + changescalar_mean_15max + Changevector_mean_15max+(1|Replicate), family = binomial, data = tab)
fm30 <- glmer(cbind(out,in.) ~ Irradiation +Temp_mean +Angle_mean_15max + Changevector_mean_15mean +(1|Replicate), family = binomial, data = tab)
fm30a <-glmer(cbind(out,in.) ~ Irradiation +Temp_mean + Duration.ms._15mean +changescalar_mean_15mean + Changevector_mean_15mean +(1|Replicate), family = binomial, data = tab)
fm29a<- glmer(cbind(out,in.) ~ Irradiation + changescalar_mean_15max + Changevector_mean_15max+ Duration.ms._15mean +scalar_mean_15mean +(1|Replicate), family = binomial, data = tab)
fm29b <- glmer(cbind(in.,out) ~ Age +Temp_mean+ Angle_mean_15max+  (1|Replicate), family = binomial, data = tab)
fm29c <- glmer(cbind(in.,out) ~ Temp_mean + Angle_mean_15max+ (1|Replicate), family = binomial, data = tab)
fm29d <- glmer(cbind(in.,out) ~ Age + Angle_mean_15max+ (1|Replicate), family = binomial, data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm1a,fm1b,fm1c,fm1d,fm1e,fm1f,fm13,fm14,fm15,fm11a,fm11b,fm11c,fm11d,fm12a,fm12b,fm12c,fm21,fm22,fm23,fm24,fm25,fm26,fm27,fm28,fm29,fm30,fm30a,fm29a,fm29b,fm29c,fm29d)
summary(fm29b)

##Threshold 20

fm31 <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean + Duration.ms._20max + changescalar_mean_20max + Angle_mean_20max + (1|Replicate), family = binomial, data = tab)
fm32 <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean +  changescalar_mean_20mean+ Angle_mean_20mean +(1|Replicate), family = binomial, data = tab)
fm31a <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean + Duration.ms._20max + changescalar_mean_20max + Angle_mean_20mean + (1|Replicate), family = binomial, data = tab)
fm32a <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean +  changescalar_mean_20max +Angle_mean_20max +(1|Replicate), family = binomial, data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm1a,fm1b,fm1c,fm1d,fm1e,fm1f,fm13,fm14,fm15,fm11a,fm11b,fm11c,fm11d,fm12a,fm12b,fm12c,fm21,fm22,fm23,fm24,fm25,fm26,fm27,fm28,fm29,fm30,fm30a,fm29a,fm29b,fm29c,fm29d,fm31,fm31a,fm32,fm32a)
summary(fm29b)


#combination of thresholds
fm41 <- glmer(cbind(out,in.) ~ Irradiation + Temp_mean + RH._Max + changescalar_mean_10max + Angle_mean_15max + Changevector_mean_15mean +(1|Replicate), family = binomial, data = tab)
fm42 <- glmer(cbind(out,in.) ~ Irradiation +Temp_mean + RH._Max + changescalar_mean_20max +Angle_mean_20max + Duration.ms._10mean +(1|Replicate), family = binomial, data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm13,fm14,fm15,fm11a,fm11b,fm11c,fm11d,fm21,fm22,fm23,fm24,fm25,fm26,fm27,fm28,fm29,fm30,fm30a,fm29a,fm29b,fm29c,fm29d,fm31,fm32,fm31a,fm32a,fm41,fm42)

#the best model remains as fm29b 

##Figure 9b

head(tab)
str(tab)
summary(fm29b)
plot((out / in.) ~ fitted(fm29b), data = tab)
FIT <- lm(tab$pcflight ~ fitted(fm29b))
abline(reg = FIT)
cor.test(tab$pcflight,fitted(fm29b))
summary(lm(fitted(fm29b)~(tab$out/tab$in.)))$r.squared
#any(is.na(tab))

## Independent variables not correlated

```




MODELS ON THE IMPACT OF SHOCK ON THE Mating abilityOF 22 AND 29 DAY OLD PUPAE (COMBINED DATA)

Figure 9c: Mating ability


```{r}


tab <- read.csv("Figure 9c.csv")
head(tab)
summary(tab)
str(tab)
tab$pcmating <- tab$pairs_formed / (tab$pairs_formed+tab$unformed_pairs)
boxplot(tab$pcmating ~ tab$Treatments, xlab ="Treatments", ylab = "mating ability")


#####SCATTER PLOTS
plot(tab$pcmating~tab$RH._Max)#*
cor.test(tab$pcmating,tab$RH._Max)
plot(tab$pcmating~tab$RH._mean)
plot(tab$pcmating~tab$Temp_mean)#*
cor.test(tab$pcmating,tab$Temp_mean)
plot(tab$pcmating~tab$Temp_Max)

#Scatter plots at threshold 5-Max
par(mfrow = c(2,2))
plot(tab$pcmating~tab$Events._5max)
plot(tab$pcmating~tab$Duration.ms._5max)
plot(tab$pcmating~tab$scalar_max_5max)
plot(tab$pcmating~tab$scalar_mean_5max)
plot(tab$pcmating~tab$changescalar_max_5max)
plot(tab$pcmating~tab$changescalar_mean_5max)#*
cor.test(tab$pcmating,tab$changescalar_mean_5max)
plot(tab$pcmating~tab$Changevector_max_5max)    
plot(tab$pcmating~tab$Changevector_mean_5max)#*
cor.test(tab$pcmating,tab$Changevector_mean_5max)
plot(tab$pcmating~tab$Angle_max_5max)        
plot(tab$pcmating~tab$Angle_mean_5max)
plot(tab$pcmating~tab$Hz_max_5max)    
plot(tab$pcmating~tab$Hz_mean_5max)
```

#Scatter plots at threshold 5-Mean
```{r}
par(mfrow = c(2,2))
plot(tab$pcmating~tab$Duration.ms._5mean)#*
plot(tab$pcmating~tab$scalar_max_5mean)
plot(tab$pcmating~tab$scalar_mean_5mean)
plot(tab$pcmating~tab$changescalar_max_5mean)
plot(tab$pcmating~tab$changescalar_mean_5mean)
plot(tab$pcmating~tab$Changevector_max_5mean)    
plot(tab$pcmating~tab$Changevector_mean_5mean)
plot(tab$pcmating~tab$Angle_max_5mean)        
plot(tab$pcmating~tab$Angle_mean_5mean)
plot(tab$pcmating~tab$Hz_max_5mean)    
plot(tab$pcmating~tab$Hz_mean_5mean)
```

#Scatter plots at threshold 10-Max
```{r}
par(mfrow = c(2,2))
plot(tab$pcmating~tab$Events_10max)#*
plot(tab$pcmating~tab$Duration.ms._10max)
plot(tab$pcmating~tab$scalar_max_10max)
plot(tab$pcmating~tab$scalar_mean_10max)
plot(tab$pcmating~tab$changescalar_max_10max)#*
plot(tab$pcmating~tab$changescalar_mean_10max)
plot(tab$pcmating~tab$Changevector_max_10max)    
plot(tab$pcmating~tab$Changevector_mean_10max)
plot(tab$pcmating~tab$Angle_max_10max)#*        
plot(tab$pcmating~tab$Angle_mean_10max)
plot(tab$pcmating~tab$Hz_max_10max)    
plot(tab$pcmating~tab$Hz_mean_10max)
```

#Scatter plots at threshold 10-Mean

```{r}
par(mfrow = c(2,2))
plot(tab$pcmating~tab$Duration.ms._10mean)
plot(tab$pcmating~tab$scalar_max_10mean)
plot(tab$pcmating~tab$scalar_mean_10mean)
plot(tab$pcmating~tab$changescalar_max_10mean)
plot(tab$pcmating~tab$changescalar_mean_10mean)
plot(tab$pcmating~tab$Changevector_max_10mean)    
plot(tab$pcmating~tab$Changevector_mean_10mean)
plot(tab$pcmating~tab$Angle_max_10mean)        
plot(tab$pcmating~tab$Angle_mean_10mean)
plot(tab$pcmating~tab$Hz_max_10mean)    
plot(tab$pcmating~tab$Hz_mean_10mean)
```

#Scatter plots at threshold 15-Max
```{r}
par(mfrow = c(2,2))
plot(tab$pcmating~tab$Events_15max)
plot(tab$pcmating~tab$Duration.ms._15max)#*
plot(tab$pcmating~tab$scalar_max_15max)
plot(tab$pcmating~tab$scalar_mean_15max)
plot(tab$pcmating~tab$changescalar_max_15max)#*
plot(tab$pcmating~tab$changescalar_mean_15max)
plot(tab$pcmating~tab$Changevector_max_15max)    
plot(tab$pcmating~tab$Changevector_mean_15max)
plot(tab$pcmating~tab$Angle_max_15max)#*        
plot(tab$pcmating~tab$Angle_mean_15max)
plot(tab$pcmating~tab$Hz_max_15max)#*   
plot(tab$pcmating~tab$Hz_mean_15max)
```

#Scatter plots at threshold 15-Mean
```{r}
par(mfrow = c(2,2))
plot(tab$pcmating~tab$Duration.ms._15mean)#*
plot(tab$pcmating~tab$scalar_max_15mean)
plot(tab$pcmating~tab$scalar_mean_15mean)
plot(tab$pcmating~tab$changescalar_max_15mean)#*
plot(tab$pcmating~tab$changescalar_mean_15mean)
plot(tab$pcmating~tab$Changevector_max_15mean)   
plot(tab$pcmating~tab$Changevector_mean_15mean)
plot(tab$pcmating~tab$Angle_max_15mean)      
plot(tab$pcmating~tab$Angle_mean_15mean)
plot(tab$pcmating~tab$Hz_max_15mean)   
plot(tab$pcmating~tab$Hz_mean_15mean)
```

#Scatter plots at threshold 20-Max
```{r}
par(mfrow = c(2,2))
plot(tab$pcmating~tab$Events_20max)
plot(tab$pcmating~tab$Duration.ms._20max)
plot(tab$pcmating~tab$scalar_max_20max)
plot(tab$pcmating~tab$scalar_mean_20max)
plot(tab$pcmating~tab$changescalar_max_20max)
plot(tab$pcmating~tab$changescalar_mean_20max)
plot(tab$pcmating~tab$Changevector_max_20max) #*   
plot(tab$pcmating~tab$Changevector_mean_20max)
plot(tab$pcmating~tab$Angle_max_20max)        
plot(tab$pcmating~tab$Angle_mean_20max)#*
plot(tab$pcmating~tab$Hz_max_20max)#*    
plot(tab$pcmating~tab$Hz_mean_20max)#*
```

#Scatter plots at threshold 20-Mean
```{r}
par(mfrow = c(2,2))
plot(tab$pcmating~tab$Duration.ms._20mean)
plot(tab$pcmating~tab$scalar_max_20mean)#*
plot(tab$pcmating~tab$scalar_mean_20mean)
plot(tab$pcmating~tab$changescalar_max_20mean)
plot(tab$pcmating~tab$changescalar_mean_20mean)
plot(tab$pcmating~tab$Changevector_max_20mean) 
plot(tab$pcmating~tab$Changevector_mean_20mean)
plot(tab$pcmating~tab$Angle_max_20mean)        
plot(tab$pcmating~tab$Angle_mean_20mean)
plot(tab$pcmating~tab$Hz_max_20mean)    
plot(tab$pcmating~tab$Hz_mean_20mean)
```



```{r}
tab$Age <- as.factor(tab$Age)



#*****************************************************************************************************************************************************

#####Caroline 
#models threshold 5
fm1 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + Temp_mean + (1|Replicate), family = binomial, data = tab)
fm2 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Changevector_mean_5max + (1|Replicate), family = binomial, data = tab)
fm3 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Duration.ms._5mean +(1|Replicate), family = binomial, data = tab)
fm4 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + (1|Replicate), family = binomial, data = tab)
fm5 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Duration.ms._5mean + Changevector_mean_5max + (1|Replicate), family = binomial, data = tab)
fm6 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + Temp_mean  + Duration.ms._5mean + (1|Replicate), family = binomial, data = tab)
fm7 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Changevector_mean_5max +Duration.ms._5mean +(1|Replicate), family = binomial, data = tab)
fm8 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Temp_mean + changescalar_mean_5max + Duration.ms._5max +(1|Replicate), family = binomial, data = tab)
fm1a<- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + Temp_mean + Duration.ms._5mean + (1|Replicate), family = binomial, data = tab)
fm1b <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + (1|Replicate), family = binomial, data = tab)
fm2a<- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Changevector_mean_5max + Age +(1|Replicate), family = binomial, data = tab)
fm1c <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + Temp_mean + Age+(1|Replicate), family = binomial, data = tab)
fm1d <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Age + RH._Max + (1|Replicate), family = binomial, data = tab)
fm1e <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Age + Temp_mean + (1|Replicate), family = binomial, data = tab)
AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm1a,fm1b,fm1c,fm1d,fm1e,fm2a)
summary(fm2)


##models for threshold 10

fm11 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + Temp_mean + Duration.ms._10max + changescalar_max_10max+Angle_max_10max+ (1|Replicate), family = binomial, data = tab)
fm12 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + Temp_mean + Duration.ms._10max + changescalar_max_10max+ (1|Replicate), family = binomial, data = tab)
fm13 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Temp_mean +RH._Max + Duration.ms._10mean +  Changevector_mean_10mean +  changescalar_max_10max+(1|Replicate), family = binomial, data = tab)
fm14 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max +Duration.ms._10max + changescalar_mean_10max+Changevector_mean_10mean + (1|Replicate), family = binomial, data = tab)
fm15 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Temp_mean +Duration.ms._10mean + changescalar_mean_10max+ (1|Replicate), family = binomial, data = tab)
fm11a <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max +Temp_mean + Duration.ms._10max + Changevector_mean_10mean +changescalar_max_10max+(1|Replicate), family = binomial, data = tab)
fm11b <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + Temp_mean + Duration.ms._10max + Changevector_max_10mean + Angle_max_10max+ (1|Replicate), family = binomial, data = tab)
fm11c <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + Temp_mean + Duration.ms._10mean + Changevector_mean_10mean +(1|Replicate), family = binomial, data = tab)
fm11d <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Temp_mean +Duration.ms._10mean +Changevector_mean_10mean + Angle_max_10max+ (1|Replicate), family = binomial, data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm1a,fm1b,fm1c,fm1d,fm1e,fm2a,fm11,fm12,fm13,fm14,fm15,fm11a,fm11b,fm11c,fm11d)
summary(fm2)
#

##threshold 15
fm21 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Temp_mean + changescalar_mean_15max + Duration.ms._15mean + Age+(1|Replicate), family = binomial, data = tab)
fm22 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Temp_mean + changescalar_mean_15max +Duration.ms._15max +(1|Replicate), family = binomial, data = tab)
fm23 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + changescalar_mean_15max + Age +(1|Replicate), family = binomial, data = tab)
fm24 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + changescalar_max_15mean  + Angle_max_15max+ Angle_mean_15max+Duration.ms._15mean +(1|Replicate), family = binomial, data = tab)
fm25 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + changescalar_mean_15max + Angle_mean_15max+ (1|Replicate), family = binomial, data = tab)
fm26 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + changescalar_mean_15max + Duration.ms._15mean +(1|Replicate), family = binomial, data = tab)
fm27 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + changescalar_max_15mean  + Age + (1|Replicate), family = binomial, data = tab)
fm28 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Temp_mean +  changescalar_mean_15mean +(1|Replicate), family = binomial, data = tab)
fm29 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + changescalar_mean_15max + Changevector_mean_15max+ (1|Replicate), family = binomial, data = tab)
fm30 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation +Temp_mean + Duration.ms._15max + Age+(1|Replicate), family = binomial, data = tab)
fm30a <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation +Temp_mean + Duration.ms._15max +changescalar_mean_15mean + Angle_mean_15max+ (1|Replicate), family = binomial, data = tab)
fm29a<- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + changescalar_max_15mean + Duration.ms._15mean +(1|Replicate), family = binomial, data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm1a,fm1b,fm1c,fm1d,fm1e,fm2a,fm11,fm12,fm13,fm14,fm15,fm11a,fm11b,fm11c,fm11d,fm21,fm22,fm23,fm24,fm25,fm26,fm27,fm28,fm29,fm29a,fm30,fm30a)
summary(fm2)

##Threshold 20

fm31 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Temp_mean + Angle_mean_20max + scalar_max_20mean +(1|Replicate), family = binomial, data = tab)
fm32 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Temp_mean + Changevector_max_20max + scalar_max_20mean +(1|Replicate), family = binomial, data = tab)
fm31a <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Age + Temp_mean +scalar_max_20mean + (1|Replicate), family = binomial, data = tab)
fm32a <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Temp_mean + Angle_mean_20max + Changevector_max_20max + (1|Replicate), family = binomial, data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm1a,fm1b,fm1c,fm1d,fm1e,fm2a,fm11,fm12,fm13,fm14,fm15,fm11a,fm11b,fm11c,fm11d,fm21,fm22,fm23,fm24,fm25,fm26,fm27,fm28,fm29,fm29a,fm30,fm30a,fm31,fm32,fm31a,fm32a)
summary(fm2)


#combination of threseholds

fm41 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Temp_mean +Duration.ms._10mean + changescalar_mean_10max+Changevector_mean_5max + Age +(1|Replicate), family = binomial, data = tab)
fm42 <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + changescalar_mean_10max+ (1|Replicate), family = binomial, data = tab)
fm41a <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Temp_mean +Duration.ms._10mean + Changevector_mean_5max +(1|Replicate), family = binomial, data = tab)
fm42a <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + changescalar_mean_10max+ (1|Replicate), family = binomial, data = tab)
fm41b <- glmer(cbind(pairs_formed,unformed_pairs) ~ Irradiation + Temp_mean +Duration.ms._10mean + changescalar_mean_10max+Changevector_mean_5max +(1|Replicate), family = binomial, data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm1a,fm1b,fm1c,fm1d,fm1e,fm2a,fm11,fm12,fm13,fm14,fm15,fm11a,fm11b,fm11c,fm11d,fm21,fm22,fm23,fm24,fm25,fm26,fm27,fm28,fm29,fm29a,fm30,fm30a,fm31,fm32,fm31a,fm32a,fm41,fm41b,fm42,fm41a,fm42a)

#the best model remains as fm2 


##Figure 9c

head(tab)
str(tab)
summary(fm2)
plot((pairs_formed / unformed_pairs) ~ fitted(fm2), data = tab)
abline(lm((pairs_formed / unformed_pairs) ~ fitted(fm2), data = tab), col = "red")
cor.test(tab$pairs_formed/tab$unformed_pairs,fitted(fm2))
summary(lm(fitted(fm2)~(tab$pairs_formed/tab$unformed_pairs)))$r.squared

########Correlation of factors from best model(fm12)
cor.test(tab$changescalar_mean_5max,tab$Changevector_mean_5max)
plot(tab$changescalar_mean_5max~tab$Changevector_mean_5max)
## Conclusion: Correlated.One can be dropped 


#######*************************************************************************************************************************************************************

```




THE IMPACT OF SHOCK ON THE INSEMINATION RATE OF 22 AND 29 DAY OLD PUPAE (COMBINED DATA)




Figure 9d:INSEMINATION RATE



```{r}



tab <- read.csv("Figure 9d.csv")
head(tab)
summary(tab)
str(tab)
tab$pcinsemination <- tab$full / (tab$full+tab$not_full)
boxplot(tab$pcinsemination ~ tab$Treatments, xlab ="Treatments", ylab = "Insemination rate")


#####SCATTER PLOTS

plot(tab$pcinsemination~tab$RH._Max)
plot(tab$pcinsemination~tab$RH._mean)#*
plot(tab$pcinsemination~tab$Temp_mean)#*
cor.test(tab$pcinsemination,tab$Temp_mean)
plot(tab$pcinsemination~tab$Temp_Max)

par(mfrow = c(2,2))

#Scatter plots at threshold 5-Max
plot(tab$pcinsemination~tab$Events._5max)#*
plot(tab$pcinsemination~tab$Duration.ms._5max)
plot(tab$pcinsemination~tab$scalar_max_5max)
plot(tab$pcinsemination~tab$scalar_mean_5max)
plot(tab$pcinsemination~tab$changescalar_max_5max)
plot(tab$pcinsemination~tab$changescalar_mean_5max)#*
plot(tab$pcinsemination~tab$Changevector_max_5max)    
plot(tab$pcinsemination~tab$Changevector_mean_5max)#*
plot(tab$pcinsemination~tab$Angle_max_5max)        
plot(tab$pcinsemination~tab$Angle_mean_5max)
plot(tab$pcinsemination~tab$Hz_max_5max)    
plot(tab$pcinsemination~tab$Hz_mean_5max)
```


#Scatter plots at threshold 5-Mean
```{r }
par(mfrow = c(2,2))
plot(tab$pcinsemination~tab$Duration.ms._5mean)#*
plot(tab$pcinsemination~tab$scalar_max_5mean)
plot(tab$pcinsemination~tab$scalar_mean_5mean)
plot(tab$pcinsemination~tab$changescalar_max_5mean)
plot(tab$pcinsemination~tab$changescalar_mean_5mean)
plot(tab$pcinsemination~tab$Changevector_max_5mean)    
plot(tab$pcinsemination~tab$Changevector_mean_5mean)
plot(tab$pcinsemination~tab$Angle_max_5mean)        
plot(tab$pcinsemination~tab$Angle_mean_5mean)
plot(tab$pcinsemination~tab$Hz_max_5mean)    
plot(tab$pcinsemination~tab$Hz_mean_5mean)
```


#Scatter plots at threshold 10-Max
```{r }
par(mfrow = c(2,2))
plot(tab$pcinsemination~tab$Duration.ms._10max)#
plot(tab$pcinsemination~tab$Duration.ms._10max)
plot(tab$pcinsemination~tab$scalar_max_10max)
plot(tab$pcinsemination~tab$scalar_mean_10max)#*
plot(tab$pcinsemination~tab$changescalar_max_10max)
plot(tab$pcinsemination~tab$changescalar_mean_10max)#*
cor.test(tab$pcinsemination,tab$changescalar_mean_10max)
plot(tab$pcinsemination~tab$Changevector_max_10max)    
plot(tab$pcinsemination~tab$Changevector_mean_10max)#*
plot(tab$pcinsemination~tab$Angle_max_10max)        
plot(tab$pcinsemination~tab$Angle_mean_10max)#*
plot(tab$pcinsemination~tab$Hz_max_10max)    
plot(tab$pcinsemination~tab$Hz_mean_10max)#*
```



#Scatter plots at threshold 10-Mean
```{r }
par(mfrow = c(2,2))
plot(tab$pcinsemination~tab$Duration.ms._10mean)#*
cor.test(tab$pcinsemination,tab$Duration.ms._10mean)
plot(tab$pcinsemination~tab$scalar_max_10mean)
plot(tab$pcinsemination~tab$scalar_mean_10mean)
plot(tab$pcinsemination~tab$changescalar_max_10mean)#
plot(tab$pcinsemination~tab$changescalar_mean_10mean)
plot(tab$pcinsemination~tab$Changevector_max_10mean)    
plot(tab$pcinsemination~tab$Changevector_mean_10mean)#
plot(tab$pcinsemination~tab$Angle_max_10mean)        
plot(tab$pcinsemination~tab$Angle_mean_10mean)
plot(tab$pcinsemination~tab$Hz_max_10mean)    
plot(tab$pcinsemination~tab$Hz_mean_10mean)
```
#Scatter plots at threshold 15-Max

```{r }

par(mfrow = c(2,2))
plot(tab$pcinsemination~tab$Events_15max)
plot(tab$pcinsemination~tab$Duration.ms._15max)#*
plot(tab$pcinsemination~tab$scalar_max_15max)
plot(tab$pcinsemination~tab$scalar_mean_15max)
plot(tab$pcinsemination~tab$changescalar_max_15max)#*
plot(tab$pcinsemination~tab$changescalar_mean_15max)
plot(tab$pcinsemination~tab$Changevector_max_15max)    
plot(tab$pcinsemination~tab$Changevector_mean_15max)
plot(tab$pcinsemination~tab$Angle_max_15max)#*        
plot(tab$pcinsemination~tab$Angle_mean_15max)
plot(tab$pcinsemination~tab$Hz_max_15max) #*   
plot(tab$pcinsemination~tab$Hz_mean_15max)
```

#Scatter plots at threshold 15-Mean

```{r }
par(mfrow = c(2,2))
plot(tab$pcinsemination~tab$Duration.ms._15mean)#*
plot(tab$pcinsemination~tab$scalar_max_15mean)#*
plot(tab$pcinsemination~tab$scalar_mean_15mean)
plot(tab$pcinsemination~tab$changescalar_max_15mean)#*
plot(tab$pcinsemination~tab$changescalar_mean_15mean)
plot(tab$pcinsemination~tab$Changevector_max_15mean)#*   
plot(tab$pcinsemination~tab$Changevector_mean_15mean)
plot(tab$pcinsemination~tab$Angle_max_15mean)        
plot(tab$pcinsemination~tab$Angle_mean_15mean)
plot(tab$pcinsemination~tab$Hz_max_15mean)    
plot(tab$pcinsemination~tab$Hz_mean_15mean)
```


#Scatter plots at threshold 20-Max
```{r }
par(mfrow = c(2,2))
plot(tab$pcinsemination~tab$Events_20max)
plot(tab$pcinsemination~tab$Duration.ms._20max)#*
plot(tab$pcinsemination~tab$scalar_max_20max)
plot(tab$pcinsemination~tab$scalar_mean_20max)
plot(tab$pcinsemination~tab$changescalar_max_20max)#*
plot(tab$pcinsemination~tab$changescalar_mean_20max)
plot(tab$pcinsemination~tab$Changevector_max_20max)    
plot(tab$pcinsemination~tab$Changevector_mean_20max)
plot(tab$pcinsemination~tab$Angle_max_20max)#*        
plot(tab$pcinsemination~tab$Angle_mean_20max)#*
plot(tab$pcinsemination~tab$Hz_max_20max)#*    
plot(tab$pcinsemination~tab$Hz_mean_20max)#*
```


#Scatter plots at threshold 20-Mean

```{r }
par(mfrow = c(2,2))
plot(tab$pcinsemination~tab$Events_20mean)#*
plot(tab$pcinsemination~tab$Duration.ms._20mean)#*
plot(tab$pcinsemination~tab$scalar_max_20mean)
plot(tab$pcinsemination~tab$scalar_mean_20mean)
plot(tab$pcinsemination~tab$changescalar_max_20mean)
plot(tab$pcinsemination~tab$changescalar_mean_20mean)
plot(tab$pcinsemination~tab$Changevector_max_20mean)#* 
plot(tab$pcinsemination~tab$Changevector_mean_20mean)
plot(tab$pcinsemination~tab$Angle_max_20mean)        
plot(tab$pcinsemination~tab$Angle_mean_20mean)
plot(tab$pcinsemination~tab$Hz_max_20mean)#*    
plot(tab$pcinsemination~tab$Hz_mean_20mean)
```
# Models of Impact of transportation on the insemination rate

```{r }
tab$Age <- as.factor(tab$Age)


#models threshold 5
fm1 <- glmer(cbind(full,not_full) ~ Irradiation +  RH._Max +Temp_mean + (1|Replicate), family = binomial, data = tab)
fm2 <- glmer(cbind(full,not_full) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Duration.ms._5mean + (1|Replicate), family = binomial, data = tab)
fm3 <- glmer(cbind(full,not_full) ~ Irradiation + RH._Max + Age + Duration.ms._5mean +(1|Replicate), family = binomial, data = tab)
fm4 <- glmer(cbind(full,not_full) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Changevector_max_5max +   (1|Replicate), family = binomial, data = tab)
fm5 <- glmer(cbind(full,not_full) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Duration.ms._5mean + Changevector_max_5max +(1|Replicate), family = binomial, data = tab)
fm6 <- glmer(cbind(full,not_full) ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + (1|Replicate), family = binomial, data = tab)
fm7 <- glmer(cbind(full,not_full) ~ Irradiation + RH._Max + Temp_mean + Changevector_max_5max +  (1|Replicate), family = binomial, data = tab)
fm8 <- glmer(cbind(full,not_full) ~ Irradiation + Temp_mean + changescalar_mean_5max + scalar_max_5mean + (1|Replicate), family = binomial, data = tab)
fm1a <- glmer(cbind(full,not_full) ~ Irradiation + Age+ Temp_mean + (1|Replicate), family = binomial, data = tab)
fm1b <- glmer(cbind(full,not_full) ~ Irradiation + Age+ RH._Max +  (1|Replicate), family = binomial, data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm1a,fm1b)
summary(fm3)

fm11 <- glmer(cbind(full,not_full) ~ Irradiation +  RH._Max + Duration.ms._10max + changescalar_mean_10max+(1|Replicate), family = binomial, data = tab)
fm12 <- glmer(cbind(full,not_full) ~ Irradiation + Temp_mean + Duration.ms._10mean + changescalar_mean_10max+ (1|Replicate), family = binomial, data = tab)
fm13 <- glmer(cbind(full,not_full) ~ Irradiation + Temp_mean + Duration.ms._10max + Changevector_mean_10max + changescalar_max_10mean+(1|Replicate), family = binomial, data = tab)
fm14 <- glmer(cbind(full,not_full) ~ Irradiation + Duration.ms._10mean + changescalar_mean_10max+ Changevector_max_10mean +(1|Replicate), family = binomial, data = tab)
fm15 <- glmer(cbind(full,not_full) ~ Irradiation + Angle_mean_15max+Duration.ms._10max + changescalar_mean_10max+ (1|Replicate), family = binomial, data = tab)
fm11a <- glmer(cbind(full,not_full) ~ Irradiation + Age + RH._Max + Duration.ms._10max + changescalar_mean_10max+(1|Replicate), family = binomial, data = tab)
fm12a <- glmer(cbind(full,not_full) ~ Irradiation + RH._Max + Duration.ms._10mean + changescalar_mean_10max+ (1|Replicate), family = binomial, data = tab)
fm11b <- glmer(cbind(full,not_full) ~ Irradiation +  RH._Max + Duration.ms._10max + scalar_mean_10max+(1|Replicate), family = binomial, data = tab)
fm11c <- glmer(cbind(full,not_full) ~ Irradiation +  RH._Max + Duration.ms._10max + changescalar_mean_10max+(1|Replicate), family = binomial, data = tab)
fm12b <- glmer(cbind(full,not_full) ~ Irradiation + Temp_mean + Duration.ms._10mean + (1|Replicate), family = binomial, data = tab)


AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm1a,fm1b,fm11,fm12,fm13,fm14,fm11a,fm11b,fm11c,fm12a,fm12b)
summary(fm12)

#threshold 15
fm21 <- glmer(cbind(full,not_full) ~ Irradiation + Temp_mean + Angle_max_15max+changescalar_max_15mean+Duration.ms._15max+(1|Replicate), family = binomial, data = tab)
fm22 <- glmer(cbind(full,not_full) ~ Irradiation + RH._mean + changescalar_max_15mean+Duration.ms._15mean+ (1|Replicate), family = binomial, data = tab)
fm23 <- glmer(cbind(full,not_full) ~ Irradiation +  Angle_max_15max+Duration.ms._15max+changescalar_max_15mean+(1|Replicate), family = binomial, data = tab)
fm24 <- glmer(cbind(full,not_full) ~ Irradiation + changescalar_max_15mean + Angle_max_15max+Duration.ms._15mean+ (1|Replicate), family = binomial, data = tab)
fm25 <- glmer(cbind(full,not_full) ~ Irradiation + changescalar_max_15mean +  Angle_max_15max+ scalar_max_10mean+(1|Replicate), family = binomial, data = tab)
fm26 <- glmer(cbind(full,not_full) ~ Irradiation + Angle_max_15max+scalar_max_15mean + Duration.ms._15max +(1|Replicate), family = binomial, data = tab)
fm27 <- glmer(cbind(full,not_full) ~ Irradiation + Angle_max_15max+Changevector_max_15mean+scalar_max_10mean+ (1|Replicate), family = binomial, data = tab)
fm28 <- glmer(cbind(full,not_full) ~ Irradiation + Temp_mean + Changevector_max_15mean+scalar_max_15mean+(1|Replicate), family = binomial, data = tab)
fm29 <- glmer(cbind(full,not_full) ~ Irradiation + changescalar_mean_15max + Changevector_mean_15mean+ Age+ (1|Replicate), family = binomial, data = tab)
fm30 <- glmer(cbind(full,not_full) ~ Irradiation + scalar_max_15mean+ Duration.ms._15mean+ Changevector_max_15mean+(1|Replicate), family = binomial, data = tab)
fm21a <- glmer(cbind(full,not_full) ~ Irradiation + Temp_mean + Angle_max_15max+changescalar_max_15mean+ (1|Replicate), family = binomial, data = tab)
fm27a <- glmer(cbind(full,not_full) ~ Irradiation +Temp_mean + Duration.ms._15mean+ Changevector_max_15mean+ (1|Replicate), family = binomial, data = tab)
fm21b <- glmer(cbind(full,not_full) ~ Irradiation + Temp_mean + Angle_max_15max+changescalar_max_15mean+ (1|Replicate), family = binomial, data = tab)
AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm1a,fm1b,fm11,fm12,fm13,fm14,fm11a,fm11b,fm11c,fm12a,fm12b,fm21,fm22,fm23,fm24,fm25,fm26,fm27,fm28,fm29,fm30,fm21a,fm27a,fm21b)
summary(fm12)

fm31 <- glmer(cbind(full,not_full) ~ Irradiation + Temp_mean +Changevector_max_20max+(1|Replicate), family = binomial, data = tab)
fm32 <- glmer(cbind(full,not_full) ~ Irradiation + Temp_mean + Changevector_max_20mean+ Duration.ms._20max +(1|Replicate), family = binomial, data = tab)
fm31a <- glmer(cbind(full,not_full) ~ Irradiation + RH._mean + Angle_max_20max +(1|Replicate), family = binomial, data = tab)
fm31b <- glmer(cbind(full,not_full) ~ Irradiation + Temp_mean +Duration.ms._20max +(1|Replicate), family = binomial, data = tab)
fm32a <- glmer(cbind(full,not_full) ~ Irradiation + Temp_mean + Changevector_max_20mean+ Duration.ms._20mean +(1|Replicate), family = binomial, data = tab)
fm31c <- glmer(cbind(full,not_full) ~ Irradiation + Temp_mean +Changevector_max_20max+Angle_mean_20max +(1|Replicate), family = binomial, data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm1a,fm1b,fm11,fm12,fm13,fm14,fm11a,fm11b,fm11c,fm12a,fm12b,fm21,fm22,fm23,fm24,fm25,fm26,fm27,fm28,fm29,fm30,fm21a,fm27a,fm21b,fm31,fm32,fm31a,fm31b,fm32a,fm31c)
summary(fm12)

#combination of threseholds

fm41 <- glmer(cbind(full,not_full) ~ Irradiation + RH._Max + changescalar_mean_10max+ Duration.ms._5mean +Angle_max_20max +(1|Replicate), family = binomial, data = tab)
fm42 <- glmer(cbind(full,not_full) ~ Irradiation + Temp_mean + Duration.ms._10mean + Angle_max_20max ++ (1|Replicate), family = binomial, data = tab)
fm43 <- glmer(cbind(full,not_full) ~ Irradiation + RH._mean + changescalar_max_15mean+Duration.ms._5mean+ Angle_max_20max +(1|Replicate), family = binomial, data = tab)
fm44 <- glmer(cbind(full,not_full) ~ Irradiation + RH._mean + Angle_max_20max +(1|Replicate), family = binomial, data = tab)
fm45 <- glmer(cbind(full,not_full) ~ Irradiation + RH._Max + Age + changescalar_mean_10max+(1|Replicate), family = binomial, data = tab)
fm46 <- glmer(cbind(full,not_full) ~ Irradiation + RH._Max + Age + changescalar_mean_10max+ Duration.ms._5mean +(1|Replicate), family = binomial, data = tab)
fm47 <- glmer(cbind(full,not_full) ~ Irradiation + RH._mean + changescalar_max_15mean+Duration.ms._5mean+ (1|Replicate), family = binomial, data = tab)
AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm1a,fm1b,fm11,fm12,fm13,fm14,fm11a,fm11b,fm11c,fm12a,fm12b,fm21,fm22,fm23,fm24,fm25,fm26,fm27,fm28,fm29,fm30,fm21a,fm27a,fm21b,fm31,fm32,fm31a,fm31b,fm32a,fm31c,fm41,fm42,fm43,fm44,fm45,fm46,fm47)
summary(fm12)                

#the best model remains as fm12 


##Figure 9d

head(tab)
str(tab)
summary(fm12)
plot((full / not_full) ~ fitted(fm12), data = tab)
abline(lm((full / not_full) ~ fitted(fm12), data = tab), col = "red")
cor.test(tab$full/tab$not_full,fitted(fm12))
summary(lm(fitted(fm12)~(tab$full/tab$not_full)))$r.squared

##### Independent variables not correlated

```







THE IMPACT OF SHOCK ON THE MEAN SPERMATHECAL VALUE(MSV) OF 22 AND 29 DAY OLD PUPAE(COMBINED DATA)
Figure 9e


```{r}



tab <- read.csv("Figure 9e.csv")
head(tab)
summary(tab)
str(tab)
boxplot(tab$MSV ~ tab$Treatments, xlab ="Treatments", ylab = "Mean Spermathecal Value(MSV)")
na.omit(tab)




#####SCATTER PLOTS
plot(tab$MSV~tab$RH._Max)#*
plot(tab$MSV~tab$RH._mean)#*
plot(tab$MSV~tab$Temp_mean)#*
plot(tab$MSV~tab$Temp_Max)


#Scatter plots at threshold 5-Max

par(mfrow = c(2,2))
plot(tab$MSV~tab$Events._5max)
plot(tab$MSV~tab$Duration.ms._5max)
plot(tab$MSV~tab$scalar_max_5max)
plot(tab$MSV~tab$scalar_mean_5max)
plot(tab$MSV~tab$changescalar_max_5max)
plot(tab$MSV~tab$changescalar_mean_5max)#*
plot(tab$MSV~tab$Changevector_max_5max)    
plot(tab$MSV~tab$Changevector_mean_5max)#*
plot(tab$MSV~tab$Angle_max_5max)        
plot(tab$MSV~tab$Angle_mean_5max)
plot(tab$MSV~tab$Hz_max_5max)    
plot(tab$MSV~tab$Hz_mean_5max)
```

#Scatter plots at threshold 5-Mean

```{r }

par(mfrow = c(2,2))
plot(tab$MSV~tab$Duration.ms._5mean)#*
plot(tab$MSV~tab$scalar_max_5mean)
plot(tab$MSV~tab$scalar_mean_5mean)
plot(tab$MSV~tab$changescalar_max_5mean)
plot(tab$MSV~tab$changescalar_mean_5mean)
plot(tab$MSV~tab$Changevector_max_5mean)    
plot(tab$MSV~tab$Changevector_mean_5mean)
plot(tab$MSV~tab$Angle_max_5mean)        
plot(tab$MSV~tab$Angle_mean_5mean)
plot(tab$MSV~tab$Hz_max_5mean)    
plot(tab$MSV~tab$Hz_mean_5mean)
```


#Scatter plots at threshold 10-Max

```{r }

par(mfrow = c(2,2))
plot(tab$MSV~tab$Events_10max)
plot(tab$MSV~tab$Duration.ms._10max)
plot(tab$MSV~tab$scalar_max_10max)
plot(tab$MSV~tab$scalar_mean_10max)#*
plot(tab$MSV~tab$changescalar_max_10max)
plot(tab$MSV~tab$changescalar_mean_10max)#*
plot(tab$MSV~tab$Changevector_max_10max)    
plot(tab$MSV~tab$Changevector_mean_10max)#*
plot(tab$MSV~tab$Angle_max_10max)        
plot(tab$MSV~tab$Angle_mean_10max)#*
plot(tab$MSV~tab$Hz_max_10max)    
plot(tab$MSV~tab$Hz_mean_10max)#*
```

#Scatter plots at threshold 10-Mean
```{r }

par(mfrow = c(2,2))
plot(tab$MSV~tab$Duration.ms._10mean)#*
plot(tab$MSV~tab$scalar_max_10mean)
plot(tab$MSV~tab$scalar_mean_10mean)
plot(tab$MSV~tab$changescalar_max_10mean)
plot(tab$MSV~tab$changescalar_mean_10mean)
plot(tab$MSV~tab$Changevector_max_10mean)    
plot(tab$MSV~tab$Changevector_mean_10mean)
plot(tab$MSV~tab$Angle_max_10mean)        
plot(tab$MSV~tab$Angle_mean_10mean)
plot(tab$MSV~tab$Hz_max_10mean)    
plot(tab$MSV~tab$Hz_mean_10mean)
```


#Scatter plots at threshold 15-Max
```{r }

par(mfrow = c(2,2))
plot(tab$MSV~tab$Events_15max)
plot(tab$MSV~tab$Duration.ms._15max)
plot(tab$MSV~tab$scalar_max_15mean )#*
plot(tab$MSV~tab$scalar_mean_15max)
plot(tab$MSV~tab$changescalar_max_15max)
plot(tab$MSV~tab$changescalar_mean_15max)
plot(tab$MSV~tab$Changevector_max_15max)    
plot(tab$MSV~tab$Changevector_mean_15max)
plot(tab$MSV~tab$Angle_max_15max)#*        
plot(tab$MSV~tab$Angle_mean_15max)
plot(tab$MSV~tab$Hz_max_15max) #*   
plot(tab$MSV~tab$Hz_mean_15max)
```

#Scatter plots at threshold 15-Mean
```{r }

par(mfrow = c(2,2))
plot(tab$MSV~tab$Duration.ms._15mean)#*
cor.test(tab$MSV,tab$Duration.ms._15mean)
plot(tab$MSV~tab$scalar_max_15mean)
plot(tab$MSV~tab$scalar_mean_15mean)
plot(tab$MSV~tab$changescalar_max_15mean)#*
plot(tab$MSV~tab$changescalar_mean_15mean)
plot(tab$MSV~tab$Changevector_max_15mean)   
plot(tab$MSV~tab$Changevector_mean_15mean)
plot(tab$MSV~tab$Angle_max_15mean)        
plot(tab$MSV~tab$Angle_mean_15mean)
plot(tab$MSV~tab$Hz_max_15mean)    
plot(tab$MSV~tab$Hz_mean_15mean)
```

#Scatter plots at threshold 20-Max

```{r }

par(mfrow = c(2,2))
plot(tab$MSV~tab$Events_20mean)
plot(tab$MSV~tab$Duration.ms._20max)#*
plot(tab$MSV~tab$scalar_max_20max)
plot(tab$MSV~tab$scalar_mean_20max)#*
plot(tab$MSV~tab$changescalar_max_20max)
plot(tab$MSV~tab$changescalar_mean_20max)
plot(tab$MSV~tab$Changevector_max_20max)#*    
plot(tab$MSV~tab$Changevector_mean_20max)
plot(tab$MSV~tab$Angle_max_20max)#*        
plot(tab$MSV~tab$Angle_mean_20max)
plot(tab$MSV~tab$Hz_max_20max)#*    
plot(tab$MSV~tab$Hz_mean_20max)#*
```

#Scatter plots at threshold 20-Mean


```{r }

par(mfrow = c(2,2))
plot(tab$MSV~tab$Duration.ms._20mean)#*
plot(tab$MSV~tab$scalar_max_20mean)
plot(tab$MSV~tab$scalar_mean_20mean)
plot(tab$MSV~tab$changescalar_max_20mean)
plot(tab$MSV~tab$changescalar_mean_20mean)
plot(tab$MSV~tab$Changevector_max_20mean)#* 
plot(tab$MSV~tab$Changevector_mean_20mean)
plot(tab$MSV~tab$Angle_max_20mean)        
plot(tab$MSV~tab$Angle_mean_20mean)
plot(tab$MSV~tab$Hz_max_20mean)#*    
plot(tab$MSV~tab$Hz_mean_20mean)
```

#Models of the impact of shock during transportation, on the Mean Spermathecal Value (MSV)



```{r }

#models threshold 5

fm1 <- lme(MSV  ~ Irradiation + RH._Max + Temp_mean, random=~1|Replicate, data = tab)
fm2 <- lme(MSV  ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Changevector_mean_5max, random=~1|Replicate,  data = tab)
fm3 <- lme(MSV  ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max +Angle_mean_5max ,  random=~1|Replicate,  data = tab)
fm4 <- lme(MSV  ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max +Duration.ms._5mean, random=~1|Replicate, data = tab)
fm5 <- lme(MSV  ~ Irradiation + RH._Max + Temp_mean + Duration.ms._5mean, random=~1|Replicate,  data = tab)
fm6 <- lme(MSV  ~ Irradiation + RH._Max + Temp_mean + changescalar_mean_5max + Angle_mean_5max + Duration.ms._5mean, random=~1|Replicate,  data = tab)
fm7 <- lme(MSV  ~ Irradiation + RH._Max + Temp_mean +Changevector_mean_5max,random=~1|Replicate,  data = tab)
fm8 <- lme(MSV  ~ Irradiation + Temp_mean + changescalar_mean_5max + Duration.ms._5max +Angle_mean_5max, random=~1|Replicate, data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8)
summary(fm1)


##models for threshold 10

fm11 <- lme(MSV  ~ Irradiation + Temp_mean +  changescalar_mean_10max+Angle_mean_10max, random=~1|Replicate, data = tab)
fm12 <- lme(MSV  ~ Irradiation + Temp_mean + Duration.ms._10mean + Changevector_mean_10max + Hz_mean_10max , random=~1|Replicate,  data = tab)
fm13 <- lme(MSV  ~ Irradiation + Temp_mean + Duration.ms._10max + Changevector_mean_10max + Hz_mean_10max , random=~1|Replicate, data = tab)
fm14 <- lme(MSV  ~ Irradiation + Duration.ms._10mean + changescalar_mean_10max+ Hz_mean_10max , random=~1|Replicate,  data = tab)
fm15 <- lme(MSV  ~ Irradiation + Temp_mean +Angle_mean_10max + changescalar_mean_10max, random=~1|Replicate,  data = tab)
fm11a <- lme(MSV  ~ Irradiation + Temp_mean + Duration.ms._10max + Changevector_mean_10max , random=~1|Replicate,  data = tab)
fm11b <- lme(MSV  ~ Irradiation + Temp_mean + Duration.ms._10mean + changescalar_mean_10max, random=~1|Replicate,  data = tab)
fm11c <- lme(MSV  ~ Irradiation + Temp_mean + Duration.ms._10mean , random=~1|Replicate, data = tab)
fm11d <- lme(MSV  ~ Irradiation + Temp_mean +Duration.ms._10mean +Hz_mean_10max, random=~1|Replicate,  data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm11,fm12,fm13,fm14,fm15,fm11a,fm11b,fm11c,fm11d)
summary(fm11c)
#

##threshold 15

fm21 <-lme(MSV  ~ Irradiation + Temp_mean + Duration.ms._15mean , random=~1|Replicate, data = tab)
fm22 <- lme(MSV  ~ Irradiation + Temp_mean +Angle_mean_15max, random=~1|Replicate, data = tab)
fm23 <- lme(MSV  ~ Irradiation + Temp_mean +Hz_max_15max , random=~1|Replicate,data = tab)
fm24 <-lme(MSV  ~ Irradiation + Temp_mean +Hz_max_15max +Duration.ms._15mean , random=~1|Replicate, data = tab)
fm25 <- lme(MSV  ~ Irradiation +Temp_mean + Duration.ms._15mean  +Angle_mean_15max, random=~1|Replicate, data = tab)
fm26 <- lme(MSV  ~ Irradiation +Temp_mean, random=~1|Replicate, data = tab)
fm27 <- lme(MSV  ~ Irradiation +Angle_mean_15max, random=~1|Replicate, data = tab)
fm28 <- lme(MSV  ~ Irradiation + Duration.ms._15mean , random=~1|Replicate, data = tab)
fm29 <- lme(MSV  ~ Irradiation +Angle_mean_15max+  Duration.ms._15mean , random=~1|Replicate, data = tab)
fm30 <- lme(MSV  ~ Irradiation +Temp_mean +Angle_mean_15max , random=~1|Replicate, data = tab)
fm30a <- lme(MSV  ~ Irradiation +Temp_mean + Duration.ms._15mean , random=~1|Replicate,data = tab)
fm29a<- lme(MSV  ~ Irradiation + Duration.ms._15mean +scalar_mean_15mean, random=~1|Replicate, data = tab)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm13,fm14,fm15,fm11a,fm11b,fm11c,fm11d,fm21,fm22,fm23,fm24,fm25,fm26,fm27,fm28,fm29,fm30,fm30a,fm29a)
summary(fm28)


##Threshold 20

fm31 <- lme(MSV  ~ Irradiation + Temp_mean + Duration.ms._20max + Angle_mean_20max +Hz_max_20mean, random=~1|Replicate,data = tab)
fm32 <- lme(MSV  ~ Irradiation + Temp_mean +Duration.ms._20mean +Hz_max_20mean, random=~1|Replicate, data = tab)
fm31a <-lme(MSV  ~ Irradiation + Temp_mean + Changevector_max_20mean + Hz_max_20max , random=~1|Replicate, data = tab)
fm32a <- lme(MSV  ~ Irradiation + Temp_mean +Changevector_max_20mean +Angle_mean_20max , random=~1|Replicate,data = tab)
fm32b <- lme(MSV  ~ Irradiation + Temp_mean +  Changevector_max_20max +Duration.ms._20mean, random=~1|Replicate,data = tab)
fm31b <- lme(MSV  ~ Irradiation + Temp_mean + Changevector_max_20mean + Hz_max_20max , random=~1|Replicate, data = tab)
fm32c <- lme(MSV  ~ Irradiation + Temp_mean + Hz_max_20max , random=~1|Replicate, data = tab)


AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm13,fm14,fm15,fm11a,fm11b,fm11c,fm11d,fm21,fm22,fm23,fm24,fm25,fm26,fm27,fm28,fm29,fm30,fm30a,fm29a,fm31,fm32,fm31a,fm32a,fm31b,fm32b,fm32c)
summary(fm28)


#combination of threseholds
fm41 <- lme(MSV  ~ Irradiation + Temp_mean + Duration.ms._10mean + changescalar_mean_10max+Angle_mean_15max , random=~1|Replicate, data = tab)
fm42 <- lme(MSV  ~ Irradiation + Temp_mean + Duration.ms._10mean + changescalar_mean_10max+Hz_mean_10max + Duration.ms._15mean , random=~1|Replicate, data = tab)
summary(fm28)

AICc(fm1,fm2,fm3,fm4,fm5,fm6,fm7,fm8,fm13,fm14,fm15,fm11a,fm11b,fm11c,fm11d,fm21,fm22,fm23,fm24,fm25,fm26,fm27,fm28,fm29,fm30,fm30a,fm29a,fm31,fm32,fm31a,fm32a,fm31b,fm32b,fm32c,fm41,fm42)
summary(fm28)

#the best model remains as fm28 

##Fig 9e

head(tab)
str(tab)
summary(fm28)
plot(MSV ~ fitted(fm28), data = tab)
abline(lm(MSV ~ fitted(fm28), data = tab), col = "red")
cor.test(tab$MSV,fitted(fm28))
summary(lm(fitted(fm28)~(tab$MSV)))$r.squared



```





ALL ANALYSIS OF SUPPORTING INFORMATION




Analysis of the statistical relationships of temperature and humidity in/between 22 and 29 day old pupae 


```{r }


###### Correlation of temperature and humidity at age 22, 29
tab2 <- read.csv("Temp_RH_22.csv")
head(tab2)
plot(tab2$Humidity,tab2$Temperature,xlab = "Humidity", ylab = " Temperature (%)")
cor.test(tab2$Humidity,tab2$Temperature, method = "pearson")


tab2 <- read.csv("Temp_RH_29_2.csv")
head(tab2)
plot(tab2$Humidity,tab2$Temperature,xlab = "Humidity", ylab = " Temperature (%)")
cor.test(tab2$Humidity,tab2$Temperature, method = "pearson")


########## differences of temperature betweeen ages  

tab2 <- read.csv("Temp_humidity .csv")
head(tab2)

boxplot(tab2$Max..T.~ tab2$age ,xlab = "Age", ylab = " Temperature (°C) ")
fm2 <- lm (Max..T.~ age, data = tab2)
summary(fm2)

boxplot(tab2$Min.T.~ tab2$age ,xlab = "Age", ylab = " Temperature (°C)")
fm2 <- lm (Min.T.~ age, data = tab2)
summary(fm2)

boxplot(tab2$Average..T.~ tab2$age ,xlab = "Age", ylab = " Temperature (°C)")
fm2 <- lm (Average..T.~ age, data = tab2)
summary(fm2)



################differences of humidity between ages

boxplot(tab2$Max..RH.~ tab2$age ,xlab = "Age", ylab = " Relative Humidity (RH)")
fm2 <- lm (Max..RH.~ age, data = tab2)
summary(fm2)

boxplot(tab2$Min.RH.~ tab2$age ,xlab = "Age", ylab = " Relative Humidity (RH)")
fm2 <- lm (Min.RH.~ age, data = tab2)
summary(fm2)

boxplot(tab2$Average.RH.~ tab2$age ,xlab = "Age", ylab = "Relative Humidity (RH)")
fm2 <- lm (Average.RH.~ age, data = tab2)
summary(fm2)


##########################################################################

data_4<- read.csv("Temp_humidity_2.csv")
head(data_4)


#### facet wrap with Pupal_age 
ggplot(data_4,aes(x=environ ,y=Environ_values,fill=factor(age)))+
  geom_boxplot(alpha=0.3) + geom_jitter(width=0.1,alpha=0.2)+
  labs(fill = "age") + 
  geom_point(position=position_jitterdodge(),alpha=0.3) + facet_wrap(~age) +
  theme_bw(base_size = 16)

#### facet wrap with Environ measure 
ggplot(data_4,aes(x=factor(age),y=Environ_values,fill=environ))+
  geom_boxplot(alpha=0.3) + geom_jitter(width=0.1,alpha=0.2)+
  labs(fill = "environ") + 
  geom_point(position=position_jitterdodge(),alpha=0.3) + facet_wrap(~environ) +
  theme_bw(base_size = 16)
```



SUPPORTING FIGURES S2-S8





ANALYSIS OF THE IMPACT OF IRRADIATION AND TRASPORTATION(COMBINED DATA FROM 22 AND 29 DAY OLD PUPAE)



Figure S2a and S3a: Emergence rate 

```{r }


tabS2a_S3a=read.csv("Figure S2a and S3a.csv")
head(tabS2a_S3a)
Emerg_rate<- tabS2a_S3a$emerged / (tabS2a_S3a$unemerged+ tabS2a_S3a$emerged)

###Figure S2a

tiff("Figure S2a.tiff", width = 7.5, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(Emerg_rate~Treatments,data = tabS2a_S3a,
        ylab = expression(bold("Emergence rate (22 and 29 day pupae)")),
        xlab = expression(bold("Treatments")))
dev.off()


### Figure S3a

tiff("Figure S3a.tiff", width = 4, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(Emerg_rate~Pupal_age,data = tabS2a_S3a,
        ylab = expression(bold("Emergence rate (22 and 29 day pupae)")),
        xlab = expression(bold("Pupal age")))

dev.off()


###Significance
###Figure s2a
fmS2a <- glmer(cbind(emerged, unemerged) ~ Treatment +(1|Replicate), family = binomial, data = tabS2a_S3a)
summary(fmS2a )

###Figure s3a
fmS3a <- glmer(cbind(emerged, unemerged) ~ Pupal_age +(1|Replicate), family = binomial, data = tabS2a_S3a)
summary(fmS3a)

###Differences among treatments- non-parametric test
kruskal.test(Emerg_rate~Treatment, tabS2a_S3a)




####################################################################################################
```

Figure S2b and S3b: Flight propensity 


```{r }

tabS2b_S3b= read.csv("Figure S2b and S3b.csv")
head(tabS2b_S3b)
Flight_rate<- tabS2b_S3b$out / (tabS2b_S3b$in.+ tabS2b_S3b$out)
tabS2b_S3b$Pupal_age<- as.factor(tabS2b_S3b$Pupal_age)

### Figure S2b

tiff("Figure S2b.tiff", width = 7.5, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(Flight_rate ~ tabS2b_S3b$Treatments,xlab = "Treatments", ylab = "flight propensity (22 and 29 day pupae)")
boxplot(Emerg_rate~Treatments,data = tabS2b_S3b,
        ylab = expression(bold("Flight propensity (22 and 29 day)")),
        xlab = expression(bold("Treatments")))

dev.off()


###Figure S3b

tiff("Figure S3b.tiff", width = 4, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(Flight_rate ~ tabS2b_S3b$Pupal_age,xlab = "Pupal age", ylab = "flight propensity (Pupal_age)")
boxplot(Emerg_rate~Pupal_age,data = tabS2b_S3b,
        ylab = expression(bold("Flight propensity (22 and 29 day)")),
        xlab = expression(bold("Pupal age")))


dev.off()


###Significance
##Figure S2b
fmS2b <- glmer(cbind(out, in.) ~ Treatment +(1|Replicate), family = binomial, data = tabS2b_S3b)
summary(fmS2b)
###Figure S3b
fmS3b <- glmer(cbind(out, in.) ~ Pupal_age +(1|Replicate), family = binomial, data = tabS2b_S3b)
summary(fmS3b)

######## non-parametric

###Differences among treatments- non-parametric test
kruskal.test(Flight_rate~Treatment, tabS2b_S3b)



```


Figure S2c and S3c: Mating ability 


```{r }

tabS2c_S3c= read.csv("Figure S2c and S3c.csv")
head(tabS2c_S3c)
mating_prop<- tabS2c_S3c$pairs_formed  / (tabS2c_S3c$unformed_pairs+ tabS2c_S3c$pairs_formed )
tabS2c_S3c$Pupal_age<- as.factor(tabS2c_S3c$Pupal_age)

### Figure S2c

tiff("Figure S2c.tiff", width = 7.5, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(mating_prop ~ tabS2c_S3c$Treatments,xlab = "Treatments", ylab = "Mating ability")
boxplot(Emerg_rate~Treatments,data = tabS2c_S3c,
        ylab = expression(bold("Mating ability (22 and 29 day)")),
        xlab = expression(bold("Treatments")))

dev.off()


### Figure S3c

tiff("Figure S3c.tiff", width = 4, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(mating_prop ~ tabS2c_S3c$Pupal_age,xlab = "Pupal age", ylab = "Mating ability")
boxplot(Emerg_rate~Pupal_age,data = tabS2c_S3c,
        ylab = expression(bold("Mating ability (22 and 29 day)")),
        xlab = expression(bold("Pupal age")))

dev.off()


###Significance 

###Figure S2c
fmS2c <- glmer(cbind(unformed_pairs, pairs_formed)~ Treatment +(1|Replicate), family = binomial, data = tabS2c_S3c)
summary(fmS2c)

###Figure S3c
fmS3c <- glmer(cbind(unformed_pairs, pairs_formed) ~ Pupal_age +(1|Replicate), family = binomial, data = tabS2c_S3c)
summary(fmS3c)

###Differences among treatments- non-parametric test
kruskal.test(mating_prop~Treatment, tabS2c_S3c)



```


Figure S2d and S3d: Insemination rate 


```{r }

tabS2d_S3d= read.csv("Figure S2d and S3d.csv")
head(tabS2d_S3d)
Insemination_rate<- tabS2d_S3d$Inseminated / (tabS2d_S3d$Empty+ tabS2d_S3d$Inseminated)
tabS2d_S3d$Pupal_age<- as.factor(tabS2d_S3d$Pupal_age)

### Figure S2d

tiff("Figure S2d.tiff", width = 7.5, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(Emerg_rate~Treatments,data = tabS2d_S3d,
        ylab = expression(bold("Insemination rate (22 and 29 day)")),
        xlab = expression(bold("Treatments")))
dev.off()


### Figure S3d

tiff("Figure S3d.tiff", width = 4, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(Emerg_rate~Pupal_age,data = tabS2d_S3d,
        ylab = expression(bold("Insemination rate (22 and 29 day)")),
        xlab = expression(bold("Pupal age")))

dev.off()



###Significance
###Figure S2d
fmS2d <- glmer(cbind(Inseminated,Empty) ~ Treatment +(1|Replicate), family = binomial, data = tabS2d_S3d)
summary(fmS2d)

###Figure S3d
fmS3d <- glmer(cbind(Inseminated, Empty) ~ Pupal_age +(1|Replicate), family = binomial, data = tabS2d_S3d)
summary(fmS3d)

###Differences among treatments- non-parametric test

kruskal.test(Insemination_rate~Treatment, tabS2d_S3d)


```


Figure S2e and S3e:Mean spermathecal value (MSV)

```{r }

tabS2e_S3e= read.csv("Figure S2e and S3e.csv")
head(tabS2e_S3e)
tabS2e_S3e$Pupal_age<- as.factor(tabS2e_S3e$Pupal_age)

###Figure S2e

tiff("Figure S2e.tiff", width = 7.5, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(MSV~Treatments,data = tabS2e_S3e,
        ylab = expression(bold("Mean Spermathecal Value (22 and day old)")),
        xlab = expression(bold("Treatments")))
dev.off()

### Figure S3e

tiff("Figure S3e.tiff", width = 4, height =4 , units = 'in', compression = 'lzw',res = 300)
boxplot(MSV~Pupal_age,data = tabS2e_S3e,
        ylab = expression(bold("Mean Spermathecal Value (22 and 29 day old)")),
        xlab = expression(bold("Pupal age")))

dev.off()

###significance
###Figure S2e
fmS2e <- lme(MSV ~ Treatment, random=~1|Replicate,, data = tabS2e_S3e)
summary(fmS2e)
###Figure S3e
fmS3e <- lme(MSV ~ Pupal_age, random=~1|Replicate,, data = tabS2e_S3e)


summary(fmS3e)

###Differences among treatments- non-parametric test
kruskal.test(MSV~Treatment, tabS2e_S3e)


```