# Computing prediction intervals for logistic regression

I would like to understand how to generate prediction intervals for logistic regression estimates.

I was advised to follow the procedures in Collett's Modelling Binary Data, 2nd Ed p.98-99. After implementing this procedure and comparing it to R's predict.glm, I actually think this book is showing the procedure for computing confidence intervals, not prediction intervals.

Implementation of the procedure from Collett, with a comparison to predict.glm, is shown below.

I would like to know: how do I go from here to producing a prediction interval instead of a confidence interval?

#Derived from Collett 'Modelling Binary Data' 2nd Edition p.98-99
#Need reproducible "random" numbers.
seed <- 67

num.students <- 1000
which.student <- 1

#Generate data frame with made-up data from students:
set.seed(seed) #reset seed
v1 <- rbinom(num.students,1,0.7)
v2 <- rnorm(length(v1),0.7,0.3)
v3 <- rpois(length(v1),1)

#Create df representing students
students <- data.frame(
intercept = rep(1,length(v1)),
outcome = v1,
score1 = v2,
score2 = v3
)

predict.and.append <- function(input){
#Create a vanilla logistic model as a function of score1 and score2
data.model <- glm(outcome ~ score1 + score2, data=input, family=binomial)

#Calculate predictions and SE.fit with the R package's internal method
# These are in logits.

predictions$actual <- input$outcome
predictions$lower <- plogis(predictions$fit - 1.96 * predictions$se.fit) predictions$prediction <- plogis(predictions$fit) predictions$upper <- plogis(predictions$fit + 1.96 * predictions$se.fit)

return (list(data.model, predictions))
}

output <- predict.and.append(students)

data.model <- output[[1]]

#summary(data.model)

#Export vcov matrix
model.vcov <- vcov(data.model)

# Now our goal is to reproduce 'predictions' and the se.fit manually using the vcov matrix
this.student.predictors <- as.matrix(students[which.student,c(1,3,4)])

#Prediction:
this.student.prediction <- sum(this.student.predictors * coef(data.model))
square.student <- t(this.student.predictors) %*% this.student.predictors
se.student <- sqrt(sum(model.vcov * square.student))

manual.prediction <- data.frame(lower = plogis(this.student.prediction - 1.96*se.student),
prediction = plogis(this.student.prediction),
upper = plogis(this.student.prediction + 1.96*se.student))

print("Data preview:")
print(paste("Point estimate of the outcome probability for student", which.student,"(2.5%, point prediction, 97.5%) by Collett's procedure:"))
manual.prediction
print(paste("Point estimate of the outcome probability for student", which.student,"(2.5%, point prediction, 97.5%) by R's predict.glm:"))
print(output[[2]][which.student,c('lower','prediction','upper')])

• A basic question, why is sqrt(sum(model.vcov * square.student)) assumed as the standard error? Isn't it the standard deviation and needs to be divided by sqrt(n)? If so, which n should be used, n used to fit the model or n of the new data frame used to predict? – Rafael Apr 4 '17 at 3:00

Prediction intervals predict where the actual response data values are predicted to fall with a given probability. Since the possible values of the response of a logistic model are restricted to 0 and 1, the 100% prediction interval is therefore $0 <= y <= 1$. No other intervals really make sense for prediction with logistic regression. Since it is always the same interval it generally is not interesting enough to generate or discuss.