#Playback cowbird 
library(ggplot2) 
library(car)     
library(emmeans) 
library(dplyr)    
library(psych)  
library(performance)
library(lme4)
library(survival)
library(coxme)
library(glmmTMB)
library(DHARMa)


#FOR BAYWING SET
rm(list=ls())

setwd("E:/CAMI/Campaña 2022-2023/MANUSCRITO 2024/Base de datos y script para BL/base datos usadas")
datos <- read.delim("baywing_orden.txt")
datos$treatment<- as.factor(datos$treatment)
datos$nest<- as.factor(datos$nest)


#FOR COWBIRD SET
# rm(list=ls())
# 
# setwd("E:/CAMI/Campaña 2022-2023/MANUSCRITO 2024/Base de datos y script para BL/base datos usadas")
# datos <- read.delim("cowbird_orden_1.txt")
# datos$treatment<- as.factor(datos$treatment)
# datos$nest<- as.factor(datos$nest)

#0: Non-mimetic, 1: Natural, 2: Synthetic
library(forcats)
datos$treatment<-fct_relevel(datos$treatment, "S")
levels(datos$treatment)



#LATENCY

#Coxme
lat<-coxme(formula=Surv(latency, event)~treatment+order+(1|nest), data=datos) #full model
cox.zph(lat)

# LR test treatment significance
lat.sin.trat<-coxme(formula=Surv(latency, event)~order+(1|nest), data=datos) #without treatment
cox.zph(lat.sin.trat)

anova(lat, lat.sin.trat, test="Chisq")

#LR test order significance
lat.sin.orden<-coxme(formula=Surv(latency, event)~treatment+(1|nest), data=datos)#without order
cox.zph(lat.sin.orden)

anova(lat, lat.sin.orden, test="Chisq")

#Comparisons
emmeans(lat.sin.orden, pairwise~treatment, type="response")



#RECRUITMENT

#Binomial
recruitment<-glmmTMB(recruitment ~ treatment+order+ (1|nest), family = binomial, weights = adults ,  data=datos) #full model

plot(simulateResiduals(recruitment)) 
testResiduals(recruitment)

# LR test treatment significance
recruit.sin.trat <-glmmTMB(recruitment ~ order+ (1|nest), family = binomial, weights = adults,  data=datos) #without treatment

plot(simulateResiduals(recruit.sin.trat))

anova(recruitment, recruit.sin.trat) 

# LR test order significance
recruit.sin.orden <-glmmTMB(recruitment ~ treatment+ (1|nest), family = binomial, weights = adults,  data=datos) #without order

plot(simulateResiduals(recruit.sin.orden))

anova(recruitment, recruit.sin.orden)

#Comparisons
emmeans(recruitment, pairwise ~ treatment, type = "response")



#DURATION

#Zero Inflated (binomial)
duration.zi<- glmmTMB(duration~treatment+order+(1|nest), zi=~treatment, data = datos, family = beta_family()) #full model

plot(simulateResiduals(duration.zi))

# LR test treatment significance (Zero inflated)
duration.null <- glmmTMB(duration ~ treatment + order + (1 | nest), zi = ~1, data = datos, family = beta_family()) #without treatment in zero inflated model

anova(duration.null, duration.zi)

#Comparisons for ZI
emmeans(duration.zi, pairwise ~ treatment, component = "zi", type = "response")

# LR test treatment significance (conditional)
duration.zi.no.treat<- glmmTMB(duration~order+(1|nest), zi=~treatment, data = datos, family = beta_family()) #without treament

anova(duration.zi.no.treat, duration.zi)

# LR test order significance (conditional)
duration.zi.no.order<- glmmTMB(duration~treatment+(1|nest), zi=~treatment, data = datos, family = beta_family()) #without order

anova(duration.zi.no.order, duration.zi)
