Computing prediction intervals for logistic regression

Question

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
)
print(head(students))

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 <- as.data.frame(predict(data.model, se.fit=TRUE, type='link'))

    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(head(students))
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')])

Answer 1

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.

Answer 2

I agree with @Greg Snow that there should not be a distinction between prediction interval and confidence interval for binary outcome models as in linear models. He noted the difference between a binomial outcome (y successes out of j trials among n groups) and a binary outcome (success or failure of one subject among n subjects). The former may have both a prediction interval for an individual proportion and a confidence interval for the population mean of proportions, while the latter has a confidence interval for the population probability.

Additionally, I discovered that @carbocation incorrectly implemented Collett's equation for calculating the variance and standard error of a predicted logit. A correct implementation should lead to identical results to glm() predictions. The equation (3.11) by Collett (2002, p. 98) is given in scalar, which might the reason to cause a confusion and result in the error.

Collett, D. (2002). Modelling binary data (2nd ed.). Chapman and Hall/CRC.

In the context of @carbocation R scripts, the correct codes to generate the standard error of a predicted logit are

# Prediction:
this.student.prediction <- sum(this.student.predictors * coef(data.model)) # linear combination of logistic model coefficients
var.prediction <- this.student.predictors %*% model.vcov %*% t(this.student.predictors) # matrix multiplication with dimensions (1*p) (p*p) (p*1) where p is the number of predictors including the intercept, which results in a scalar (vector of length 1)
se.prediction <- sqrt(var.prediction) # if not a scalar, se = sqrt(diag(vcov)) where vcov is a square matrix of variance-covariance. 

This procedure of matrix multiplication to generate standard error of a predicted value applies to many other types of models. To address the question in the comment by @Rafael, a standard error is simply the standard deviation of an estimate for a population parameter. When using individual scores to estimate a population mean, the standard error of a population mean estimate is created by dividing the standard deviation of individual scores by sqrt(n) where n is the sample size of individual scores. For an estimated model, however, model estimates, such as coefficients and predicted logits of probabilities according to a logistic regression, are already of the population parameters and are not a sample of individual scores. Thus, the standard deviation of model estimates, such as coefficients and predicted outcomes, are their standard error. Such uncertainty is from sampling errors of population parameters instead of from variability among individual subjects. This is also why a bootstrapped estimate of the standard error of a statistic with an unknown distribution uses the standard deviation of the observed statistics among bootstrapped samples, which does not shrink when taking more bootstrapped samples.