# Cargar paquetes necesarios
library(car)
library(multcomp)
library(emmeans)
library(Rmisc)
library(ggplot2)

# Cargar datos
data <- read.csv("dataHR_2.csv", sep=";")

# Descripcion de los datos
summary <- summarySE(data, measurevar = "survival",
                     groupvars = c("species"))
summary

data_albo$total_eggs <- data_albo$collapsed + data_albo$total.viable + data_albo$bleached

# Agrupar por especie y RH y calcular totales y porcentajes
library(dplyr)
summary2 <- data_albo %>%
  group_by(RH, month) %>%
  summarise(
    Total_Collapsed = sum(collapsed),
    Total_Eggs = sum(total_eggs),
    Total_Viable = sum(total.viable),
    #Collapsed_Percent = 100 * sum(collapsed) / sum(total_eggs)
  )

print(summary2)


# Convertir la columna a factor 
data$RH <- as.factor(data$RH)
data$species <- as.factor(data$species)

# Combine factors if you want to analyze species × RH 
data$group <- interaction(data$species, data$RH, sep = "_")

# Modelo GLM para supervivencia
plot(data$survival)

modelo_supervivencia <- glm(
  survival ~ group,
  data = data,
  family = binomial(link = "logit")
)
summary(modelo_supervivencia)

# Calcular la desviación residual y los grados de libertad
deviance <- modelo_supervivencia$deviance
df_residual <- modelo_supervivencia$df.residual

# Test de sobredispersión
dispersion_test <- deviance / df_residual
print(dispersion_test) #>1 = sobredispersión

modelo_supervivencia <- glm(
  survival ~ species*RH,
  data = data,
  family = quasibinomial(link = "logit")
)
summary(modelo_supervivencia)

# Análisis post-hoc para supervivencia
# Obtener los valores estimados marginales (EMMs)
emmeans_obj <- emmeans(modelo_supervivencia, ~ species * RH)

# Comparaciones pareadas con Fisher's LSD
lsd_result <- cld(emmeans_obj, adjust = "none", Letters = letters)

# Convertir a dataframe y transformar del logit a proporción
emmeans_df <- as.data.frame(lsd_result)
emmeans_df$survival <- plogis(emmeans_df$emmean)  # logit to probability

ggplot(emmeans_df, aes(x = factor(RH), y = survival, fill = species)) +
  geom_bar(stat = "identity", position = position_dodge(0.8),
           colour = "black", width = 0.7) +
  geom_errorbar(aes(ymin = survival - 0.2*SE, ymax = survival + 0.2*SE),
                width = 0.1, position = position_dodge(0.8)) +
  geom_text(aes(label = .group), 
            position = position_dodge(0.8), 
            vjust = -2, size = 4, colour = "gray25") +
  ylab("Egg survival (proportion)") +
  xlab("Relative Humidity (%)") +
  labs(fill = "Species") +
  scale_fill_manual(values = c("#c1c1c1", "#505050")) +
  theme_classic() +
  theme(legend.text = element_text(face = "italic"))


# Cargar datos
data_aeg <- read.csv("RH_aegypti.csv", sep=";")

summary <- summarySE(data_aeg, measurevar = "survival",
                     groupvars = c("month", "RH"))
summary

# Convertir la columna a factor 
data_aeg$RH <- as.factor(data_aeg$RH)
data_aeg$month <- as.factor(data_aeg$month)

# Modelo GLM para supervivencia
modelo_aeg <- glm(
  survival ~ month * RH,
  data = data_aeg,
  family = quasibinomial(link = "logit")
)
summary(modelo_aeg)

# Análisis post-hoc para supervivencia
emmeans_aeg <- emmeans(modelo_aeg, ~ month * RH)
lsd_aeg <- cld(emmeans_aeg, adjust = "none", Letters = letters)

# Convertir a dataframe y transformar del logit a proporción
emmeans_aeg <- as.data.frame(lsd_aeg)
emmeans_aeg$survival <- plogis(emmeans_aeg$emmean)  # logit to probability

ggplot(emmeans_aeg, aes(x = factor(RH), y = survival, fill = month)) +
  geom_bar(stat = "identity", position = position_dodge(0.8),
           colour = "black", width = 0.7) +
  geom_errorbar(aes(ymin = survival - 0.1*SE, ymax = survival + 0.1*SE),
                width = 0.1, position = position_dodge(0.8)) +
  geom_text(aes(label = .group), 
            position = position_dodge(0.8), 
            vjust = -1.5, size = 4, colour = "gray25") +
  ylab("Egg survival (proportion)") +
  xlab("Relative Humidity (%)") +
  labs(fill = "Exposure time 
       (month)")+
  theme_minimal() +
  scale_fill_manual(values = c("#c1c1c1", "#505050", "#242424"))

# Cargar datos
data_albo <- read.csv("RH_albopictus.csv", sep=";")

summary <- summarySE(data_albo, measurevar = "survival",
                     groupvars = c("month", "RH"))
summary

# Convertir la columna a factor 
data_albo$RH <- as.factor(data_albo$RH)
data_albo$month <- as.factor(data_albo$month)

# Modelo GLM para supervivencia
modelo_albo <- glm(
  survival ~ month * RH,
  data = data_albo,
  family = quasibinomial(link = "logit")
)
summary(modelo_albo)

# Análisis post-hoc para supervivencia
emmeans_albo <- emmeans(modelo_albo, ~ month * RH)
lsd_albo <- cld(emmeans_albo, adjust = "none", Letters = letters)

# Convertir a dataframe y transformar del logit a proporción
emmeans_albo <- as.data.frame(lsd_albo)
emmeans_albo$survival <- plogis(emmeans_albo$emmean)  # logit to probability

ggplot(emmeans_albo, aes(x = factor(RH), y = survival, fill = month)) +
  geom_bar(stat = "identity", position = position_dodge(0.8),
           colour = "black", width = 0.7) +
  geom_errorbar(aes(ymin = survival - 0.1*SE, ymax = survival + 0.1*SE),
                width = 0.1, position = position_dodge(0.8)) +
  geom_text(aes(label = .group), 
            position = position_dodge(0.8), 
            vjust = -1, size = 4, colour = "gray25") +
  ylab("Egg survival (proportion)") +
  xlab("Relative Humidity (%)") +
  labs(fill = "Exposure time 
       (month)")+
  theme_minimal() +
  scale_fill_manual(values = c("#c1c1c1", "#505050", "#242424"))