gistfile1.r

t <- 10 # number of times
n <- 100 # number of individuals 'n'
x <- rnorm(t,0,1) # environmental variable 'x'
mu <- 4 # average survival response to 'x'
s <- 2 # variation among individuals in survival response to 'x'
b0 <- 1.5
beta <- rnorm(n,mu,s) # survival response for individual to 'x'

# Deterministic survival probabilty for each individual at each time
theta <- exp(b0+x%*%t(beta))/(1+exp(b0+x%*%t(beta))) 
# Simulated data for the population
y <- matrix(rbinom(n*t,1,theta),nrow=t)

# response matrix for GLM: cbind(number of successes, number of failures)
response <- cbind(matrix(y,ncol=1),rep(1,n*t)-matrix(y,ncol=1))
model <- summary(glm(response~rep(x,n),binomial(link=logit)))

# Bootstrapped parameter predictions for model and individual effects
bpred <- rnorm(10000,model$coef[2,1],model$coef[2,2])
bind <- rnorm(10000,mu,s)

# Survival predictions
predx <- seq(-5,5,0.1)
# "probability of an average individual" prediction
predya <- exp(model$coef[1,1]+predx*model$coef[2,1])/(1+exp(model$coef[1,1]+predx*model$coef[2,1]))
# same as above, prediction from data
predy2 <- exp(model$coef[1,1]+x*model$coef[2,1])/(1+exp(model$coef[1,1]+x*model$coef[2,1]))
# "average probability" prediction
predyi <- exp(b0+predx*mu)/(1+exp(b0+predx*mu))

# Plot 1. Simulated data, model prediction (pink) and deterministic average (blue)
par(mar=c(4.5,4.5,0.5,0.5),mgp=c(2,0.7,0))
plot(jitter(y)~matrix(rep(x,n),nrow=t),yaxt="n",ylab="Survive",xlab="environment 'x'",col="darkgrey")
axis(2,yaxp=c(0,1,1),las=1)
abline(h=c(0,1))
points(exp(model$coef[1,1]+x*model$coef[2,1])/(1+exp(model$coef[1,1]+x*model$coef[2,1]))
 ~x,pch=19,cex=2,col="deeppink")
points(exp(b0+x*mu)/(1+exp(b0+x*mu))~x,ylim=c(0,1),pch=19,cex=2,col="royalblue")

# Plot 2a. Bootstrapped parameter estimates from model fit (pink) and simulated data (blue)
# If I was projecting results for population growth, this is what I would use
par(mfcol=c(1,2),mar=c(4.5,4.5,0.5,0.5),mgp=c(2,0.7,0))
plot(density(bpred)$x,density(bpred)$y,lwd=2,col="deeppink",type="l",xlim=c(0,8),xlab="Bootstrapped beta (n=10000)",ylab="density")
lines(density(bind)$x,density(bind)$y,lwd=2,col="royalblue")
legend("topright",c("GLM fit","'Real' response"),lwd=2,col=c("deeppink","royalblue"),bty="n")

# Plot 2b. Predictions are not good because the GLM is predicting 
# the probability that an average individual will survive (x-axis solid line)
# not the average probability that an individual will survive (y-axis solid line)
plot(predya,predyi
 ,xlim=c(0,1),ylim=c(0,1),xlab="Predicted survival",ylab="'Actual' survival",type="l")
points(rep(predy2,n),matrix(theta,ncol=1),col="darkgrey")
abline(0,1,lty=2)
legend("bottomleft",c(""),bty="n",title="Overpredicting survival at low 'x'",text.col="deeppink")
legend("topright",c(""),bty="n",title="Underpredicting survival at high 'x'",text.col="deeppink")