####################################### (A) Packages importation ###########################################
library(ggplot2)
library(lme4)
library(lmerTest)
library(AICcmodavg)
library(plyr)
library(emmeans)
library(tidyverse)
library(car)
library(ggpubr)
library(performance) 
library(see)
library(sjPlot)
library(writexl)
library(sjtable2df)

####################################### (B) Dataset importation and preparation ############################

data_respiro_acoustic <- read.csv("_Data_respiro_acoustic.csv", header = T, sep = ";", dec = ".")
summary(data_respiro_acoustic)
data_respiro_acoustic$Hiss <- as.factor(data_respiro_acoustic$Hiss)
str(data_respiro_acoustic)


Stimu_only <- data_respiro_acoustic %>%
  dplyr::filter(Traitement == "S")

hist(data_respiro_acoustic$EWL..mg.h.)
hist(data_respiro_acoustic$VO2..ml.h.)



# Function to Create Q-Q Plots and Histograms
plot_diagnostics <- function(variable, title) {
  par(mfrow=c(2,2))
  hist(variable, main=paste("Histogram of", title), xlab=title)
  qqnorm(variable, main=paste("Q-Q Plot of", title))
  qqline(variable)
  hist(log(variable), main=paste("Histogram of log(", title, ")"), xlab=paste("log(", title, ")"))
  qqnorm(log(variable), main=paste("Q-Q Plot of log(", title, ")"))
  qqline(log(variable))
}

#Graph
plot_diagnostics(data_respiro_acoustic$EWL..mg.h., "EWL")
plot_diagnostics(data_respiro_acoustic$VO2..ml.h., "VO2")

## Log transformation for EWL and VO2



##################################### (C) Analysis  ########################################################


##### Model hissing probability #####
model_Proba_Hiss<- glmer(cbind(Hiss_Class, Hiss_Class_neg) ~ Sex*Session + Mass + (1|Individu), data = Stimu_only, family= "binomial")


## Surdipersion test
overdisp_fun <- function(model) {
  rdf <- df.residual(model)
  rp <- residuals(model,type="pearson")
  Pearson.chisq <- sum(rp^2)
  prat <- Pearson.chisq/rdf
  pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
  c(chisq=Pearson.chisq,ratio=prat,rdf=rdf,p=pval)
}


overdisp_fun(model_Proba_Hiss)
## No surdispersion (p > 0.05)



## Model output
summary(model_Proba_Hiss)
Anova(model_Proba_Hiss, type = "III")
vif(model_Proba_Hiss)



emmeans(model_Proba_Hiss, list(pairwise ~ Sex * Session), adjust = "tukey")

## Hissing probability
Stimu_only$predicted_prob<- predict(model_Proba_Hiss, type = "response")

## Summary table hissing probability according to sex and session
summary_table_sex_session <- ddply(Stimu_only, c("Sex", "Session"), summarise,
                                   N = length(predicted_prob),
                                   mean_prob = round(mean(predicted_prob), 2),
                                   sd = round(sd(predicted_prob), 2),
                                   se = round(sd/ sqrt(N), 2))

## Summary table hissing probability according to sex
summary_table_sex <- ddply(Stimu_only, c("Sex"), summarise,
                           N = length(predicted_prob),
                           mean_prob = round(mean(predicted_prob), 2),
                           sd = round(sd(predicted_prob), 2),
                           se = round(sd/ sqrt(N), 2))


summary_table_session <- ddply(Stimu_only, c("Session"), summarise,
                               N = length(predicted_prob),
                               mean_prob = round(mean(predicted_prob), 2),
                               sd = round(sd(predicted_prob), 2),
                               se = round(sd/ sqrt(N), 2))









##### Model 1  : influence treatment on VO2 #####
model_treatment_VO2 <- lmer(log(VO2..ml.h.) ~ Traitement*Sex + Session + Mass + (1|Individu), data = data_respiro_acoustic)
summary(model_treatment_VO2)
step(model_treatment_VO2)

## Check model asumptions
check_model(model_treatment_VO2)
#OK

## Model 1 simplified
model_treatment_VO2_step <- lmer(log(VO2..ml.h.) ~ Traitement*Sex + (1|Individu), data = data_respiro_acoustic)
summary(model_treatment_VO2_step)
anova(model_treatment_VO2_step, type = "marginal")

## Check model asumptions
check_model(model_treatment_VO2_step)
#OK


## Model output
tab1 <- tab_model(model_treatment_VO2_step, show.se = TRUE, show.stat = TRUE)
tab_df1 <- mtab2df(tab1, n_models = 1, output = "data.frame")
# Exportation
write_xlsx(tab_df1, "tab1.xlsx")


## Summary treament effect on VO2
summary_treatment_O2 <- ddply(data_respiro_acoustic, c("Traitement"), summarise,
                              N = length(VO2..ml.h.),
                              mean_O2 = round(mean(VO2..ml.h.), 2),
                              sd = round(sd(VO2..ml.h.), 2),
                              se = round(sd/ sqrt(N), 2))
summary_treatment_O2


## Summary treament and sex effect on VO2
summary_treatment_sex_O2 <- ddply(data_respiro_acoustic, c("Traitement", "Sex"), summarise,
                                  N = length(VO2..ml.h.),
                                  mean_O2 = round(mean(VO2..ml.h.), 2),
                                  sd = round(sd(VO2..ml.h.), 2),
                                  se = round(sd/ sqrt(N), 2))
summary_treatment_sex_O2






##### Model 2  : treatment on EWL #####
model_treatment_EWL <- lmer(log(EWL..mg.h.)~ Traitement*Sex + Session + SVL + (1|Individu), data = data_respiro_acoustic)
summary(model_treatment_EWL)
step(model_treatment_EWL)

## Check model asumptions
check_model(model_treatment_EWL)
#OK

## Model 2 Simplified 
model_treatment_EWL_step <- lmer(log(EWL..mg.h.) ~  Traitement*Sex + (1 | Individu), data = data_respiro_acoustic)
summary(model_treatment_EWL_step)
anova(model_treatment_EWL_step, type = "marginal")

## Check model asumptions
check_model(model_treatment_EWL_step)
#OK



## Model output
tab2 <- tab_model(model_treatment_EWL_step, show.se = TRUE, show.stat = TRUE)
tab_df2 <- mtab2df(tab2, n_models = 1, output = "data.frame")
# Exportation
write_xlsx(tab_df2, "tab2.xlsx")



## Summary treatment on EWL
summary_Traitement_EWL <- ddply(data_respiro_acoustic, c("Traitement"), summarise,
                                    N = length(EWL..mg.h.),
                                    mean_EWL = round(mean(EWL..mg.h.), 2),
                                    sd = round(sd(EWL..mg.h.), 2),
                                    se = round(sd/ sqrt(N), 2))

summary_Traitement_EWL


## Summary treatment and sex effect on EWL
summary_Traitement_sex_EWL <- ddply(data_respiro_acoustic, c("Traitement", "Sex"), summarise,
                         N = length(EWL..mg.h.),
                         mean_EWL = round(mean(EWL..mg.h.), 2),
                         sd = round(sd(EWL..mg.h.), 2),
                         se = round(sd/ sqrt(N), 2))

summary_Traitement_sex_EWL







##### Model 3  : influence hiss occurence on VO2 #####
model_hiss_occur_VO2 <- lmer(log(VO2..ml.h.) ~ Hiss*Sex + Session + Mass + (1|Individu), data = data_respiro_acoustic)
summary(model_hiss_occur_VO2)
step(model_hiss_occur_VO2)

## Check model asumptions
check_model(model_hiss_occur_VO2)
#OK

## Model 3 simplified
model_hiss_occur_VO2_step <- lmer(log(VO2..ml.h.) ~ Hiss*Sex + (1|Individu), data = data_respiro_acoustic)
summary(model_hiss_occur_VO2_step)
anova(model_hiss_occur_VO2_step, type = "marginal")

## Check model asumptions
check_model(model_hiss_occur_VO2_step)
#OK

emmeans(model_hiss_occur_VO2_step, list(pairwise ~ Hiss * Sex), adjust = "tukey")


## Model output
tab3 <- tab_model(model_hiss_occur_VO2_step, show.se = TRUE, show.stat = TRUE)
tab_df3 <- mtab2df(tab3, n_models = 1, output = "data.frame")
# Exportation
write_xlsx(tab_df3, "tab3.xlsx")


## Summary hiss occur on VO2
summary_hiss_occur_O2 <- ddply(data_respiro_acoustic, c("Hiss"), summarise,
                                  N = length(VO2..ml.h.),
                                  mean_O2 = round(mean(VO2..ml.h.), 2),
                                  sd = round(sd(VO2..ml.h.), 2),
                                  se = round(sd/ sqrt(N), 2))
summary_hiss_occur_O2


## Summary hiss occur on VO2
summary_hiss_occur_sex_O2 <- ddply(data_respiro_acoustic, c("Hiss","Sex"), summarise,
                               N = length(VO2..ml.h.),
                               mean_O2 = round(mean(VO2..ml.h.), 2),
                               sd = round(sd(VO2..ml.h.), 2),
                               se = round(sd/ sqrt(N), 2))
summary_hiss_occur_sex_O2









##### Model 4  : influence hiss occurence on EWL #####
model_hiss_occur_EWL <- lmer(log(EWL..mg.h.) ~ Hiss*Sex + Session + SVL + (1|Individu), data = data_respiro_acoustic)
summary(model_hiss_occur_EWL)
step(model_hiss_occur_EWL)

## Check model asumptions
check_model(model_hiss_occur_EWL)
#OK

## Model 4 simplified 
model_hiss_occur_EWL_step <- lmer(log(EWL..mg.h.) ~  Hiss + Sex + (1|Individu), data = data_respiro_acoustic)
summary(model_hiss_occur_EWL_step)
anova(model_hiss_occur_EWL_step, type = "marginal")

## Check model asumptions
check_model(model_hiss_occur_EWL_step)
#OK

emmeans(model_hiss_occur_EWL_step, list(pairwise ~ Hiss + Sex), adjust = "tukey")


## Model output
tab4 <- tab_model(model_hiss_occur_EWL_step, show.se = TRUE, show.stat = TRUE)
tab_df4 <- mtab2df(tab4, n_models = 1, output = "data.frame")
# Exportation
write_xlsx(tab_df4, "tab4.xlsx")


## Summary hiss occurence on EWL 
summary_Hiss_occurence_EWL <- ddply(data_respiro_acoustic, c("Hiss"), summarise,
                                    N = length(EWL..mg.h.),
                                    mean_EWL = round(mean(EWL..mg.h.), 2),
                                    sd = round(sd(EWL..mg.h.), 2),
                                    se = round(sd/ sqrt(N), 2))
summary_Hiss_occurence_EWL



## Summary hiss occurence and sex effect on EWL 
summary_Hiss_Sex_EWL <- ddply(data_respiro_acoustic, c("Sex"), summarise,
                              N = length(EWL..mg.h.),
                              mean_EWL = round(mean(EWL..mg.h.), 2),
                              sd = round(sd(EWL..mg.h.), 2),
                              se = round(sd/ sqrt(N), 2))
summary_Hiss_Sex_EWL


## Summary hiss occurence and session effect on EWL 
summary_Hiss_Session_EWL <- ddply(data_respiro_acoustic, c("Hiss", "Session"), summarise,
                                  N = length(EWL..mg.h.),
                                  mean_EWL = round(mean(EWL..mg.h.), 2),
                                  sd = round(sd(EWL..mg.h.), 2),
                                  se = round(sd/ sqrt(N), 2))
summary_Hiss_Session_EWL








##### Model selection ####

## EWL Model
model_treatment_EWL_step_ML <- lmer(log(EWL..mg.h.) ~  Traitement*Sex + (1 | Individu), data = data_respiro_acoustic, REML = F)
model_hiss_occur_EWL_step_ML <- lmer(log(EWL..mg.h.) ~  Hiss + Sex + (1|Individu), data = data_respiro_acoustic, REML = F)
model_null_EWL <- lmer(log(EWL..mg.h.) ~ 1 + (1 | Individu), data = data_respiro_acoustic, REML = F)

## VO2 Model
model_treatment_VO2_step_ML <- lmer(log(VO2..ml.h.) ~ Traitement*Sex + (1|Individu), data = data_respiro_acoustic, REML = F)
model_hiss_occur_VO2_step_ML <- lmer(log(VO2..ml.h.) ~ Hiss*Sex + (1|Individu), data = data_respiro_acoustic, REML = F)
model_null_VO2 <- lmer(log(VO2..ml.h.) ~ 1 + (1 | Individu), data = data_respiro_acoustic, REML = F)

## best model selection EWL  
model_EWL <- list(model_treatment_EWL_step_ML, model_hiss_occur_EWL_step_ML, model_null_EWL)

aic_table <- aictab(model_EWL)
aic_table
#### best EWL model : model_hiss_occur_EWL_step_ML


## best model selection VO2
model_VO2 <- list(model_treatment_VO2_step_ML, model_hiss_occur_VO2_step_ML, model_null_VO2)

aic_table <- aictab(model_VO2)
aic_table
#### best VO2 model : model_hiss_occur_VO2_step_ML 





##### Model 5 : Covariation between EWL and VO2 according to sex and hiss occurence #####
model_covar_hiss_occur<- lmer(log(EWL..mg.h.) ~ Hiss*log(VO2..ml.h.) + Sex + (1|Individu), data = Stimu_only)
summary(model_covar_hiss_occur)
step(model_covar_hiss_occur)


## Check model asumptions
check_model(model_covar_hiss_occur)
#OK


## Model 5 simplified
model_covar_hiss_occur_step <- lmer(log(EWL..mg.h.) ~ Hiss*log(VO2..ml.h.) + (1|Individu), data = Stimu_only)
summary(model_covar_hiss_occur_step)
anova(model_covar_hiss_occur_step, type = "marginal")


## Check model asumptions
check_model(model_covar_hiss_occur_step)
#OK



## Model output
tab5 <- tab_model(model_covar_hiss_occur_step, show.se = TRUE, show.stat = TRUE)
tab_df5 <- mtab2df(tab5, n_models = 1, output = "data.frame")
# Exportation
write_xlsx(tab_df5, "tab5.xlsx")





##### Model 6 : influence hissing duration on VO2 #####
model_hiss_duration_VO2 <- lmer(log(VO2..ml.h.) ~ scale(Hiss_sum)*Sex + Session + scale(Mass) + (1|Individu), data = Stimu_only, na.action = na.omit)
summary(model_hiss_duration_VO2)
step(model_hiss_duration_VO2)

## Check model asumptions
check_model(model_hiss_duration_VO2)
#OK

## Model 6 simplified
model_hiss_duration_VO2_step <- lmer(log(VO2..ml.h.) ~ scale(Hiss_sum) + Sex + scale(Mass) + (1 | Individu), data = Stimu_only, na.action = na.omit)
summary(model_hiss_duration_VO2_step)
anova(model_hiss_duration_VO2_step, type = "marginal")

## Check model asumptions
check_model(model_hiss_duration_VO2_step)
#OK



## Model output
tab6 <- tab_model(model_hiss_duration_VO2_step, show.se = TRUE, show.stat = TRUE)
tab_df6 <- mtab2df(tab6, n_models = 1, output = "data.frame")
# Exportation
write_xlsx(tab_df6, "tab6.xlsx")


## Beta extraction ##
model_hiss_duration_VO2_step_beta <- lmer(scale(VO2..ml.h.) ~ scale(Hiss_sum) + Sex + scale(Mass) + (1 | Individu), data = Stimu_only, na.action = na.omit)
summary(model_hiss_duration_VO2_step_beta)



## Summary sex on VO2
summary_sex_O2 <- ddply(data_respiro_acoustic, c("Sex"), summarise,
                                   N = length(VO2..ml.h.),
                                   mean_O2 = round(mean(VO2..ml.h.), 2),
                                   sd = round(sd(VO2..ml.h.), 2),
                                   se = round(sd/ sqrt(N), 2))
summary_sex_O2






##### Model 7 : influence hissing duration on EWL #####
model_hiss_duration_EWL <- lmer(log(EWL..mg.h.) ~ scale(Hiss_sum)*Sex + Session + scale(SVL) + (1|Individu), data = Stimu_only, na.action = na.omit)
summary(model_hiss_duration_EWL)
step(model_hiss_duration_EWL)

## Check model asumptions
check_model(model_hiss_duration_EWL)
#OK


## Model 7 simplified
model_hiss_duration_EWL_step <- lmer(log(EWL..mg.h.) ~ scale(Hiss_sum) + Sex + (1 | Individu),  data = Stimu_only, na.action = na.omit)
summary(model_hiss_duration_EWL_step)
anova(model_hiss_duration_EWL_step, type = "marginal")

## Check model asumptions
check_model(model_hiss_duration_EWL_step)
#OK


## Model output
tab7 <- tab_model(model_hiss_duration_EWL_step, show.se = TRUE, show.stat = TRUE)
tab_df7 <- mtab2df(tab7, n_models = 1, output = "data.frame")
# Exportation
write_xlsx(tab_df7, "tab7.xlsx")

## Summary sex on EWL
summary_sex_EWL <- ddply(data_respiro_acoustic, c("Sex"), summarise,
                         N = length(EWL..mg.h.),
                         mean_EWL = round(mean(EWL..mg.h.), 2),
                         sd = round(sd(EWL..mg.h.), 2),
                         se = round(sd/ sqrt(N), 2))
summary_sex_EWL


## Beta extraction ##
model_hiss_duration_EWL_step_beta <- lmer(scale(EWL..mg.h.) ~ scale(Hiss_sum) + Sex + (1 | Individu),  data = Stimu_only, na.action = na.omit)
summary(model_hiss_duration_EWL_step_beta)



#################################### (D) Figures ###########################################################

summary_table_sex_session$Sex <- factor(summary_table_sex_session$Sex, levels = c("M", "F"))

##### Figure 1 hissing probability according to session and sex ##### 
Figure_1 <- ggplot(summary_table_sex_session, aes(x = Session, y = mean_prob, color = Sex, fill = Sex)) +
  geom_errorbar(aes(ymin = mean_prob - se, ymax = mean_prob + se),
                width = 0.3, position = position_dodge(width = 0.9), size = 1) +
  geom_point(stat = "summary", position = position_dodge(width = 0.9), size = 5) +
  labs(x = "Session", y = "Hissing probability") +
  scale_shape_manual(values = c("Session_1" = 16, "Session_2" = 17)) +
  scale_color_manual(values = c("M" = "#00BFC4", "F" = "#F8766D")) +
  scale_fill_manual(values = c("M" = "#00BFC4", "F" = "#F8766D")) +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 16, hjust = 0.5),
        axis.text.y = element_text(size = 16),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20, vjust = 1.5),
        legend.text = element_text(size = 18),
        legend.title = element_text(size = 18),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(),
        axis.ticks = element_line(size = 1, color = "black"),
        legend.position = c(0.8, 0.8),
        legend.key.width = unit(1, "cm"), legend.key.height = unit(1, "cm"),
        strip.text = element_blank()) +
  scale_y_continuous(limits = c(0, 1), expand = c(0, 0)) + 
  facet_wrap(~Sex)

Figure_1 <- Figure_1 + theme(plot.margin = margin(t = 6, r = 6, b = 6, l = 6, unit = "mm"))
ggsave("Figure_1.svg",Figure_1,width=250,height=200,units=c("mm"),dpi=900,bg="transparent",limitsize = FALSE)


##### Figure 3a Inlfuence hiss occurence on VO2 ##### 
Figure_3a <- ggplot(data_respiro_acoustic, aes(x = Hiss, y = log(VO2..ml.h.), color = Sex)) +
  geom_boxplot(fill = "transparent", size = 0.7, outlier.shape = NA) +
  stat_summary(geom = "point", shape = 17, size = 5, 
               mapping = aes(color = factor(Sex)),
               position = position_dodge(width = 0.75)) +
  labs(x = "Hiss Occurrence", y = expression(log ~ VO2 ~ (ml.~ h^{-1}))) +
  scale_color_manual(values = c("M" = "#00BFC4", "F" = "#F8766D")) +
  scale_x_discrete(labels = c("No", "Yes")) +
  scale_y_continuous(limits = c(1, 5)) +
  geom_jitter(position = position_jitterdodge(jitter.width = 0.15, dodge.width = 0.75), alpha = 0.5) +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 16, hjust = 0.5, color = "black"),
        axis.text.y = element_text(size = 16),
        axis.title.x = element_text(size = 20, vjust = 0.1),
        legend.text = element_text(size = 14),
        legend.title = element_text(size = 18, hjust = 0.3),
        strip.text = element_text(size = 12),
        axis.title.y = element_text(size = 20, vjust = 1.5),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(),
        axis.ticks = element_line(size = 1, color = "black"), 
        legend.position = c(0.09, 0.8)) +
  guides(fill = guide_legend(override.aes = list(size = 5))) +
  theme(legend.key.width = unit(2.5, "cm"), legend.key.height = unit(2.5, "cm"))



##### Figure 3b Inlfuence hiss occurence on EWL ##### 
Figure_3b <- ggplot(data_respiro_acoustic, aes(x = Hiss, y = log(EWL..mg.h.), color = Sex)) +
  geom_boxplot(fill = "transparent", size = 0.7, outlier.shape = NA) +
  stat_summary(geom = "point", shape = 17, size = 5, 
               mapping = aes(color = factor(Sex)),
               position = position_dodge(width = 0.75)) +
  labs(x = "Hiss Occurrence", y = expression(log ~ EWL ~ (mg.~ h^{-1}))) +
  scale_color_manual(values = c("M" = "#00BFC4", "F" = "#F8766D")) +
  scale_x_discrete(labels = c("No", "Yes")) +
  scale_y_continuous(limits = c(3, 5), breaks = seq(2, 5, by = 1)) +
  geom_jitter(position = position_jitterdodge(jitter.width = 0.15, dodge.width = 0.75), alpha = 0.5) +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 16, hjust = 0.5, color = "black"),
        axis.text.y = element_text(size = 16),
        axis.title.x = element_text(size = 20, vjust = 0.1),
        legend.text = element_text(size = 14),
        legend.title = element_text(size = 18, hjust = 0.3),
        strip.text = element_text(size = 12),
        axis.title.y = element_text(size = 20, vjust = 1.5),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(),
        axis.ticks = element_line(size = 1, color = "black"), 
        legend.position = c(0.09, 0.8)) +
  guides(fill = guide_legend(override.aes = list(size = 5))) +
  theme(legend.key.width = unit(2.5, "cm"), legend.key.height = unit(2.5, "cm"))




figure_3a <- Figure_3a + theme(plot.margin = margin(t = 6, r = 5, b = 10, l = 7, unit = "mm"))
figure_3b <- Figure_3b + theme(plot.margin = margin(t = 10, r = 5, b = 6, l = 5, unit = "mm"))


## Figure fusion
Figure_3 <- ggarrange(figure_3a, figure_3b, ncol = 1, nrow = 2)

ggsave("Figure_3.svg",Figure_3,width=250,height=400,units=c("mm"),dpi=900,limitsize = FALSE)



#### Figure 4 Covariation between EWL and VO2 according to sex ##### 
Figure_4 <- ggplot(Stimu_only, aes(x = log(VO2..ml.h.), y = log(EWL..mg.h.), color = Hiss)) +
  geom_point(shape = 1, size = 4, stroke = 2) +
  geom_smooth(method = "lm", se = T, formula = y ~ x, fill = "#d0d0d0") +
  labs(x = expression(log ~ VO2 ~ (ml.~ h^{-1})), y = expression(log ~ EWL ~ (mg.~ h^{-1}))) +
  scale_color_manual(values = c("0" = "#0CB702", "1" = "#E68613"), breaks = c("1", "0"),labels = c("Yes", "No")) + 
  scale_y_continuous(limits = c(3, 5), breaks = seq(2, 5, by = 1)) +
  guides(linetype = guide_legend(title = "Hissing Occurence"), color = guide_legend(title = "Hiss Occur")) +
  theme(legend.key.size = unit(1.5, "lines")) + theme_minimal() +
  theme(axis.text.x = element_text(size = 16, hjust = 0.5),
        axis.text.y = element_text(size = 16),
        axis.title.x = element_text(size = 20, vjust = 0.1),
        legend.text = element_text(size = 14), 
        legend.title = element_text(size = 16),
        strip.text = element_text(size = 12),
        axis.title.y = element_text(size = 20, vjust = 1.5), 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(), 
        axis.ticks = element_line(size = 1, color = "black"), 
        legend.position = c(0.8, 0.2)) +
  theme(legend.key.width = unit(1.5, "cm"), legend.key.height = unit(1.5, "cm"))

Figure_4 <- Figure_4 + theme(plot.margin = margin(t = 4, r = 4, b = 4, l = 4, unit = "mm"))


ggsave("Figure_4.svg",Figure_4,width=250,height=200,units=c("mm"),dpi=900,bg="transparent",limitsize = FALSE)


##### Figure 5a Inlfuence hiss duration on VO2 ##### 
Stimu_only$Sex <- factor(Stimu_only$Sex, levels = c("M", "F"))

Figure_5a <- ggplot(Stimu_only, aes(x = Hiss_sum, y = log(VO2..ml.h.), color = Sex)) +
  geom_point(shape = 1, size = 4, stroke = 2) +
  geom_smooth(method = "lm", se = T, formula = y ~ x, fill = "#d0d0d0") +
  labs(x = "Hiss duration (s)", y = expression(log ~ VO2 ~ (ml.~ h^{-1}))) +
  scale_color_manual(values = c("M" = "#00BFC4", "F" = "#F8766D")) + 
  guides(linetype = guide_legend(title = "Sex"), shape = guide_legend(title = "Sex"), color = guide_legend(title = "Sex")) +
  theme(legend.key.size = unit(1.5, "lines")) + theme_minimal() +
  theme(axis.text.x = element_text(size = 16, hjust = 0.5),
        axis.text.y = element_text(size = 16),
        axis.title.x = element_text(size = 20, vjust = 0.1),
        legend.text = element_text(size = 14), 
        legend.title = element_text(size = 16),
        strip.text = element_text(size = 12),
        axis.title.y = element_text(size = 20, vjust = 1.5), 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(), 
        axis.ticks = element_line(size = 1, color = "black"), 
        legend.position = c(0.8, 0.2)) +
  guides(fill = guide_legend(override.aes = list(size = 5))) +
  theme(legend.key.width = unit(1.5, "cm"), legend.key.height = unit(1.5, "cm"))



##### Figure 5b Influence hiss duration on EWL #####
Stimu_only$Sex <- factor(Stimu_only$Sex, levels = c("M", "F"))

Figure_5b <- ggplot(Stimu_only, aes(x = Hiss_sum, y = log(EWL..mg.h.), color = Sex)) +
  geom_point(shape = 1, size = 4, stroke = 2) +
  geom_smooth(method = "lm", se = T, formula = y ~ x, fill = "#d0d0d0") +
  labs(x = "Hiss duration (s)", y = expression(log ~ EWL ~ (mg.~ h^{-1}))) +
  scale_color_manual(values = c("M" = "#00BFC4", "F" = "#F8766D")) + 
  scale_y_continuous(limits = c(3, 5), breaks = seq(2, 5, by = 1)) +
  guides(linetype = guide_legend(title = "Sex"), shape = guide_legend(title = "Sex"), color = guide_legend(title = "Sex")) +
  theme(legend.key.size = unit(1.5, "lines")) + theme_minimal() +
  theme(axis.text.x = element_text(size = 16, hjust = 0.5),
        axis.text.y = element_text(size = 16),
        axis.title.x = element_text(size = 20, vjust = 0.1),
        legend.text = element_text(size = 14), 
        legend.title = element_text(size = 16),
        strip.text = element_text(size = 12),
        axis.title.y = element_text(size = 20, vjust = 1.5), 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(), 
        axis.ticks = element_line(size = 1, color = "black"), 
        legend.position = c(0.8, 0.2)) +
  guides(fill = guide_legend(override.aes = list(size = 5))) +
  theme(legend.key.width = unit(1.5, "cm"), legend.key.height = unit(1.5, "cm"))



figure_5a <- Figure_5a + theme(plot.margin = margin(t = 5, r = 6, b = 10, l = 8, unit = "mm"))
figure_5b <- Figure_5b + theme(plot.margin = margin(t = 10, r = 6, b = 2, l = 5, unit = "mm"))



## Figure fusion
Figure_5 <- ggarrange(figure_5a, figure_5b, ncol = 1, nrow = 2)
ggsave("Figure_5.svg",Figure_5,width=250,height=400,units=c("mm"),dpi=900,limitsize = FALSE)