Confidence interval for difference between two predicted probabilities in R In R, I have estimated a logistic regression and calculated two predicted probabilities (with 95% confidence intervals) using the code shown:
set.seed(1234)

x = runif(100, 0, 1)
y = rbinom(100, size=1, prob = x)

model = glm(y ~ x, family = binomial("logit"))

newdata = data.frame(x = c(.25, .75))

predicted.probs = predict(model, newdata, type="response", se.fit = T)

upper.interval  = predicted.probs$fit + 1.96 * predicted.probs$se.fit
lower.interval  = predicted.probs$fit - 1.96 * predicted.probs$se.fit

Hmisc::errbar(newdata$x, predicted.probs$fit, upper.interval, lower.interval)

Rather than plotting these two predicted probabilities separately, I want to plot the difference between them, along with a 95% confidence interval for that difference. I know that the point estimate for that difference is simply the difference between the two predicted probabilities, but I do not know how to determine the confidence interval. I could simply add widths of the two original confidence intervals, but this seems wrong to me given that the two predicted probabilities are based on the same model and not independent.
Is there a way in R to calculate a confidence interval for this difference between two predicted probabilities?
 A: One way to do this will be profile likelihood. If we have a parameter vector $\psi$, profile likelihood is usually calculated for one of the components of $\psi$, but it can be defined for any parametric function of $\psi$. Below is a definition, suppose $L(\psi)$ is the likelihood function and interest (or focus) is on a scalar function $\theta = \theta(\psi)$, then
$$ L_P(\theta) = \max_{\{\psi\colon \theta(\psi)=\theta \}} L(\psi)$$
The implementations of profile likelihood in R (elsewhere?) is not of these generality, so let us make it "by hand".
The model is
$$ \DeclareMathOperator{\P}{\mathbb{P}} 
p_x= \P(Y=1 \mid X=x)= \frac1{1+e^{-\beta_0 - \beta_1 x}} $$ and the interest parameter $\theta$ is
$$ \theta = p_{0.75} - p_{0.25} $$
It doesn't look promising to try to solve the optimization symbolically, so we try numerically. This is a first attempt, so maybe we can do better. First, a plot of the (negative) profile likelihood for $\theta$, using the data simulated in the question:

the two blue lines are cutoffs for confidence intervals of 95 and 99%, respectively, based on quantiles from the reference chi-square distribution with 1 df. R code is below:
### First run code from question
library(bbmle)

make_negloglik <- function(y, x) {
   n <- length(y)
   stopifnot( n == length(x) )
   Vectorize( function(beta0, beta1) 
       sum(ifelse(y==0, log1p(exp(beta0  +  beta1*x)),
                                   log1p(exp(-beta0 - beta1*x)))) )
    }

negloglik <- make_negloglik(y, x)

mod.bb <- bbmle::mle2(negloglik,  start=list(beta0=-2, beta1=4))

mod.prof <- bbmle::profile(mod.bb)

plot(mod.prof) # Not shown 

grid <- expand.grid(beta0=seq(-2.8, -0.5, len=100),
                    beta1=seq(1.8, 7.1, len=100))
grid$negloglik <-  with(grid, negloglik(beta0, beta1)) 

P <- function(beta0, beta1, x) 1/( 1  + exp( -beta0 -beta1 * x))

theta <- function(beta0, beta1) P(beta0, beta1, 0.75) - P(beta0, beta1, 0.25)

### Adding theta as a column to data.frame grid:

grid$theta <- with(grid, theta(beta0, beta1))

profile_negloglik <- function(grid) {
    rt <- with(grid,  range(theta))
    seq_theta <- seq(rt[1], rt[2], len=201)
    delta <- diff(seq_theta[1:2])
    npl <- numeric(length=length(seq_theta))
    for (t in seq_along(seq_theta)) {
        tt <- seq_theta[t]
        npl[t] <- with(grid, min(grid[ (tt-delta/2 <= theta) & (theta <= tt + delta/2),
                                      "negloglik" ]))
        }
    return(data.frame(theta=seq_theta, npl=npl))
    }

npl_frame <- profile_negloglik(grid)

npl_min <- with(npl_frame, min(npl))

library(ggplot2)

ggplot(npl_frame, aes(theta, npl))  +  geom_line(color="red")  +
    ggtitle("Profile negative loglikelihood for theta")  +
    geom_hline(yintercept=npl_min)  +
    geom_hline(yintercept=npl_min +  
    qchisq(0.95, 1)/2, color="blue")  +
    geom_hline(yintercept=npl_min +  
    qchisq(0.99, 1)/2, color="blue")   +  ylim(52, 70)   

The idea of the code is:

*

*Define a rectangle in parameter space given by individual 99% confidence intervals (calculated by profiling with the R package bbmle)

*use expand.grid to cover the rectangle

*add to the grid data frame a column with the negative loglikelihood, another column with $\theta$

*Find the range of $\theta$ and subdivide it in many small intervals

*For each of the intervals, find the minimum negative log likelihood over the interval, and associate that with the midpoint

*finally, plot this as an approximation of the negative profile loglikelihood function of $\theta$.

As a comparison, let us also calculate an approximate 95% confidence interval using the delta method. Calculations in R:
 theta_grad <- deriv(expression(  1/( 1  + exp( -beta0 -beta1 * 0.75)) 
              -  1/( 1  + exp( -beta0 -beta1 * 0.25))), 
              c("beta0", "beta1"), function.arg=TRUE)

grad <- theta_grad(coef(model)[1], coef(model)[2])

grad
(Intercept) 
  0.4880566 
attr(,"gradient")
           beta0      beta1
[1,] -0.05555914 0.06582565

grad <- attr(grad, "gradient")
V <- vcov(model)

theta.se <- sqrt( grad %*% V %*% t(grad)  )

 ( CI <- c(0.4881 -2*theta.se, 0.4881  + 2*theta.se ) )
[1] 0.3154351 0.6607649

which is quite close to the profile interval.
A: The emmeans package provides an easy and reliable way to calculated an asymptotic confidence interval for the difference in probabilities via the multivariate delta method (see also @kjetilbhalvorsen's answer for computational details):
library(emmeans)

x_pred <- c(0.75, 0.25)
emm_resp <- emmeans(model, ~ x, at = list(x = x_pred), regrid = "response")
# infer = c(TRUE, TRUE) gives test statistic w/ corresponding p-value, and CI
pairs(emm_resp, infer = c(TRUE, TRUE))
# contrast      estimate     SE  df asymp.LCL asymp.UCL z.ratio p.value
# x0.75 - x0.25    0.488 0.0863 Inf     0.319     0.657   5.653  <.0001
# 
# Confidence level used: 0.95 

A: As briefly discussed in my answer to the other poster, you can simply calculate the overall standard error from the two standard errors and use the normal approximation to get the confidence interval. The overall standard error is simply sqrt(se_0^2 + se_1^2). The following function implements this (in a very crude fashion):
get_diff_ci <- function(pred) {
  diff <- pred$fit[2] - pred$fit[1]
  diff_se <- sqrt(pred$se.fit[1]^2 + pred$se.fit[2]^2)
  upper <- diff + 1.96 * diff_se
  lower <- diff - 1.96 * diff_se
  out_dat <- data.frame(diff=diff, lower=lower, upper=upper)
  return(out_dat)
}

So in your example it would be:
set.seed(1234)

x = runif(100, 0, 1)
y = rbinom(100, size=1, prob = x)

model = glm(y ~ x, family = binomial("logit"))

newdata = data.frame(x = c(.25, .75))

predicted.probs = predict(model, newdata, type="response", se.fit = T)

get_diff(predicted.probs)

To "prove" that this is appropriate, a small monte-carlo simulation using the original example given:
set.seed(1234)

n_repeats <- 10000

out <- vector(mode="list", length=n_repeats)
for (i in 1:n_repeats) {
  x = runif(100, 0, 1)
  y = rbinom(100, size=1, prob = x)
  
  model = glm(y ~ x, family = binomial("logit"))
  
  newdata = data.frame(x = c(.25, .75))
  pred <- predict(model, newdata, type="response", se.fit=T)
  
  out[[i]] <- get_diff_ci(pred)
}
out <- dplyr::bind_rows(out)

true_diff <- mean(out$diff)
out$true_in_ci <- true_diff <= out$upper & true_diff >= out$lower
mean(out$true_in_ci)

Here I simply repeated your simulation 10000 times and calculated the difference and confidence interval of the two predicted probabilities for each repetition. The mean of the individual differences of the probabilities can be used as an estimate for the true underlying difference (because I am too lazy to derive the actual value). By simply checking what proportion of the estimated CIs contain this true value we can judge whether the confidence interval is correct. In this case the confidence intervals contain the true value 94.59% of the time, which is almost exactly equal to the desired 95%. The discrepancy is probably due to the rounding of the 1.96 z-value and simulation error.
EDIT (15.06.2022):
As pointed out by @Sextus Empiricus, this method does not seem to work as well as it seemy by just looking at the single monte-carlo study and should probably not be used.
