# 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? Apr 4, 2017 at 3:00
• Jan 16 at 11:22

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.

• I'm looking for a 95% prediction interval of a prediction which is in log-odds space. Later I transform that to probability space. A 100% prediction interval would never be interesting for any procedure, right? For example, a 100% prediction interval for linear regression would include -Inf to Inf... At any rate, as you can see in my code, the prediction interval is calculated in log odds space, which is then transformed into probability space later. So I don't think my question is pointless. Apr 17, 2012 at 0:34
• The log-odds can be converted to a probability and you can compute a confidence interval on the probability (or the log-odds). But a prediction interval is on the response variable which is 0 or 1. If your outcome is survival with 0=dead and 1=alive, then you can predict the probability of being alive for a given set of covariates and compute a confidence interval on that probability. But the outcome is 0/1, you can't have a patient that is 62% alive it has to be 0 or 1, so the only possible prediction intervals are 0-0, 0-1, and 1-1 (which is why most people stick to confidence intervals). Apr 17, 2012 at 0:57
• If you have a situation where the response is binomial (which could be an aggregate of 0-1s under the same conditions), then a prediction interval may make sense. May 24, 2012 at 5:49
• Logistic regression is regression of a probability, trying to model the probability of some event as a function of regressor variables. Prediction intervals in this setting is taken as intervals on the probability scale, or the log-odds scale, so makes perfect senes. Jun 24, 2015 at 14:23
• @Cesar, the prediction interval formula is derived by assuming that Y is normally distributed about the line, but in logistic regression we don't have a normal distribution, we have a Bernoulli or Binomial. Applying the formulas on that page would lead to either a confidence interval (can already do this) or an artificially widened confidence interval that does not meet the definition of a prediction interval (predicting actual outcomes on the original outcome scale). Like Glen_b mentioned, a prediction interval may make sense if the outcome is truly binomial. Sep 16, 2016 at 15:27

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.

• This is fantastic. Thank you for your close read and for the corrections. We should probably edit my question to add a disclaimer and point to your answer. (Or possibly correct it in-line, but my guess is that it will make more sense for posterity to leave it incorrect there with a warning.) Dec 2, 2022 at 20:51
• Sounds good to me. The use of 1.96 instead a more accurate qnorm(0.975) might be another source of disparity from other's results. The {ggeffects} package will calculate these intervals via ggeffects::ggpredict(), ggeffects::ggemmeans(), geffects::ggeffect(). The definition of each statistic varies. Similar packages include {marginaleffects} and {margins} with {prediction}. Dec 9, 2022 at 20:13
• Disclaimer: I have not tested the correction I made. So I encourage you to do so and see whether the resulting CI matches with those built from glm() and predict(). Dec 9, 2022 at 20:19