#Script for analysis of data in Caves et al 2018

long<-read.csv("LongFormZebraFinchData.csv") #Read the Data from CSV file

#Build The Labelling Model 
require(lme4)
#Calculate average success for each bird for each category of trial.
aggmod<-aggregate(cbind(long$pass,long$across,long$edist,long$col1,long$col2,long$contrast)
                  ,list(long$cat,long$bird.ID),FUN=mean) 
#Count number of trials for each bird for each trial
aggmod[,9]<-aggregate(cbind(long$pass,long$across,long$edist,long$col1,long$col2,long$contrast)
                      ,list(long$cat,long$bird.ID),FUN=length)[,8]
#Name columns in dataframef
names(aggmod)<-c("cat","bird.ID","pass","across","dist","col1","col2","contrast","n")
#Limit dataframe to labelling data and rows with 5+ trials per bird
aggmod<-subset(aggmod,(col1==1|col2==8)&n>4)

#Build a mixed linear model of pass rate with the binary indicator & perceptual distance as fixed effects, bird ID as a random effect
model1<-lmer(aggmod$pass~aggmod$dist+aggmod$across+(1|aggmod$bird.ID))

#Build a mixed linear model of pass rate with contrast & perceptual distance as fixed effects, bird ID as a random effect
model2<-lmer(aggmod$pass~aggmod$dist+aggmod$contrast+(1|aggmod$bird.ID))

#Compare the AIC values of the two models because the number of parameters is the same
anova(model2,model1)
#Model with the binary indicator is much better than the model with the contrast variable

#Add a random slope of the binary variable to model 1
model3<-lmer(aggmod$pass~aggmod$across+aggmod$dist+(aggmod$across|aggmod$bird.ID))

#compare models 1 & 3
anova(model1,model3)
#Model 3 much better by either likelihood ratio test or AIC

#Summary of the best model
summary(model3)

#Obtain chi squared values and p values for Table 1
#Chromatic distance (delta S):
model4<-lmer(aggmod$pass~aggmod$across+(aggmod$across|aggmod$bird.ID))
anova(model3,model4)
#Model 3 (with distance) performs better than model 4 (without distance)

#binary indicator
model5<-lmer(aggmod$pass~aggmod$dist+(aggmod$across|aggmod$bird.ID))
anova(model3,model5)
#Model 3 (with binary) performs better than model 5 (without binary)

#Visual inspection of model residuals
qqnorm(resid(model3))
qqline(resid(model3))
hist(resid(model3),10,main="",xlab="Residual Value")
abline(v=0,col="Darkgreen",lwd=3)
plot(resid(model3),ylab="Residual Value")
abline(h=0,lwd=3,col="darkgreen")

#Calculate Relative size of $Across effect as compared to delta s of 1
fixef(model3)[2]/fixef(model3)[3]

library(MuMIn)
#Calculate marginal and conditional R^2 values for LMM
r.squaredGLMM(model3)

###############Model that includes all data, not just labelling (Extended Data Table 3)#################
#Methodology is identical to that used above, but with full dataset.

#Build The Full Model
require(lmerTest)
#Calculate average success for each bird for each category of trial.
aggmod<-aggregate(cbind(long$pass,long$across,long$edist,long$col1,long$col2,long$contrast)
                  ,list(long$cat,long$bird.ID),FUN=mean) 
#Count number of trials for each bird for each trial
aggmod[,9]<-aggregate(cbind(long$pass,long$across,long$edist,long$col1,long$col2,long$contrast)
                      ,list(long$cat,long$bird.ID),FUN=length)[,8]
#Name columns in dataframe
names(aggmod)<-c("cat","bird.ID","pass","across","dist","col1","col2","contrast","n")
#Limit dataframe rows with 5+ trials per bird
aggmod<-subset(aggmod,n>4)

#Build a mixed linear model of pass rate with the binary indicator & perceptual distance as fixed effects, bird ID as a random effect
model1<-lmer(aggmod$pass~aggmod$dist+aggmod$across+(1|aggmod$bird.ID))

#Build a mixed linear model of pass rate with contrast & perceptual distance as fixed effects, bird ID as a random effect
model2<-lmer(aggmod$pass~aggmod$dist+aggmod$contrast+(1|aggmod$bird.ID))

#Compare the AIC values of the two models because the number of parameters is the same
anova(model2,model1)
#Model with the binary indicator is much better than the model with the contrast variable

#Add a random slope of the binary variable to model 1
model3<-lmer(aggmod$pass~aggmod$across+aggmod$dist+(aggmod$across|aggmod$bird.ID))

#compare models 1 & 3
anova(model1,model3)
#Model 3 much better by either likelihood ratio test or AIC

#Summary of the best model
summary(model3)

#Note that the coefficients from the model including just labelling data 
#and the model including both labelling and discrimination data
#are extremely similar

#Obtain chi squared values and p values for ED Table 3
#Chromatic distance (delta S):
model4<-lmer(aggmod$pass~aggmod$across+(aggmod$across|aggmod$bird.ID))
anova(model3,model4)
#Model 3 (with distance) performs better than model 4 (without distance)

#binary indicator
model5<-lmer(aggmod$pass~aggmod$dist+(aggmod$across|aggmod$bird.ID))
anova(model3,model5)
#Model 3 (with binary) performs better than model 5 (without binary)

#Visual inspection of model residuals
qqnorm(resid(model3))
qqline(resid(model3))
hist(resid(model3),10,main="",xlab="Residual Value")
abline(v=0,col="Darkgreen",lwd=3)
plot(resid(model3),ylab="Residual Value")
abline(h=0,lwd=3,col="darkgreen")

#Calculate Relative size of $Across effect as compared to delta s of 1
fixef(model3)[2]/fixef(model3)[3]

library(MuMIn)
#Calculate marginal and conditional R^2 values for LMM
r.squaredGLMM(model3)


#############Grayscale Discrimination Data (Table 1) #################
gray<-read.csv("GrayLongFormData.csv")

#Build a mixed model with contrast as a fixed effect and bird ID as a random effect
gmod1<-lmer(gray$pass~gray$cont+(1|gray$bird.ID))
summary(gmod1)

#Obtain chi squared and p values for Table 1, gray trials
gmod2<-lmer(gray$pass~(1|gray$bird.ID))
anova(gmod1,gmod2)

#Increased brightness increases discrimination.
#Yes, suggesting variation between birds in the magnitude of the effect of brightness on performance