# Standard error for the difference between two predicted probabilities

Imagine a logit model with a continuous independent variable, a binary independent variable (the treatment) and an interaction between the two. In a first step, I want to predict probabilities at extreme values of the continuous variable (i.e. at the two ends of the scale) for both treatment levels. Standard errors for these estimates can theoretically be derived with the delta method, as explained here, and implemented in R with the modmarg::marg command. Just considering one value of the continuous scale for now, let us call the standard errors $$\sigma_{t=0}$$ for the prediction of one treatment level, and $$\sigma_{t=1}$$ for the prediction at the other treatment level.

In a second step, I am interested in the treatment difference at the same values of the continuous independent variable. Estimation is straight-forward, but I am unsure about how to calculate the standard error of the difference. According to the variance sum law, if random variables $$X$$ and $$Y$$ are independent, then

$$Var(X \pm Y) = Var(X) + Var(Y).$$

So if the predicted probabilities per treatment were independent, the standard error of the difference $$\sigma_\Delta$$ would be:

$$\sqrt{\sigma_{t=0}^2+\sigma_{t=1}^2}$$

This is indeed what Stata's margins command seems (did not find an equivalent way in R to calculate difference plus SE directly). But is it justified to assume independence?