#setting the working directory setwd("") #GLMM_ attempts - 2013-14 - both - nestling binomKPattempt<-read.csv("1i - attempts - 2013-14 - both - nestling.csv") names(binomKPattempt) head(binomKPattempt) #setting the factors from the data sheet binomKPattempt$colony<-factor(binomKPattempt$Colony) binomKPattempt$chick<-factor(binomKPattempt$NEWID_Rep) binomKPattempt$year<-factor(binomKPattempt$Year) binomKPattempt$preytype<-factor(binomKPattempt$PREYTYPE) binomKPattempt$preysize<-factor(binomKPattempt$PREYSIZE) binomKPattempt$klepattempt<-factor(binomKPattempt$KPFACTOR) table(binomKPattempt$preytype,binomKPattempt$colony) table(binomKPattempt$preysize,binomKPattempt$colony) table(binomKPattempt$colony,binomKPattempt$KP.ATTEMPT) table(binomKPattempt$year,binomKPattempt$chick) library(lme4) library(multcomp) library(LMERConvenienceFunctions) library(lmtest) library(glmmADMB) library(lsmeans) #writing the global model to investigate overall differences in KP attempts in the two colonies binomKPattemptPS<-subset(binomKPattempt,!preysize=="NA") ## To eliminate NAs in preysize from the csv file binKPattempt1a<-glmer(KP.ATTEMPT~colony+preysize+(1|year/chick), data=binomKPattemptPS, family=binomial) ##DATA=attempt"PS" summary(binKPattempt1a) binKPattempt1b<-glmer(KP.ATTEMPT~colony+(1|year/chick), data=binomKPattemptPS, family=binomial) summary(binKPattempt1b) # model simplification lrtest(binKPattempt1a, binKPattempt1b) ##no significant change### with prey size removed binKPattempt2<-glmer(KP.ATTEMPT~colony+(1|year/chick), data=binomKPattempt, family=binomial) ##BEST summary(binKPattempt2) # binKPattempt3<-glmer(KP.ATTEMPT~colony+TideFactor+(1|year/chick), data=binomKPattempt, family=binomial) # summary(binKPattempt3) #Generate our predicted means modpredict<-glmer(KP.ATTEMPT~colony+(1|year/chick), data=binomKPattempt, family=binomial) predict<- expand.grid(colony=c("Jetty", "Quarry"), KP.ATTEMPT=0) mm<- model.matrix(terms(modpredict),predict) predict$KP.ATTEMPT<-mm %*% fixef(modpredict) #multiply matrix with fixed effect to get predicted response (transformed value predict$se1<-sqrt(diag(mm %*% vcov(modpredict) %*% t(mm))) #fancy SE of predicted response ###Back transformed values (glmer only) predict$KP.ATTEMPT.bt<-plogis(predict$KP.ATTEMPT) # back transform predicted response and SE predict$se1.uppbt <-plogis(predict$KP.ATTEMPT + predict$se1) predict$se1.lowbt <-plogis(predict$KP.ATTEMPT - predict$se1) predict ##PLOT stolenprobi <- lsmeans (modpredict, "colony") stolenprobi sp_means<-matrix(c(plogis( -0.2281474),plogis( -2.5453534)),dimnames=list(c("Jetty","Quarry"))) sp_upp<-matrix(c(plogis(-0.2281474+0.1818350),plogis(-2.393696+0.2538055)),dimnames=list(c("Jetty","Quarry"))) sp_low<-matrix(c(plogis(-0.2281474-0.1818350),plogis(-2.393696-0.2538055)),dimnames=list(c("Jetty","Quarry"))) par(mar=c(5,4,1,2),mgp=c(3,1,0)) bp1i<-barplot(sp_means,beside=TRUE,ylim=c(0,0.5),space=0.1,names.arg=c("Jetty","Quarry"),xlab="Colony",ylab="Likelihood of Kleptoparasitism attempt") arrows(bp1i,sp_upp,bp1i,sp_low,code=3,angle=90) GLMM_ attempts - 2013-15 - Jetty - mobnest JETTYKPattempt<-read.csv("2 - attempts - 2013-15 - Jetty - mobnest.csv") names(JETTYKPattempt) head(JETTYKPattempt) #setting the factors from the data sheet JETTYKPattempt$colony<-factor(JETTYKPattempt$Colony) JETTYKPattempt$chick<-factor(JETTYKPattempt$NEWID_Rep) JETTYKPattempt$year<-factor(JETTYKPattempt$Year) JETTYKPattempt$preytype<-factor(JETTYKPattempt$PREYTYPE) JETTYKPattempt$preysize<-factor(JETTYKPattempt$PREYSIZE) JETTYKPattempt$age<-factor(JETTYKPattempt$BINARYAGE) table(JETTYKPattempt$preytype,JETTYKPattempt$colony) table(JETTYKPattempt$preysize,JETTYKPattempt$colony) table(JETTYKPattempt$colony,JETTYKPattempt$KP.ATTEMPT) table(JETTYKPattempt$year,JETTYKPattempt$chick) JETTYKPattempt$windscale<-scale(JETTYKPattempt$Wind, center = TRUE, scale = TRUE) library(lme4) library(multcomp) library(LMERConvenienceFunctions) library(lmtest) library(glmmADMB) #writing the global model to investigate H gull differences in KP attempt with CHICK AGE, tide, and food size. JETTYhgullKPattempt1<-glmer(BINOMIAL.HGULL.ATTEMPT~age+TideFactor+preysize+(1|year/chick), data=JETTYKPattempt, family=binomial) summary(JETTYhgullKPattempt1) #using the backfit function for glmers bfFixefLMER_t.fnc(JETTYhgullKPattempt1, item = FALSE, method = c("llrt"), threshold = NULL, t.threshold = NULL, alphaitem = NULL, prune.ranefs = TRUE, set.REML.FALSE = TRUE, keep.single.factors=FALSE, reset.REML.TRUE = TRUE, log.file = NULL) JETTYhgullKPattemptsize<-glmer(BINOMIAL.HGULL.ATTEMPT~preysize+(1|year/chick), data=JETTYKPattempt, family=binomial) summary(JETTYhgullKPattemptsize) JETTYhgullKPattempttideage1<-glmer(BINOMIAL.HGULL.ATTEMPT~age+TideFactor+(1|year/chick), data=JETTYKPattempt, family=binomial) summary(JETTYhgullKPattempttideage1) bfFixefLMER_t.fnc(JETTYhgullKPattempttideage1, item = FALSE, method = c("llrt"), threshold = NULL, t.threshold = NULL, alphaitem = NULL, prune.ranefs = TRUE, set.REML.FALSE = TRUE, keep.single.factors=FALSE, reset.REML.TRUE = TRUE, log.file = NULL) #Tern kp attempts JETTYTERNKPattempt1<-glmer(BINOMIAL.TERN.ATTEMPT~age+TideFactor+preysize+(1|year/chick), data=JETTYKPattempt, family=binomial) summary(JETTYTERNKPattempt1) #using the backfit function for glmers bfFixefLMER_t.fnc(JETTYTERNKPattempt1, item = FALSE, method = c("llrt"), threshold = NULL, t.threshold = NULL, alphaitem = NULL, prune.ranefs = TRUE, set.REML.FALSE = TRUE, keep.single.factors=FALSE, reset.REML.TRUE = TRUE, log.file = NULL) JETTYTERNKPattemptsize<-glmer(BINOMIAL.TERN.ATTEMPT~preysize+(1|year/chick), data=JETTYKPattempt, family=binomial) summary(JETTYTERNKPattemptsize) JETTYTERNKPattempttideage1<-glmer(BINOMIAL.TERN.ATTEMPT~age+TideFactor+(1|year/chick), data=JETTYKPattempt, family=binomial) summary(JETTYTERNKPattempttideage1) bfFixefLMER_t.fnc(JETTYTERNKPattempttideage1, item = FALSE, method = c("llrt"), threshold = NULL, t.threshold = NULL, alphaitem = NULL, prune.ranefs = TRUE, set.REML.FALSE = TRUE, keep.single.factors=FALSE, reset.REML.TRUE = TRUE, log.file = NULL) #Items stolen, H gull first JETTYhgullKPstolen1<-glmer(BINOMIAL.HGULL.STOLEN~age+TideFactor+preysize+(1|year/chick), data=JETTYKPattempt, family=binomial) summary(JETTYhgullKPstolen1) bfFixefLMER_t.fnc(JETTYhgullKPstolen1, item = FALSE, method = c("llrt"), threshold = NULL, t.threshold = NULL, alphaitem = NULL, prune.ranefs = TRUE, set.REML.FALSE = TRUE, keep.single.factors=FALSE, reset.REML.TRUE = TRUE, log.file = NULL) JETTYhgullKPstolensize<-glmer(BINOMIAL.HGULL.STOLEN~preysize+(1|year/chick), data=JETTYKPattempt, family=binomial) summary(JETTYhgullKPstolensize) JETTYhgullKPstolenagetide<-glmer(BINOMIAL.HGULL.STOLEN~age+TideFactor+(1|year/chick), data=JETTYKPattempt, family=binomial) summary(JETTYhgullKPstolenagetide) bfFixefLMER_t.fnc(JETTYhgullKPstolenagetide, item = FALSE, method = c("llrt"), threshold = NULL, t.threshold = NULL, alphaitem = NULL, prune.ranefs = TRUE, set.REML.FALSE = TRUE, keep.single.factors=FALSE, reset.REML.TRUE = TRUE, log.file = NULL) JETTYhgullKPstolenage<-glmer(BINOMIAL.HGULL.STOLEN~age+(1|year/chick), data=JETTYKPattempt, family=binomial) summary(JETTYhgullKPstolenage) JETTYternKPstolen1<-glmer(BINOMIAL.TERN.STOLEN~age+TideFactor+preysize+(1|year/chick), data=JETTYKPattempt, family=binomial) summary(JETTYternKPstolen1) bfFixefLMER_t.fnc(JETTYternKPstolen1, item = FALSE, method = c("llrt"), threshold = NULL, t.threshold = NULL, alphaitem = NULL, prune.ranefs = TRUE, set.REML.FALSE = TRUE, keep.single.factors=FALSE, reset.REML.TRUE = TRUE, log.file = NULL) JETTYternKPstolensize<-glmer(BINOMIAL.TERN.STOLEN~preysize+(1|year/chick), data=JETTYKPattempt, family=binomial) summary(JETTYternKPstolensize) JETTYternKPstolenagetide<-glmer(BINOMIAL.TERN.STOLEN~age+TideFactor+(1|year/chick), data=JETTYKPattempt, family=binomial) summary(JETTYternKPstolenagetide) bfFixefLMER_t.fnc(JETTYternKPstolenagetide, item = FALSE, method = c("llrt"), threshold = NULL, t.threshold = NULL, alphaitem = NULL, prune.ranefs = TRUE, set.REML.FALSE = TRUE, keep.single.factors=FALSE, reset.REML.TRUE = TRUE, log.file = NULL) JETTYternKPstolenage<-glmer(BINOMIAL.TERN.STOLEN~age+(1|year/chick), data=JETTYKPattempt, family=binomial) summary(JETTYternKPstolenage) #########################################################################################GLMM_ success - 2013-15 - Jetty - mobnest######################### JETTYKPsuccess<-read.csv("3 - success - 2013-15 - Jetty - mobnest.csv") names(JETTYKPsuccess) head(JETTYKPsuccess) #setting the factors from the data sheet JETTYKPsuccess$colony<-factor(JETTYKPsuccess$Colony) JETTYKPsuccess$chick<-factor(JETTYKPsuccess$NEWID_Rep) JETTYKPsuccess$year<-factor(JETTYKPsuccess$Year) JETTYKPsuccess$preytype<-factor(JETTYKPsuccess$PREYTYPE) JETTYKPsuccess$preysize<-factor(JETTYKPsuccess$PREYSIZE) JETTYKPsuccess$age<-factor(JETTYKPsuccess$BINARYAGE) JETTYKPsuccess$species<-factor(JETTYKPsuccess$SPECIESKP) JETTYKPsuccess$tide<-factor(JETTYKPsuccess$TideFactor) JETTYKPsuccess$species<-relevel(JETTYKPsuccess$species, ref="HARTLAUB") table(JETTYKPsuccess$preytype,JETTYKPsuccess$colony) JETTYKPsuccess$windscale<-scale(JETTYKPsuccess$Wind, center = TRUE, scale = TRUE) library(lme4) library(multcomp) library(LMERConvenienceFunctions) library(lmtest) library(glmmADMB) Overallkpsuccessall<-glmer(KP.SUCCESS~species+age+tide+preysize+(1|year/chick), data=JETTYKPsuccess, family=binomial) summary(Overallkpsuccessall) #using the backfit function for glmers bfFixefLMER_t.fnc(Overallkpsuccessall, item = FALSE, method = c("llrt"), threshold = NULL, t.threshold = NULL, alphaitem = NULL, prune.ranefs = TRUE, set.REML.FALSE = TRUE, keep.single.factors=FALSE, reset.REML.TRUE = TRUE, log.file = NULL) OverallkpsuccessSIZE<-glmer(KP.SUCCESS~age+preysize+(1|year/chick), data=JETTYKPsuccess, family=binomial) summary(OverallkpsuccessSIZE) Overallkpsuccess<-glmer(KP.SUCCESS~species+age+tide+(1|year/chick), data=JETTYKPsuccess, family=binomial) summary(Overallkpsuccess) #using the backfit function for glmers bfFixefLMER_t.fnc(Overallkpsuccess, item = FALSE, method = c("llrt"), threshold = NULL, t.threshold = NULL, alphaitem = NULL, prune.ranefs = TRUE, set.REML.FALSE = TRUE, keep.single.factors=FALSE, reset.REML.TRUE = TRUE, log.file = NULL) modpredict4<-glmer(BINOMIAL.TERN.STOLEN~colony+(1|year/chick), data=binomKPattempt, family=binomial) predict4<- expand.grid(colony=c("Jetty", "Quarry"), BINOMIAL.TERN.STOLEN=0) mm4<- model.matrix(terms(modpredict4),predict4) predict4$BINOMIAL.TERN.STOLEN<-mm4 %*% fixef(modpredict4) #multiply matrix with fixed effect to get predicted response (transformed value predict4$se1<-sqrt(diag(mm4 %*% vcov(modpredict4) %*% t(mm4))) #fancy SE of predicted response ###Back transformed values (glmer only) predict4$BINOMIAL.TERN.STOLEN.bt<-plogis(predict4$BINOMIAL.TERN.STOLEN) # back transform predicted response and SE predict4$se1.uppbt <-plogis(predict4$BINOMIAL.TERN.STOLEN + predict4$se1) predict4$se1.lowbt <-plogis(predict4$BINOMIAL.TERN.STOLEN - predict4$se1) predict4 ##################################################################################################GLMM_ parent swallow self ############################# BinomSELF<-read.csv("5xiv - parent swallow self NEW.csv") names(BinomSELF) head(BinomSELF) # #setting the factors from the data sheet with y= swallow sefl (binomial) BinomSELF$chick<-factor(BinomSELF$NEWID_Rep) BinomSELF$year<-factor(BinomSELF$Year) BinomSELF$KPAtt<-factor(BinomSELF$KP.ATTEMPT) BinomSELF$preysize<-factor(BinomSELF$PREYSIZE) BinomSELF$age<-factor(BinomSELF$BINARYAGE) BinomSELF$species<-factor(BinomSELF$SPECIESKP) BinomSELF$tide<-factor(BinomSELF$TideFactor) par(mfrow=c(1,2)) library(lme4) library(multcomp) library(lmtest) library(lsmeans) BSelfJQ <- BinomSELF[BinomSELF$Year %in% c("_2013","_2014"),] ############################################## OverallBINSELF1<-glmer(BINSELF~age+KPAtt+preysize+(1|year/chick), data=BSelfJQ, family=binomial) summary(OverallBINSELF1) #so first remove preysize as previously described, then you can proceed OverallBINSELF2<-glmer(BINSELF~age+KPAtt+(1|year/chick), data=BSelfJQ, family=binomial) OverallBINSELF3<-glmer(BINSELF~KPAtt+(1|year/chick), data=BSelfJQ, family=binomial) summary(OverallBINSELF3) modpredictXIV<-glmer(BINSELF~KPAtt+(1|year/chick), data=BSelfJQ, family=binomial) predictXIV<- expand.grid(KPAtt=c("0", "1"), BINSELF=0) mmXIV<- model.matrix(terms(modpredictXIV),predictXIV) predictXIV$BINSELF<-mmXIV %*% fixef(modpredictXIV) #multiply matrix with fixed effect to get predicted response (transformed value predictXIV$se1<-sqrt(diag(mmXIV %*% vcov(modpredictXIV) %*% t(mmXIV))) #fancy SE of predicted response # ###Back transformed values (glmer only) predictXIV$BINSELF.bt <-plogis(predictXIV$BINSELF) # back transform predicted response and SE predictXIV$se1.uppbt <-plogis(predictXIV$BINSELF + predictXIV$se1) predictXIV$se1.lowbt <-plogis(predictXIV$BINSELF - predictXIV$se1) predictXIV ##PLOT KPATT Swallowself <- lsmeans (modpredictXIV, "KPAtt") Swallowself sp_means<-matrix(c(plogis(-5.615537),plogis(-4.979677)),dimnames=list(c("0", "1"))) sp_upp<-matrix(c(plogis(-5.615537+1.252325),plogis(-4.979677+1.155123)),dimnames=list(c("0", "1"))) sp_low<-matrix(c(plogis(-5.615537-1.252325),plogis(-4.979677-1.155123)),dimnames=list(c("0", "1"))) par(mar=c(5,4,1,2),mgp=c(3,1,0)) bpss<-barplot(sp_means,beside=TRUE,ylim=c(0,0.06),space=0.1,names.arg=c("NOATT", "ATT"),xlab="Both colonies",ylab="Likelihood of parent swallow self") arrows(bpss,sp_upp,bpss,sp_low,code=3,angle=90) BSelfJALL <- BinomSELF[BinomSELF$Colony %in% c("Jetty"),] OverallALLSELF1<-glmer(BINSELF~age+KPAtt+preysize+(1|year/chick), data=BSelfJALL, family=binomial) summary(OverallALLSELF1) OverallALLSELF2<-glmer(BINSELF~age+KPAtt+(1|year/chick), data=BSelfJALL, family=binomial) ## NO Converge??########## OverallALLSELF3<-glmer(BINSELF~KPAtt+(1|year/chick), data=BSelfJALL, family=binomial) ## Same as model 3 above? summary(OverallALLSELF3) modpredictXIVb<-glmer(BINSELF~KPAtt+(1|year/chick), data=BSelfJALL, family=binomial) predictXIVb<- expand.grid(KPAtt=c("0", "1"), BINSELF=0) mmXIVb<- model.matrix(terms(modpredictXIVb),predictXIVb) predictXIVb$BINSELF<-mmXIVb %*% fixef(modpredictXIVb) #multiply matrix with fixed effect to get predicted response (transformed value predictXIVb$se1<-sqrt(diag(mmXIVb %*% vcov(modpredictXIVb) %*% t(mmXIVb))) #fancy SE of predicted response # ###Back transformed values (glmer only) # Your code here in the next line was wrong for the predicted mean backtransformed, you deleted the '.bt <' somehow predictXIVb$BINSELF.bt <-plogis(predictXIVb$BINSELF) # back transform predicted response and SE predictXIVb$se1.uppbt <-plogis(predictXIVb$BINSELF + predictXIVb$se1) predictXIVb$se1.lowbt <-plogis(predictXIVb$BINSELF - predictXIVb$se1) predictXIVb ########################################################################### #################GLMM_ factors affecting feedign attempts####################### JettyFeedAtt<-read.csv("5i - factors affecting feedign attempts.csv") names(JettyFeedAtt) head(JettyFeedAtt) # # #setting the factors from the data sheet with y= swallow sefl (binomial) # JettyFeedAtt$chick<-factor(JettyFeedAtt$NEWID_Rep) JettyFeedAtt$year<-factor(JettyFeedAtt$Year) JettyFeedAtt$KPAtt<-factor(JettyFeedAtt$KP.ATTEMPT) JettyFeedAtt$preysize<-factor(JettyFeedAtt$PREYSIZE) JettyFeedAtt$age<-factor(JettyFeedAtt$BINARYAGE) library(lme4) library(multcomp) library(LMERConvenienceFunctions) library(lmtest) library(glmmADMB) library(lsmeans) poissonJETTYfeedattempt1<-glmer(Feeding.attempts~age+KPAtt+preysize+(1|year/chick), data=JettyFeedAtt, family=poisson) ##NO CONVERGENCE### summary(poissonJETTYfeedattempt1) poissonJETTYfeedattempt3<-glmer(Feeding.attempts~KPAtt+preysize+(1|year/chick), data=JettyFeedAtt, family=poisson) ##NO CONVERGENCE### summary(poissonJETTYfeedattempt3) poissonJETTYfeedattempt4<-glmer(Feeding.attempts~preysize+KPAtt+(1|year/chick), data=JettyFeedAtt, family=poisson) ##NO CONVERGENCE### summary(poissonJETTYfeedattempt4) poissonJETTYfeedattempt2<-glmer(Feeding.attempts~age+KPAtt+(1|year/chick), data=JettyFeedAtt, family=poisson) ##BEST## summary(poissonJETTYfeedattempt2) poissonJETTYfeedattempt5<-glmer(Feeding.attempts~preysize+(1|year/chick), data=JettyFeedAtt, family=poisson) ##maybe separately prey size only??## summary(poissonJETTYfeedattempt5) par(mfrow=c(1,2)) ## AGE ## feedattemptsJETTYAge<- lsmeans (poissonJETTYfeedattempt2, "age") feedattemptsJETTYAge sp_means<-matrix(c(exp(0.4978229),exp(0.6572812)),dimnames=list(c("Mobile","Nestling"))) sp_upp<-matrix(c(exp(0.4978229+0.05335659),exp(0.6572812+0.04631172)),dimnames=list(c("Mobile","Nestling"))) sp_low<-matrix(c(exp(0.4978229-0.05335659),exp(0.6572812-0.04631172)),dimnames=list(c("Mobile","Nestling"))) par(mar=c(5,4,1,2),mgp=c(3,1,0)) bp5xv<-barplot(sp_means,beside=TRUE,ylim=c(0,3),space=0.1,names.arg=c("Mobile","Nestling"),xlab="Age",ylab="Number of Feeding attempt") arrows(bp5xv,sp_upp,bp5xv,sp_low,code=3,angle=90) ## KP ATTEMPTS ## feedattemptsJETTYKP<- lsmeans (poissonJETTYfeedattempt2, "KPAtt") feedattemptsJETTYKP sp_means<-matrix(c(exp(0.2435837),exp(0.9115204)),dimnames=list(c("0","1"))) sp_upp<-matrix(c(exp(0.2435837+0.05098647),exp(0.9115204+0.04594070)),dimnames=list(c("0","1"))) sp_low<-matrix(c(exp(0.2435837-0.05098647),exp(0.9115204-0.04594070)),dimnames=list(c("0","1"))) par(mar=c(5,4,1,2),mgp=c(3,1,0)) bp5xvb<-barplot(sp_means,beside=TRUE,ylim=c(0,3),space=0.1,names.arg=c("NO ATTEMPT","ATTEMPT"),ylab="") arrows(bp5xvb,sp_upp,bp5xvb,sp_low,code=3,angle=90) par(mfrow=c(1,1)) ## PREY SIZE### NO BIG EFFCT!! feedattemptsJETTYKPps<- lsmeans (poissonJETTYfeedattempt5, "preysize") feedattemptsJETTYKPps sp_means<-matrix(c(exp(0.6493000),exp(0.6120286)),dimnames=list(c("LARGE","SMALL"))) sp_upp<-matrix(c(exp(0.6493000+0.10004666),exp(0.6120286+0.08207987)),dimnames=list(c("LARGE","SMALL"))) sp_low<-matrix(c(exp(0.6493000-0.10004666),exp(0.6120286-0.08207987)),dimnames=list(c("LARGE","SMALL"))) par(mar=c(5,4,1,2),mgp=c(3,1,0)) bp5xvC<-barplot(sp_means,beside=TRUE,ylim=c(0,3),space=0.1,names.arg=c("LARGE","SMALL"),xlab="Prey size",ylab="") arrows(bp5xvC,sp_upp,bp5xvC,sp_low,code=3,angle=90) #############################################################################################GLMM_ factors affecting feedign attempts by Species############## #setting an abbreviation for the data file JFeedAttSp<-read.csv("5xvi- factors affecting feedign attempts by Species.csv") # #setting the factors from the data sheet JFeedAttSp$chick<-factor(JFeedAttSp$NEWID_Rep) JFeedAttSp$year<-factor(JFeedAttSp$Year) JFeedAttSp$KPAtt<-factor(JFeedAttSp$KP.ATTEMPT) JFeedAttSp$preysize<-factor(JFeedAttSp$PREYSIZE) JFeedAttSp$age<-factor(JFeedAttSp$BINARYAGE) JFeedAttSp$species<-factor(JFeedAttSp$SPECIES) ### Relevel species JFeedAttSp$species<-relevel(JFeedAttSp$SPECIES, ref="GULL") library(lme4) library(multcomp) library(lmtest) library(lsmeans) poissJFeedAttSP1<-glmer(Feeding.attempts~KPAtt+species+preysize+(1|year/chick/new.event), data=JFeedAttSp, family=poisson) summary(poissJFeedAttSP1) poissJFeedAttSP2<-glmer(Feeding.attempts~species+(1|year/chick/new.event), data=JFeedAttSp, family=poisson) ##Best?? summary(poissJFeedAttSP2) lrtest(poissJFeedAttSP1, poissJFeedAttSP2) poissJFeedAttSP3<-glmer(Feeding.attempts~KPAtt+(1|year/chick/eventID), data=JFeedAttSp, family=poisson) summary(poissJFeedAttSP3) lrtest(poissJFeedAttSP1, poissJFeedAttSP3) par(mfrow=c(1,3)) ##By species JfeedattemptsSP<- lsmeans (poissJFeedAttSP1, "species") JfeedattemptsSP sp_means<-matrix(c(exp(0.4878659),exp(0.4878655)),dimnames=list(c("Gull","Tern"))) sp_upp<-matrix(c(exp(0.4878659+0.07127151),exp(0.4878655+0.07127151)),dimnames=list(c("Gull","Tern"))) sp_low<-matrix(c(exp(0.4878659-0.07127151),exp(0.4878655-0.07127151)),dimnames=list(c("Gull","Tern"))) par(mar=c(5,4,1,2),mgp=c(3,1,0)) bp5xvi<-barplot(sp_means,beside=TRUE,ylim=c(0,3),space=0.1,names.arg=c("Gull","Tern"),xlab="",ylab="Number of feeding passes") arrows(bp5xvi,sp_upp,bp5xvi,sp_low,code=3,angle=90) ##By KPAtt JfeedattemptsKP<- lsmeans (poissJFeedAttSP1, "KPAtt") JfeedattemptsKP sp_means<-matrix(c(exp(0.1647132),exp(0.8110182)),dimnames=list(c("0","1"))) sp_upp<-matrix(c(exp(0.1647132+0.07703198),exp(0.8110182+0.06963378)),dimnames=list(c("0","1"))) sp_low<-matrix(c(exp(0.1647132-0.07703198),exp(0.8110182-0.06963378)),dimnames=list(c("0","1"))) par(mar=c(5,4,1,2),mgp=c(3,1,0)) bp5xvikp<-barplot(sp_means,beside=TRUE,ylim=c(0,3),space=0.1,names.arg=c("No KP Attempt","KP Attempt"),ylab="") arrows(bp5xvikp,sp_upp,bp5xvikp,sp_low,code=3,angle=90) ##By Prey size JfeedattemptsPS<- lsmeans (poissJFeedAttSP1, "preysize") JfeedattemptsPS sp_means<-matrix(c(exp(0.4652090),exp(0.5105223)),dimnames=list(c("LARGE","SMALL"))) sp_upp<-matrix(c(exp(0.4652090+0.08261830),exp(0.5105223+0.06658481)),dimnames=list(c("LARGE","SMALL"))) sp_low<-matrix(c(exp(0.4652090-0.08261830),exp(0.5105223-0.06658481)),dimnames=list(c("LARGE","SMALL"))) par(mar=c(5,4,1,2),mgp=c(3,1,0)) bp5xvips<-barplot(sp_means,beside=TRUE,ylim=c(0,3),space=0.1,names.arg=c("Large","Small"),ylab="") arrows(bp5xvips,sp_upp,bp5xvips,sp_low,code=3,angle=90) ############## Factors inlfuence feeding passes when KP occur par(mfrow=c(1,2)) dir() ##MOD # 2 JFeedAttKP<-read.csv("5xvia- factors affecting feedign attempts by SP.csv") head(JFeedAttKP) JFeedAttKP$chick<-factor(JFeedAttKP$NEWID_Rep) JFeedAttKP$year<-factor(JFeedAttKP$Year) JFeedAttKP$KPAtt<-factor(JFeedAttKP$KP.ATTEMPT) JFeedAttKP$preysize<-factor(JFeedAttKP$PREYSIZE) JFeedAttKP$age<-factor(JFeedAttKP$BINARYAGE) poissJFeedAttKP1<-glmer(Feeding.attempts~preysize+age+(1|year/chick), data=JFeedAttKP, family=poisson) summary(poissJFeedAttSP1) ##PREY SIZE JfeedattemptsKPps<- lsmeans (poissJFeedAttKP1, "preysize") JfeedattemptsPSps sp_means<-matrix(c(exp(0.4652090),exp(0.5105223)),dimnames=list(c("LARGE","SMALL"))) sp_upp<-matrix(c(exp(0.4652090+0.08261830),exp(0.5105223+0.06658481)),dimnames=list(c("LARGE","SMALL"))) sp_low<-matrix(c(exp(0.4652090-0.08261830),exp(0.5105223-0.06658481)),dimnames=list(c("LARGE","SMALL"))) par(mar=c(5,4,1,2),mgp=c(3,1,0)) bp5xvips<-barplot(sp_means,beside=TRUE,ylim=c(0,3),space=0.1,names.arg=c("Large","Small"),ylab="Feeding passess") arrows(bp5xvips,sp_upp,bp5xvips,sp_low,code=3,angle=90) ##AGE JfeedattemptsKPage<- lsmeans (poissJFeedAttKP1, "age") JfeedattemptsKPage sp_means<-matrix(c(exp(0.7311760),exp(0.8415436)),dimnames=list(c("MOBILE","NESTLING"))) sp_upp<-matrix(c(exp(0.7311760+0.09908247),exp(0.8415436+0.08917856)),dimnames=list(c("MOBILE","NESTLING"))) sp_low<-matrix(c(exp(0.7311760-0.09908247),exp(0.8415436-0.08917856)),dimnames=list(c("MOBILE","NESTLING"))) par(mar=c(5,4,1,2),mgp=c(3,1,0)) bp5xviage<-barplot(sp_means,beside=TRUE,ylim=c(0,3),space=0.1,names.arg=c("Mobile","Nestling"),ylab="") arrows(bp5xviage,sp_upp,bp5xviage,sp_low,code=3,angle=90)