Logistic regression: how to compute a prediction interval

Question

Suppose I have a simple logistic regression model:

$log(p/(1-p)) = \beta_0 + \beta_1x$

Then I know:

$p/(1-p) = e^{\beta_0 + \beta_1x}$

and

$p = e^{\beta_0 + \beta_1x}/(1 + e^{\beta_0 + \beta_1x})$

How can I obtain a confidence interval for p?

That is, how can I obtain a 95% confidence interval on the probability of success (would that be called a prediction interval? Maybe I'm using the wrong term)?

I've looked in a few places, but I can't find the formula. For example, I've looked:

• here (which indicates such an interval is not of interest because each value in the response can only be 0 or 1, but that's not what I'm interested in),

• and here (and again, no formula, but similar question)

Answer 1

would that be called a prediction interval? Maybe I'm using the wrong term)

That would not be a prediction interval. A prediction interval would incorporate uncertainty in the data generation. Its a bit useless for a binary logistic regression since we know the outcome will either be 0 or 1. A prediction interval may be more useful when you have trial data (e.g. I predict between 8 and 12 events out of the 20).

Onto calculating the confidence interval. First, start by computing $\operatorname{Var}(\operatorname{logit}(p))$. Since $\operatorname{logit}(p) = \beta_0 + \beta_1x$

$$\operatorname{Var}(\operatorname{logit}(p)) = \operatorname{Var}(\beta_0) + x^2 \operatorname{Var}(\beta_1) + 2x\operatorname{Cov}(\beta_0, \beta_1)$$

You can find these quantities from the covariance matrix estimated from the model. Take the square root of this quantity and you have the standard deviation of $\operatorname{logit}(p)$. Let's call this quantity $\sigma$ for now.

A confidence interval for $\operatorname{logit}(p)$ is $(\theta_L, \theta_U) = (\operatorname{logit}(p) - 1.96 \sigma,\operatorname{logit}(p) + 1.96 \sigma )$. Invert each endpoint to obtain a confidence interval for $p$.

Here is how one might perform this in R.

library(tidyverse)

#Simulate some data to run a regression on
x = rnorm(100)
eta = -0.8 + 0.2*x
p = plogis(eta)
y = rbinom(100, 1, p)

# Fit a model
model = glm(y~x, family = binomial())
Bigma = vcov(model)

# Formula above for standard deviation of logit p.  Vectorize for ease of computation
sig =Vectorize(function(x) sqrt(Bigma[1,1] + x^2*Bigma[2,2] + 2*x*Bigma[1,2]))

# New xs to make prediction on.
new_x = seq(-2, 2, 0.01)
logit_p = predict(model, newdata=list(x=new_x))
theta_L = logit_p - 1.96*sig(new_x)
theta_U = logit_p + 1.96*sig(new_x)

# Invert the estimates
p_L = plogis(theta_L)
p = plogis(logit_p)
p_U = plogis(theta_U)

d = data.frame(x=new_x, p_L, p, p_U)

# Now compare with builtin tools

fits = predict(model, newdata = list(x=new_x), se.fit = T) %>% 
       bind_cols() %>% 
       mutate(x=new_x) %>% 
       mutate(p = plogis(fit),
              p_L = plogis(fit - 1.96*se.fit),
              p_U = plogis(fit + 1.96*se.fit))



ggplot()+
  # R's computation of the confidence interval in Black
  geom_line(data=fits, aes(x, p))+
  geom_ribbon(data=fits, aes(x=x, ymin=p_L, ymax=p_U), alpha = 0.5)+
  # Out computation of the confidence interval in Red
  geom_line(data=d, aes(x, p), color = 'red')+
  geom_ribbon(data=d, aes(x=x, ymin=p_L, ymax=p_U), alpha = 0.5, fill='red')

Here is an example output

The red and the black are right overtop one another, meaning the manual procedure reproduces the output of R's built in methods for calculating the standard error of the fit. You can run this code for multiple random seeds and verify they produce the same estimates (at least up to a few decimal places I imagine).

Answer 2

There are three or four options for confidence intervals.

You have a non-linear function of coefficients in your third equation, and you can use the delta method to calculate the approximate variance of that function. This may get you values outside $[0,1]$ interval, but that's not as big a deal as people make it out to be.

You can also predict the index function part ($\beta_0 + \beta_1\cdot x)$, which is a linear function of coefficients, get the lower and upper bounds of the CIs for that, and then use the inverse logit to get the CIs on the probability scale. These should lie inside $[0,1]$. This is what Demetri Pananos suggested above.

Third, you can bootstrap the estimation and the prediction part together.

But since your $x$ is binary, you can just fit a linear probability model with heteroskedasticity-robust variance and be done with it. This is arguably the easiest path of all and will give you the same point estimates and very similar CIs.

Here's annotated Stata output showing how to do all four:

. set seed 9191945

. sysuse auto
(1978 automobile data)

. gen high_mpg = mpg>22

. logit foreign i.high_mpg, nolog

Logistic regression                                     Number of obs =     74
                                                        LR chi2(1)    =  14.76
                                                        Prob > chi2   = 0.0001
Log likelihood = -37.652732                             Pseudo R2     = 0.1639

------------------------------------------------------------------------------
     foreign | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
  1.high_mpg |   2.077817   .5699326     3.65   0.000     .9607695    3.194864
       _cons |  -1.767662   .4089589    -4.32   0.000    -2.569207   -.9661172
------------------------------------------------------------------------------

. /* Delta Method */
. margins high_mpg

Adjusted predictions                                        Number of obs = 74
Model VCE: OIM

Expression: Pr(foreign), predict()

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
    high_mpg |
          0  |   .1458333   .0509424     2.86   0.004     .0459881    .2456785
          1  |   .5769231   .0968907     5.95   0.000     .3870209    .7668253
------------------------------------------------------------------------------

. /* inverse logit */
. margins high_mpg, predict(xb)

Adjusted predictions                                        Number of obs = 74
Model VCE: OIM

Expression: Linear prediction (log odds), predict(xb)

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
    high_mpg |
          0  |  -1.767662   .4089589    -4.32   0.000    -2.569207   -.9661172
          1  |   .3101549   .3969581     0.78   0.435    -.4678687    1.088179
------------------------------------------------------------------------------

. transform_margins invlogit(@)
----------------------------------------------
             |         b         ll         ul
-------------+--------------------------------
    high_mpg |
          0  |  .1458333   .0711467   .2756551
          1  |  .5769231   .3851208   .7480386
----------------------------------------------

. /* Bootstrap */
. capture program drop savemargins

. program savemargins, rclass
  1.         logit foreign i.high_mpg, nolog
  2.         margins high_mpg, post
  3. end

. bootstrap _b, reps(200): savemargins
(running savemargins on estimation sample)

Bootstrap replications (200)
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 
..................................................    50
..................................................   100
..................................................   150
..................................................   200

Adjusted predictions                                       Number of obs =  74
                                                           Replications  = 200

------------------------------------------------------------------------------
             |   Observed   Bootstrap                         Normal-based
             | coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
    high_mpg |
          0  |   .1458333   .0482008     3.03   0.002     .0513615    .2403052
          1  |   .5769231   .0963147     5.99   0.000     .3881497    .7656965
------------------------------------------------------------------------------

. /* Het-Robust LPM */
. regress foreign ibn.high_mpg, noconstant robust

Linear regression                               Number of obs     =         74
                                                F(2, 72)          =      21.23
                                                Prob > F          =     0.0000
                                                R-squared         =     0.4398
                                                Root MSE          =     .41375

------------------------------------------------------------------------------
             |               Robust
     foreign | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
    high_mpg |
          0  |   .1458333   .0516451     2.82   0.006     .0428808    .2487859
          1  |   .5769231   .0982272     5.87   0.000     .3811108    .7727353
------------------------------------------------------------------------------

All the point estimates are identical, and the CIs are generally very similar across all four approaches, even with 74 observations.