##LOAD MCMC OUTPUT FOR CHICKS

#sum the capture-recapture matrix by rows to get the number of times each individual was re-sighted
obs_byind <- rowSums(obsmatrix)

#get the posterior mode or mean from mcmc output for each re-sight probability in the capture-recapture matrix
predicMCMC_byind = mat.or.vec(4788, 1)
for(i in 1:4788){
dpredicMCMC_byind <- density(REMC.mcmc[,i+1])
#to get posterior mode
predicMCMC_byind[i] <- max(dpredicMCMC_byind$x[which.max(dpredicMCMC_byind$y)], 0)
#to get posterior mean
predicMCMC_byind[i] <- mean(REMC.mcmc[,i+1])
}

#turn predicMCMC_byind vector into an i (individual) by z (sampling event) matrix
predic_byind <- matrix(predicMCMC_byind, nrow=266, ncol=18, byrow=FALSE)

#sum rows to get parameters for the Poisson-binomial distribution
predic_byind_lambda <- rowSums(predic_byind)

##LOAD PACKAGE 'poibin'

#get 95% prediction intervals from Poisson-binomial dist
quantiles_byind = mat.or.vec(266, 2)
for(i in 1:266){
quantiles_byind[i, ] <- qpoibin(c(.025, .975), predic_byind[i, ])
}

#put the posterior mode/mean, 95% prediction intervals, and observed data into the same matrix
forplot_byind <- cbind(predic_byind_lambda, quantiles_byind, obs_byind)
#order the matrix by mode/mean
forplot_byind <- forplot_byind[order(predic_byind_lambda),]

#for each individual, record a 1 if number of re-sights falls outside of the prediction interval, 0 otherwise
postcheck_byind = mat.or.vec(266, 1)
a = mat.or.vec(266, 1)
b = mat.or.vec(266, 1)
for(i in 1:266){
if(forplot_byind[i,4] > forplot_byind[i, 3]){
a[i] <- 1
}
if(forplot_byind[i, 4] < forplot_byind[i, 2]){
b[i] <- 1
}
postcheck_byind[i] <- a[i] + b[i]
}

#set plotting parameters
quartz.options(width=5.5, height=5)
par(mfcol=c(2, 2), mar=c(3, 4, 1, 1))

#plot for figure 4b
plot(forplot_byind[,1], type='l', ylim=c(0, 12), xlab=" ", xaxt="n", ylab=" ", cex.axis=.7, cex=.7)
segments(1:266, forplot_byind[,2], 1:266, forplot_byind[,3],  lwd=6, col="gray", lend="butt")
mtext(side=1, line=1, "Each dot is an individual", cex=.7)
mtext(side=2, line=2, "Number of times chick was resighted", cex=.7)
mtext(side=2, line=2, at=14.5, "b)", cex=.75, las=2)
points(forplot_byind[,4], pch=16)