ospreyData<-read.csv("OspreyFeedingByEvent.csv",header=T) #Analysis of the effect of a parent being present table(ospreyData$Dead, ospreyData$Parent) fisher.test(table(ospreyData$Dead, ospreyData$Parent)) #Install lme4 package if necessary chooseCRANmirror() install.packages("lme4",dependencies=T) library("lme4") #Fit model looking for an effect of chick age on the probability of prey being delivered alive model1<-glmer(Dead~ChickAge + (1|Pair)+ (ChickAge|Pair),family=binomial,data= ospreyData) summary(model1) #Calculate effect on odds per day with back-transformed Wald 95% confidence intervals exp(-0.06752) exp(-0.06752+c(-1.96,1.96)* 0.02531) #Now fit a model looking for an effect of chick age on the probability of the non-delivering parent being present at the nest model2<-glmer(Parent ~ChickAge +(1|Pair)+ (ChickAge|Pair),family=binomial,data= ospreyData) summary(model2) #Now create plots #Read in data pooled by day OspreyFeedingByDay<-read.csv("OspreyFeedingByDay.csv",header=T) #Use this to create a data frame with data pooled in bins of five days for each pair totals<-dead<-day<-parent<-NULL pair<-unique(ospreyData$Pair)[1] for (i in 1:20){ for (j in 1:length(unique(ospreyData$Pair))){ tempData<-OspreyFeedingByDay[as.numeric(OspreyFeedingByDay $Pair)==j,] tempData2<-tempData[tempData$ChickAge>((i-1)*5)&tempData$ChickAge<=(i*5),] totals<-c(totals,sum(tempData2$NumberOfFish)) dead <-c(dead,sum(tempData2$NumberDead)) pair<-as.factor(c(pair,unique(ospreyData$Pair)[j])) day<-c(day,i*5-2.5) parent<-c(parent,sum(tempData2$ParentPresent)) levels(pair)<-levels(ospreyData$Pair) } } pair<-pair[-1] plotData<-data.frame(day,totals, dead, pair,parent) #Produce the plots layout(matrix(c(1,2),2,1),widths=4,heights=c(3.75,4)) par(mar=c(2,4,1,1)) plot(jitter(plotData$day),plotData$dead/plotData$totals*100,col=0, ylab="% prey dead on arrival at nest", xlab="",ylim=c(0,110)) rect(44, -20, 66, 140, density = 20, col = "grey", border = "transparent") points(jitter(plotData$day),plotData$dead/plotData$totals*100,col=plotData$pair) lines(seq(0,100,length=1000),100-100*exp(-0.74322+seq(0,100,length=1000)*-0.06486)/(1+exp(-0.74322+seq(0,100,length=1000)*-0.06486)),lwd=2) lines(unique(plotData$day),100*tapply(plotData$dead,plotData$day,sum)/tapply(plotData$totals,plotData$day,sum),lty=2,lwd=2,col=2) text(2,108,"a)") par(mar=c(4,4,1,1)) plot(jitter(plotData$day),plotData$parent/plotData$totals*100,col=0, ylab="% other parent present at nest", xlab="Youngest chick age (days)",ylim=c(0,110)) rect(44, -20, 66, 140, density = 20, col = "grey", border = "transparent") points(jitter(plotData$day),plotData$parent/plotData$totals*100,col=plotData$pair) lines(seq(0,100,length=1000),100*exp(8.54079 +seq(0,100,length=1000)*-0.13983)/(1+exp(8.54079 +seq(0,100,length=1000)*-0.13983)),lwd=2) lines(unique(plotData$day),100*tapply(plotData$parent,plotData$day,sum)/tapply(plotData$totals,plotData$day,sum),lty=2,lwd=2,col=2) text(2,108,"b)")