3
$\begingroup$

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)

$\endgroup$
2
  • 3
    $\begingroup$ A confidence interval for p is not a prediction interval. That's why you're finding answers that are not relevant. $\endgroup$
    – Glen_b
    Commented Aug 10, 2021 at 0:10
  • 2
    $\begingroup$ Are you after a CI at some specific x or some kind of confidence band across an interval of x? If a band, are you after a pointwise interval or some overall (simultaneous/interval-wise) coverage? $\endgroup$
    – Glen_b
    Commented Aug 10, 2021 at 5:00

2 Answers 2

5
$\begingroup$

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 = TRUE) %>% 
       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

enter image description here

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).

$\endgroup$
8
  • $\begingroup$ This is great! Although you've assumed x is continuous here - what if it's binary as well? $\endgroup$ Commented Aug 10, 2021 at 1:45
  • $\begingroup$ @StatsSorceress Formulae apply still. $\endgroup$ Commented Aug 10, 2021 at 1:45
  • $\begingroup$ Okay, so I'm guessing it means I have two confidence intervals: the case when x=0 and the case when x=1. If x=0, then Var(logit(p)) would be Var(beta_0), which is a constant. I think then new_x would be 0, and sig(new_x) reduces to sqrt(Bigma[1,1]). Then I would still invert the estimates to get p_L and p_U, and those would be my lower and upper bounds, respectively. $\endgroup$ Commented Aug 10, 2021 at 1:54
  • $\begingroup$ And if x=1, then Var(logit(p)) = Var(beta_0) + Var(beta_1) + 2Cov(beta_0, beta_1), and if I square root that, I get sig(new_x), and again I would invert the estimates to get p_L and p_U. Does that make sense? $\endgroup$ Commented Aug 10, 2021 at 1:56
  • 2
    $\begingroup$ @swag2198 Yes, this is frequentist. $\endgroup$ Commented Aug 10, 2021 at 13:45
2
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ Thank you Dimitriy, but a formula has not been provided here. $\endgroup$ Commented Aug 10, 2021 at 1:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.