#Temporal response pattern screaming cowbird

library(dplyr)
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/Rufo final")
datos<-read.xlsx("respuesta_RU_longitudinal_sinS.xlsx")

# Binomial model with interaction
m_rufo<-glmer(evento~deciles*tratamiento+(1|nido), family=binomial, datos) 

sim<-simulateResiduals(m_rufo) 
plot(sim) 
testDispersion(sim, "two.sided", plot = T) #el test da ok

summary(m_rufo)

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

#Fig 3. PLOT with predicted values and CI

library (ggeffects)
library(effects)

predichos<-ggeffect(m_rufo, terms=c("tratamiento","deciles [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)) +  # Removí el mapeo de `colour = x` aquí
  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("RN" = "#009E73", "RS" = "#E69F00"),
    labels = c("RN" = "Natural", "RS" = "Synthetic")) +
  scale_shape_manual(name = "Stimuli", values = c("RN" = 16, "RS" = 17),labels = c("RN" = "Natural", "RS" = "Synthetic")) +
  scale_fill_manual(values = c("RN" = "#009E73", "RS" = "#E69F00"),labels = c("RN" = "Natural", "RS" = "Synthetic")) +
  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),  
    axis.title.y = element_text(size = 26),  
    axis.text.x = element_text(size = 22),   
    axis.text.y = element_text(size = 22))
