#10 September 2024
#Luke D. Fannin
#GLMM code for Rosien et al., 2024, "Food-washing monkeys recognize the law of diminishing returns", elife
#Script Revised from 12 April 2023, in response to reviewer feedback

#General Note: Analyses run in RStudio using R version 4.4.1 (2024-06-14) -- "Race for Your Life"

#Notes on datasets and analyses:

#Revised Script composes two GLMM analyses
#1 ---> Macaque brushing data: dependent variable (brushing time, seconds)
#       Independent variables (fixed) - ordinal rank, sex, grit treatment (high, medium, low), and the interaction
#       Individual ID will be a random variable
#       Sexes combined for the analysis, all under the same ordinal rank hierarchy 

#2 ---> Macaque washing data: dependent variable (washing time, seconds)
#       Independent variables (fixed) - ordinal rank, sex, grit treatment (high, medium, low), and the interaction
#       Individual ID will be a random variable
#       Sexes combined for the analysis, all under the same ordinal rank hierarchy 

#Both datasets were modeled with a Poisson error distribution.
#Lots of biologically-meaningful zeros, so zero-inflation also applied.

#Load requisite packages;
library(ggplot2) #for plots 
library(glmmTMB)#generalized linear model package for running the zero-inflated models
library(car)#for generating fixed effects model tests 
library(DHARMa) # creating residual diagnostic plots for assessing model fits

#####################################################################################
#Analysis #1 - brushing data (n = 575 actions, updated 22 Aug 2024)
#Note 1- Half seconds rounded up to the nearest second in this dataset.
#Note 2 - Both sexes combined into a single ordinal rank hierarchy.
#Note 3 - The working directory below will need to be changed to match where the user has placed the relevant .csv file
#Note 4 - relevant .csv file in zenodo is: “2024.Sept.10.brushing.GLMM.csv”

#Load dataset:
brush <- read.csv(“***insert working directory file path here***”)

#check data frame
head(brush)

#Making a visualization
brush_plot <- ggplot(brush, aes(x=o.rank, y=duration.sec, color=grit.treatment,shape=grit.treatment))+
  geom_smooth(method=lm,fullrange=FALSE,se=FALSE)+
  geom_point()+
  scale_y_continuous(limits=c(0,13),breaks=seq(0,13,3))+
  scale_x_continuous(limits=c(1,18),breaks=seq(1,18,3))+
  theme(legend.position="left")+
  scale_color_manual(values=pal)

#Display plot:
brush_plot

#Generate the model below
#Generating a zero-inflated poisson distribution
fit_zi_poisson_brush <- glmmTMB(duration.sec~o.rank+ 
                                  grit.treatment + 
                                  sex + 
                                  o.rank*grit.treatment +
                                  (1|ind),
                                  data=brush,ziformula =~1,
                                  family="poisson")

#Summary:
summary(fit_zi_poisson_brush)

#Confidence intervals of all the fixed effects:
confint(fit_zi_poisson_brush)

#Analysis of Deviance Table for Parameter Estimates
if (requireNamespace("car") && getRversion() >= "3.6.0") {
  Anova(fit_zi_poisson_brush) ## default type II
  Anova(fit_zi_poisson_brush,type="II")
}

#Model Diagnostics, DHARMa Simulated Standardized Residuals and Plots
#Generate simulated standardized residuals:
simulationOutput <- simulateResiduals(fittedModel = fit_zi_poisson_brush, plot = F)

#Generate QQ and quantile residual vs. predicted plots:
plot(simulationOutput)

#Generate just a QQ plot:
plotQQunif(simulationOutput,testUniformity = FALSE,testDispersion = FALSE,testOutliers = FALSE)

#Generate a dispersion statistic for the residuals:
testDispersion(simulationOutput)

#Generate a residual vs. predicted plot, without quantiles:
plotResiduals(simulationOutput,quantreg=FALSE)

######################################################################################
#Analysis #2 - washing data (n = 362 actions, updated 22 Aug 2024)
#Note 1- Half seconds rounded up to the nearest second in this dataset.
#Note 2 - Both sexes combined into a single ordinal rank hierarchy.
#Note 3 - The working directory below will need to be changed to match where the user has placed the relevant .csv file
#Note 4 - The relevant .csv file in Zenodo is “2024.Sept.10.washing.GLMM.csv”

#Load dataset:
wash <- read.csv(“***insert working directory file path here***”)

#Check data:
head(wash)

#Create aa data vizualization:
pal<- c("#3B9AB2","#E5BD15","#F21A00")
wash_plot <- ggplot(wash, aes(x=o.rank, y=duration.sec, color=grit.treatment,shape=grit.treatment))+
  geom_smooth(method=lm,fullrange=FALSE,se=FALSE)+
  geom_point()+
  scale_y_continuous(limits=c(0,20),breaks=seq(0,20,4))+
  scale_x_continuous(limits=c(1,24),breaks=seq(1,24,3))+
  theme(legend.position="left")+
  scale_color_manual(values=pal)

#display plot
wash_plot

#Fitting zero-inflation poisson model with a quadratic term for ordinal rank: 
fit_zipoisson_wash.quad <- glmmTMB(duration.sec~o.rank+ 
                                        I(o.rank^2)+
                                        grit.treatment + 
                                        sex+
                                        I(o.rank^2)*grit.treatment+
                                        (1|ind),data=wash,ziformula =~1,
                                        family="compois") 

#Model Summary
summary(fit_zipoisson_wash.quad)

#Confidence Intervals
confint(fit_zipoisson_wash.quad)

#ANOVA table for poisson washing
if (requireNamespace("car") && getRversion() >= "3.6.0") {
  Anova(fit_zipoisson_wash.quad) ## default type II
  Anova(fit_zipoisson_wash.quad,type="II")
}

#Diagnostic Plots via DHARMa:
#Generate standardized simulated residuals:
simulationOutput.washing.quad <- simulateResiduals(fittedModel = fit_zipoisson_wash.quad, plot = F)

#Plot QQ plot plus residual vs. predicted plot with quantiles:
plot(simulationOutput.washing.quad)

#Generate Dispersion Statistics:
testDispersion(simulationOutput.washing.quad)

#Generate just a QQ plot:
plotQQunif(simulationOutput.washing.quad,testUniformity = FALSE,testDispersion = FALSE,testOutliers = FALSE)

#Generate a residual vs. predicted plot, without quantiles:
plotResiduals(simulationOutput.washing.quad,quantreg=FALSE)

#####################################################################################
#Alternative washing model
#Fitting the washing model, but without the quadratic:
fit_zipoisson_wash_non_quadratic <- glmmTMB(duration.sec~o.rank+ 
                                grit.treatment + 
                                sex+
                                o.rank*grit.treatment+
                                (1|ind),data=wash,ziformula =~1,
                                family="compois") 

#Full Model Summary:
summary(fit_zipoisson_wash_non_quadratic)

#Model Conidence Intervals:
confint(fit_zipoisson_wash_non_quadratic)

#Generate the Analysis of Deviance Table:
if (requireNamespace("car") && getRversion() >= "3.6.0") {
  Anova(fit_zipoisson_wash_non_quadratic) ## default type II
  Anova(fit_zipoisson_wash_non_quadratic,type="II")
}

#Model Diagnostic Plots Generated with DHARMa:

#Generate model residuals
simulationOutput.washing <- simulateResiduals(fittedModel = fit_zipoisson_wash_non_quadratic, plot = F) 

#Produce dual QQ plot and residual vs. predicted plot with quantiles:
plot(simulationOutput.washing)

#Generate just a QQ plot:
plotQQunif(simulationOutput.washing,testUniformity = FALSE,testDispersion = FALSE,testOutliers = FALSE)

#Generate residual vs. predicted plot, no quantiles:
plotResiduals(simulationOutput.washing,quantreg=FALSE)

#Generate a dispersion statistic for the model:
testDispersion(simulationOutput.washing)