#Temporal response pattern baywing

library(openxlsx)
library(lme4)
library(DHARMa)
library(ggplot2)
library(emmeans)
library(multcomp)

rm(list=ls())

setwd("E:/CAMI/Campaña 2022-2023/MANUSCRITO 2024/bases de datos grafico gabo/archivos R/Musico final")
datos <- read.xlsx("respuesta_ab_longitudinal_sinS.xlsx")
names(datos)[4]<-paste("periodo")

# Binomial model with interaction
m2<-glmer(respuesta~periodo*tratamiento+(periodo|nido), family=binomial, datos)

sim<-simulateResiduals(m2) 
plot(sim)
summary(m2)

#Post HOC Comparisons between treatments by period
m2 %>%
  emmeans::ref_grid(at = list(periodo = c(1, 2, 3, 4, 5, 6, 7, 8, 9))) %>%
  emmeans::contrast(method = "pairwise", by="periodo", type="response") 

#Fig 3. PLOT with predicted values and CI

library (ggeffects)
library(effects)

predichos<-ggeffect(m_ab, terms=c("tratamiento","periodo [1:10]"), ci_level=0.95, type="sim")
class(predichos$group)
predichos$deciles = as.numeric(predichos$group) 

ggplot(predichos, aes(x = deciles, y = predicted, color = x, shape = x)) + 
  geom_point(size = 4) + 
  geom_line(aes(y = predicted, x = deciles)) + 
  geom_ribbon(aes(x = deciles, ymin = conf.low, ymax = conf.high, fill = x), alpha = 0.3, show.legend = FALSE) + 
  scale_x_continuous(name = "Time (deciles)", breaks = seq(1, 10, 1), limits = c(1, 10)) +   scale_y_continuous(name = "Proportion of nests") + 
  scale_color_manual(
    name = "Stimuli", 
    values = c("MN" = "#009E73", "MS" = "#E69F00"),
    labels = c("MN" = "Natural", "MS" = "Synthetic")
  ) +
  scale_shape_manual(
    values = c("MN" = 16, "MS" = 17),
    labels = c("MN" = "Natural", "MS" = "Synthetic")
  ) +
  scale_fill_manual(values = c("MN" = "#009E73", "MS" = "#E69F00")) + 
  theme_classic()+
  theme(
    legend.position = "none", 
    legend.text = element_text(size = 22), 
    legend.title = element_text(size = 20, face = "bold"), 
    axis.title.x = element_text(size = 26),  # Tamaño del título del eje x
    axis.title.y = element_text(size = 26),  # Tamaño del título del eje y
    axis.text.x = element_text(size = 22),   # Tamaño de las etiquetas del eje x
    axis.text.y = element_text(size = 22)    # Tamaño de las etiquetas del eje y
  )

