##Manip analysis R File from
##EM Caves, LE Schweikert, PA Green, MN Zipple, C Taboada, S Peters, S Nowicki, and S Johnsen. “Variation in carotenoid-containing retinal oil droplets correlates with variation in perception of carotenoid coloration.” Behavioral Ecology and Sociobiology. 
##Author(s) for correspondence: Stephen Nowicki (snowicki@duke.edu) and Sönke Johnsen (sjohnsen@duke.edu),  Duke University


#Set your working directory to source location
library(lme4)
library(lmerTest)
library(afex)

manip<-read.csv("upload_file.csv")
lambda<-read.csv("lambda2.csv")
lambda$lambdacent<-lambda$lambda-mean(lambda$lambda)


#manip contains data from Caves et al, 2018 and the pass rate from the manipulation experiment "pass"#

modNoRand<-lm(pass~across1+across2+across3+across4+across5+
                across6+across7+treatment,data=subset(manip,bird.ID!="RED108"&
                                                        col2-col1<=3&is.na(lambda)==F))
modNull<-lmer(pass~across1+across2+across3+across4+across5+
                across6+across7+treatment+(1|bird.ID),data=manip)
afexNull<-mixed(pass~across1+across2+across3+across4+across5+
        across6+across7+relevel(treatment,ref=1)+(1|bird.ID),data=manip
      , method = "KR", control = lmerControl(optCtrl = list(maxfun = 1e6)), expand_re = TRUE)
summary(modNull)
summary(afexNull)

modNullDist<-lmer(pass~edist+treatment+(1|bird.ID),data=manip)
afexNullDist<-mixed(pass~edist+relevel(treatment,ref=1)+(1|bird.ID),data=manip
                    , method = "KR", control = lmerControl(optCtrl = list(maxfun = 1e6)), expand_re = TRUE)
summary(modNullDist)
summary(afexNullDist)

modNullContrast<-lmer(pass~treatment+contrast+(1|bird.ID),data=manip)
afexNullContrast<-mixed(pass~relevel(treatment,ref=1)+contrast+(1|bird.ID),data=manip
                    , method = "KR", control = lmerControl(optCtrl = list(maxfun = 1e6)), expand_re = TRUE)
summary(modNullContrast)
summary(afexNullContrast)
anova(modNull,modNullContrast,modNullDist)

modLamb<-lmer(pass~lambcent+across1+across2+across3+across4+across5+
                 across6+across7+treatment+(1|bird.ID),data=manip)
afexLamb<-mixed(pass~lambcent+across1+across2+across3+across4+across5+
                  across6+across7+relevel(treatment,ref=1)+(1|bird.ID),data=manip,
                method = "KR", control = lmerControl(optCtrl = list(maxfun = 1e6)), expand_re = TRUE)
summary(modLamb)
summary(afexLamb)

anova(modLamb,modNull)
anova(afexNull,afexNullDist,afexLamb)

modLamb5<-lmer(pass~lambcent+across1+across2+across3+across4+across5+
                across6+across7+treatment+(across5|bird.ID),data=manip)
afexLamb5<-mixed(pass~lambcent+across1+across2+across3+across4+across5+
        across6+across7+relevel(treatment,ref=1)+(across5|bird.ID),data=manip,
      method = "KR", control = lmerControl(optCtrl = list(maxfun = 1e6)), expand_re = TRUE)
summary(modLamb5)
summary(afexLamb5)

modLambInt<-lmer(pass~lambcent*across5+across1+across2+across3+across4+
                    across6+across7+treatment+(1|bird.ID),data=manip)
afexLambInt<-mixed(pass~lambcent+across1+across2+across3+across4+lambcent*across5+
                           across6+across7+treatment+(1|bird.ID),data=manip,
                   method = "KR", control = lmerControl(optCtrl = list(maxfun = 1e6)), expand_re = TRUE)
summary(modLambInt)
summary(afexLambInt)

anova(afexNull,afexLamb,afexLambInt)
anova(modNull,modLamb,modLambInt)

modLambInt5<-lmer(pass~treatment+across1+across2+across3+across4+
                   across6+across7+lambcent*across5+(across5|bird.ID),data=manip)
afexLambInt5<-mixed(pass~lambcent+across1+across2+across3+across4+lambcent*across5+
                     across6+across7+treatment+(across5|bird.ID),data=manip,
                   method = "KR", control = lmerControl(optCtrl = list(maxfun = 1e6)), expand_re = TRUE)
summary(modLambInt5)
summary(afexLambInt5)

qqnorm(resid(modLambInt5))
qqline(resid(modLambInt5))

anova(modNull,modLamb,modLambInt)

modLambIntAll<-lmer(pass~across5*lambcent+across2*lambcent+across1*lambcent+lambcent*across7+across3*lambcent+
                      across6*lambcent+across4*lambcent+treatment+(1|bird.ID),data=subset(manip,bird.ID!="RED108"&      col2-col1<=3&is.na(lambda)==F),REML=F,control=lmerControl(optCtrl=list(maxfun=20000) ))
afexLambIntAll<-mixed(pass~across5*lambcent+across2*lambcent+across1*lambcent+lambcent*across7+across3*lambcent+
                     across6*lambcent+across4*lambcent+treatment+(1|bird.ID),data=manip,
                   method = "KR", control = lmerControl(optCtrl = list(maxfun = 1e6)), expand_re = TRUE)

summary(modLambIntAll)
summary(afexLambIntAll)

modLambIntAllSlopes<-lmer(pass~across1*lambcent+across2*lambcent+across3*lambcent+
                           across4*lambcent+across5*lambcent+across6*lambcent+lambcent*across7+treatment+(across2+across4+across5+across6+across7|bird.ID),data=manip)
afexLambIntAllSlopes<-mixed (pass~across1*lambcent+across2*lambcent+across3*lambcent+across4*lambcent+across5*lambcent+across6*lambcent+across7*lambcent+treatment+(across2+across4+across5+across6+across7|bird.ID), 
                       data = manip, method = "KR", control = lmerControl(optCtrl = list(maxfun = 1e6)), expand_re = TRUE)
summary(modLambIntAllSlopes)
summary(afexLambIntAllSlopes)

####################################
#Figure that Shows contribution of Across 5 versus mean lambda mid
across5coef<-lmer(pass~across1+across2+across3+across4+across5+
                    across6+across7+treatment+(across5|bird.ID),data=manip)
across5_random_effects<-coef(across5coef)
across5_random_effects
lambda[,c(2,4)] #Note that we have added the random slope estimate of the
#effect of crossing the 5-6 boundary to the "lambda" file, which contains
#our average measures of lambda for each bird

plot(lambda$across5~lambda$lambda,ylab="Across 5 Coeficient",xlab="Lambda",pch=19,cex=2)
lm_lambda_across5<-lm(lambda$across5~lambda$lambda)
summary(lm_lambda_across5)