Logistic regression: how to compute a prediction interval 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)
 A: 
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).
