#10 September 2024
#Luke D Fannin
#Code for figures for Rosien et al., 2024 “Food-washing monkeys recognize the law of diminishing returns”
#Figures made in RStudio in R version 4.4.1
######################################################################################
#Figure 1 - Sand Particle Size Distributions and Compositions 
#Note - Working directory will have to change

#Loading packages for Figure 1
library(ggplot2)
library(ggsci) #For color pallete styles
library(wesanderson)
library(gridExtra) #Aligning multiple plots on the same grid

#Loading sand dataset (file name: “2024.Sept.10.sand.particle.dataset.csv” on Zenodo)
sand <- read.csv(“***Insert working directory file path here***”)

#Using the wes anderson color gradient:
pal <- wes_palette("Zissou1",8,type="continuous")
pal2 <- c("#3B9AB2","#78B7C5","#B1C177",
          "#EBCC2A","#E5BD15","#E1AF00","#E96400",
          "#F21A00")###Tan Aquatic color pallate

######################################################################################
#Panel A - Histogram of particle size on the x axis:
#Density plot with histogram version A
xdensity <- ggplot(data=sand,aes(x=mean.diameter,fill=classification.narrow)) +
  geom_histogram(aes(y=..count..),binwidth = 30) +
  scale_fill_manual(values=pal2)+
  geom_density(aes(y = ..density.. * nrow(sand) * 30), position = "identity",
               colour = "black", fill = "grey", alpha = 0.0)+
  theme(legend.position="none")+
  geom_vline(xintercept = 25,linetype="dashed")

#Export Figure
ggsave(filename = "fig1a_particle_distribution_w_scale.pdf",width =14, height = 14,units = "cm",dpi = 300,useDingbats = FALSE)

######################################################################################
#Panel #2 - scatter plot, size on the x axis, angularity on the y axis, same dataset as panel A
scatterPlot <- ggplot(sand,aes(x=mean.diameter,y=circularity,color=classification.narrow,group=1))+
  geom_point()+
  scale_color_manual(values=pal2)+
  geom_smooth(span=1.0,se=FALSE)+
  theme(legend.position="none")+
  geom_vline(xintercept = 25,linetype="dashed")
#If line is desired:
stat_smooth(method="lm",formula=y~poly(x,2),fullrange = FALSE)

#Arranging panels A and B on the same  plot
grid.arrange(xdensity,scatterPlot, 
             ncol=1, nrow=2, widths=c(4), heights=c(3, 4))

g <- arrangeGrob(xdensity,scatterPlot, 
ncol=1, nrow=2, widths=c(4), heights=c(3, 4))

#Save figure file:
ggsave(filename = "fig1a_particle_angularity_size.pdf",g,width =12, height = 14,units = "cm",dpi = 300,useDingbats = FALSE)

######################################################################################
#Code for Figure 2, in text (dot plots with standard errors, three category ordinal rank score)
#Ordinal ranks divided into three evenly spaced categories: “High” (1-8); “Middle”: (9-16); and “Low”: (17-24)
#Note - working directories will differ depending on where the author places the requisite .csv file

library(ggplot2)

#####
#Panel B——>Brushing dot plot, three-rank category
#Note - dataset is titled “2024.Sept.10.intake.dot.chart.figure2.csv”

#Load dataset:
brushing.dot.plot<-read.csv(“***insert working directory here***”,stringsAsFactors = TRUE)

#Check dataset
head(brushing.dot.plot)

#Plot code
brushing.plot.dot<- ggplot(brushing.dot.plot, aes(x=rank, y=ave.duration, group=treatment, color=treatment,shape=treatment)) + 
  geom_pointrange(aes(ymin=ave.duration-se.duration, ymax=ave.duration+se.duration))+
  scale_shape_manual(values = c(17, 15, 16))+
  scale_y_continuous(limits=c(0,9),breaks=seq(0,9,3))+
  theme(legend.position="none")

#With a different scale on the y axis:
brushing.plot.dot<- ggplot(brushing.dot.plot, aes(x=rank, y=ave.duration, group=treatment, color=treatment,shape=treatment)) + 
  geom_pointrange(aes(ymin=ave.duration-se.duration, ymax=ave.duration+se.duration))+
  scale_shape_manual(values = c(17, 15, 16))+
  theme(legend.position="none")

#View plot:
brushing.plot.dot

#Save plot:
ggsave(filename = "brushing.dot.plot.rank.categories.fig.2.pdf",width =6, height = 7,units = "cm",dpi = 300,useDingbats = FALSE)

#####
#Panel C——>Washing dot plot, three-rank category
#Note: dataset is titled “2024.Sept.10.washing.dot.chart.figure2.csv”

#Load dataset:
washing.dot.plot<-read.csv("***insert working directory here***",stringsAsFactors = TRUE)
head(washing.dot.plot)

#Plot code
washing.plot.dot<- ggplot(washing.dot.plot, aes(x=rank, y=ave.duration, group=treatment, color=treatment,shape=treatment)) + 
  geom_pointrange(aes(ymin=ave.duration-se.duration, ymax=ave.duration+se.duration))+
  scale_shape_manual(values = c(17, 15, 16))+
  scale_y_continuous(limits=c(0,9),breaks=seq(0,9,3))+
  theme(legend.position="none")

#With an alternative scale on the y-axis:
washing.plot.dot<- ggplot(washing.dot.plot, aes(x=rank, y=ave.duration, group=treatment, color=treatment,shape=treatment)) + 
  geom_pointrange(aes(ymin=ave.duration-se.duration, ymax=ave.duration+se.duration))+
  scale_shape_manual(values = c(17, 15, 16))+
  theme(legend.position="none")

#Display plot:
washing.plot.dot

#Save plot:
ggsave(filename = "washing.dot.plot.rank.categories.fig.2.pdf",width =6, height = 7,units = "cm",dpi = 300,useDingbats = FALSE)

#####
#Panel D——>Intake rate dot plot, three-rank category
#Note: dataset is titled: “2024.Sept.10.intake.dot.chart.figure2.csv”

#Load dataset:
intake.rate.dot.plot<-read.csv("***insert working directory here***",stringsAsFactors = TRUE)

#Display intake rate data
head(intake.rate.dot.plot)

intake.plot.dot<- ggplot(intake.rate.dot.plot, aes(x=rank, y=intake)) + 
  geom_pointrange(aes(ymin=intake-se, ymax=intake+se))+
  theme(legend.position="none")

#Display plot:
intake.plot.dot

#Export Plot
ggsave(filename = "intake.dot.plot.rank.categories.fig.2.pdf",width =6, height = 7,units = "cm",dpi = 300,useDingbats = FALSE)

#####################################################################################
#Supplemental Figure 2 
#Note - working directories used below will not be the same as the ones used by new users
library(ggplot2)

#######
#Panel A ——> brushing, all individuals (listed in the same, single ordinal rank hierarchy)
#Note - dataset titled: “2024.Sept.10.brushing.bar.chart.figureS2.csv”

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

#Rank category - convert to a factor
brushing.figure$rank=as.factor(brushing.figure$rank)

#Bar Plot:
ggplot(data=brushing.figure, aes(x=rank, y=ave.duration, fill=treatment)) +
  geom_bar(stat="identity")+
  theme(legend.position="none")+
  scale_y_continuous(limits=c(0,14),breaks=seq(0,14,2))

#Export Plot:
ggsave(filename = "brushing.v.optimality.all.individuals.o.ranks.bar.chart.pdf",width =10, height = 4,units = "cm",dpi = 300,useDingbats = FALSE)

#######
#PanelB——> washing, all individuals (listed in the same, single ordinal rank hierarchy)
#Note - dataset titled: “2024.Sept.10.washing.bar.chart.figureS2.csv”

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

#Rank category - convert to a factor
washing.figure$rank=as.factor(washing.figure$rank)

#Bar Plot:
ggplot(data=washing.figure, aes(x=rank, y=ave.duration, fill=treatment)) +
  geom_bar(stat="identity")+
  theme(legend.position="none")+
  scale_y_continuous(limits=c(0,14),breaks=seq(0,14,2))

#Export Plot:
ggsave(filename = "washing.v.optimality.all.individuals.o.ranks.bar.chart.pdf",width =10, height = 4,units = "cm",dpi = 300,useDingbats = FALSE)

#######
#Panel C——> intake rate (in KJoules/min), all individuals (listed in the same, single ordinal rank hierarchy)
#Note - dataset titled: “2024.Sept.10.intake.bar.chart.figureS2”

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

#Rank category - convert to a factor
intake.figure$rank=as.factor(intake.figure$rank)

#Bar Plot:
ggplot(data=intake.figure, aes(x=rank, y=intake)) +
  geom_bar(stat="identity")+
  theme(legend.position="none")

#Export plot
ggsave(filename = "intake.v.optimality.all.individuals.o.ranks.bar.chart.pdf",width =10, height = 4,units = "cm",dpi = 300,useDingbats = FALSE)