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

• You would likely need to use the delta method for this. Is there some reason why examining the coefficient of $\beta_1$ is not sufficient? Aug 10, 2021 at 16:44
• Okay thanks for starting this off! What is the delta method and how do I use it? And why would examining the coefficient for beta_1 give me what I want? Aug 10, 2021 at 16:49

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)

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.

• Thank you Demetri, this looks wonderful. I'm curious about the statement "Caution, this approach will yield answers below 0 or above 1.". Is there another approach that would not have this problem? Aug 10, 2021 at 18:55
• Why use this method when glm already provides the covariance matrix of the parameter estimates (and the OP has already obtained it)?
– whuber
Aug 10, 2021 at 21:09
• @whuber I don’t think this is a particularly good approach, but it is what OP asked. I’m not sure what you mean when you say glm provided the covariance matrix for parameters. I extract that via vcov but the variance OP seeks can’t be obtained directly from glm Aug 10, 2021 at 21:53
• @StatsSorceress That initial diff command is taking the difference in the probabilities. Honestly, this approach is overkill. You have a binary outcome with a binary group indicator. The easiest thing to do would be to do a test of proportions and get a confidence interval. I'm all for "all tests are just regressions" but this seems particularly overkill to me. Aug 10, 2021 at 22:29
• This rule is also (better?) known as "Gauss' law of error propagaton". It should be noted, though, that this is not an equality, but an approximation. I would thus suggest to replace the equality sign with an approximate sign (\approx in LaTeX). Aug 11, 2021 at 6:50

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

• A third possibility could be to use profile likelihood Nov 12, 2021 at 14:08