How can I obtain a confidence interval for the difference between two probabilities from a binary logistic regression model?

Question

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?$

Answer 1

The delta method states

$$ \operatorname{Var}(g(X)) = [g'(X)]^2 \operatorname{Var}(X)$$

Because this problem involves two parameters, we can extend this to the multivariate delta method

$$ =\nabla g^T \, \Sigma \, \nabla g $$

Here,

$$ g = \left[e^{\beta_{0}+\beta_{1}}-e^{\beta_{0}}\right] /\left[\left(1+e^{\beta_{0}+\beta_{1}}\right)\left(1+e^{\beta_{0}}\right)\right] $$

and $\Sigma$ is the variance covariance matrix from your model. $\nabla g$ is...gross. I'm not going to do that by hand, and computer algebra while fast yields a mess of symbols. You can however use autodifferentiation compute the gradient. Once you calculate the variance, then its simply your estimate of the difference in probs plus/minus 1.96 times the standard deviation (root of the variance). Caution, this approach will yield answers below 0 or above 1.

We can do this in R in the following way (note you need to install the autodiffr package).

library(autodiffr)

g = function(b)  (exp(b[1] + b[2]) - exp(b[1])) / ((1+ exp(b[1] + b[2]))*(1+exp(b[1])))

x = rbinom(100, 1, 0.5)
eta = -0.8 + 0.2*x
p = plogis(eta)
y = rbinom(100, 1, p)

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

grad_g = makeGradFunc(g)
nabla_g = grad_g(coef(model))


se = as.numeric(sqrt(nabla_g %*% Bigma %*% nabla_g))


estimate = diff(predict(model, newdata=list(x=c(0, 1)), type='response'))

estimate + c(-1, 1)*1.96*se

Repeating this procedure for this modest example shows that the resulting confidence interval has near nominal coverage, which is a good thing, but I imagine things would become worse as the probabilities approach 0 or 1.

Answer 2

There are two principal approaches:

1. 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$).
2. 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).