How can I obtain a confidence interval for the difference between two probabilities from a binary logistic regression model? I saw this question here, but it doesn't have a clear answer.
Suppose I have a simple logistic regression model with binary $x:$
$$\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})$$
So if $x=0,$ then the model becomes:
$$\log(p/(1-p)) = \beta_0$$
and if $x=1,$ then the model becomes:
$$\log(p/(1-p)) = \beta_0 + \beta_1$$
To obtain a confidence interval for $p$ when $x=0,$ I did this:
model <- glm(y~., family=binomial(), data)

#For x=0
Bigma = vcov(model)
sig = sqrt(Bigma[1,1])

logit_p = coef(model)[1][[1]] #The intercept

theta_L = logit_p - 1.96*sig
theta_U = logit_p + 1.96*sig

p_L = plogis(theta_L) 
p_U = plogis(theta_U)

The confidence interval is (p_L, p_U).
Then a confidence interval for p when x=1:
sig_x1 = sqrt(Bigma[1,1] + Bigma[2,2] + 2*Bigma[1,2])

logit_p_x1 = coef(model)[1][[1]] + coef(model)[2][[1]] #beta_0 + beta_1

theta_L_x1 = logit_p_x1 - 1.96*sig_x1
theta_U_x1 = logit_p_x1 + 1.96*sig_x1

p_L_x1 = plogis(theta_L_x1) 
p_U_x1 = plogis(theta_U_x1)

The confidence interval is (p_L_x1, p_U_x1).
Now I would like a confidence interval for the difference in probability of success when $x=0$ and $x=1.$
I can obtain the point estimate of the difference:
$$p_{x_1} - p_{x_0} = [e^{\beta_0 + \beta_1} - e^{\beta_0}]/[(1+e^{\beta_0 + \beta_1})(1+e^{\beta_0})]$$
I know the next step is to compute the standard error of the difference, but I don't know how to do this.
Question: What is the formula and R code for a confidence interval for the difference in the two probabilities when $x=0$ and when $x=1?$

There are two principal approaches:

*

*Estimate the variance $\hat{\sigma}^2$ of $\theta=p_{x1}-p_{x0}$ and assume that $\theta$ is normally distributed. Then the confidence interval is $\pm z_{1-\alpha/2} \hat{\sigma}$, where $z_{1-\alpha/2}$ is the quantile of the standard normal distribution (qnorm(1-alpha/2) in R), which is 1.96 for a 95% interval ($\alpha=0.05$).

*Directly estimate a non-parametric (and possibly non-symmetrc) confidence interval with the bootstrap method.

The variance for method 1. can be estimated in different ways. One is the solution based on Gauss' error propagation law suggested by @demetri-pananos. Another method would be the Jackknife, which consists in estimating the parameter $\theta_{(i)}$ $n$ times, each time with one observable left out, and compute the Jackknife Variance therefrom:
\begin{align}
\sigma_{JK} &= \sqrt{\frac{n-1}{n}\sum_{i=1}^n \Big(\theta_{(i)}-\theta_{(.)}\Big)^2}\\
 \mbox{ with }&\quad \theta_{(.)}=\frac{1}{n}\sum_{i=1}^n\theta_{(i)} \nonumber
\end{align}
Method 2. is similar to the Jackknife, but the parameter is estimated several times from observables drawn with replacement. From these estimates, the confidence interval can be estimated in different ways. In comparative studies, the "Bias Corrected Accelerated Bootstrap" had the best coverage probability. See section 6 of this report for R code how to compute it and different comparative studies (section 5.2 of the same report explains the Jackknife and lists R code how to compute it).
