Closed form for the variance of a sum of two estimates in logistic regression?

Question

In logistic regression with an intercept term and with at least one dependent variable which is categorical, is there a closed form for the variance of the sum of the intercept and the coefficient of the categorical variable, or do you have to sample from a multivariate distribution with the means and variances of the intercept and the coefficient to get a reliable measure of the variance of this sum?

$P(y) = \frac{exp(y)}{1 + exp(y)}$

$y = \beta_{0} + \beta_{1}x + \epsilon$

where x is categorical.

Would the formula for the variance of a sum (of two random variables) be applicable here?

$Var(\beta_{0} + \beta_{1}) = Var(\beta_{0}) + Var(\beta_{1}) + 2Cov(\beta_{0}, \beta_{1})$

The reason I ask, is that in this comment I got the impression that no closed form existed for the variance in question, and the advice was to sample from a multivariate distribution with the means and variances of $\beta_{0}$ and $\beta_{1}$

set.seed(1)
dependent.var <- sample(c(TRUE, FALSE), 100, replace = TRUE, prob = c(0.3, 0.7))
independent.var <- ifelse(dependent.var, sample(c("Red", "Blue"), replace = TRUE,
  size = 10, prob = c(0.8, 0.2)), sample(c("Red", "Blue"), size = 10,
  replace = TRUE, prob = c(0.4, 0.6)))
table(dependent.var, independent.var)
##               independent.var
## dependent.var Blue Red
##         FALSE   42  26
##         TRUE     7  25
my.fit <- glm(dependent.var ~ 1 + independent.var, family = binomial(logit))
coef(summary(my.fit))
                    Estimate Std. Error   z value     Pr(>|z|)
(Intercept)        -1.791759  0.4082482 -4.388897 0.0000113927
independent.varRed  1.752539  0.4951042  3.539737 0.0004005255
> vcov(my.fit)
                   (Intercept) independent.varRed
(Intercept)          0.1666666         -0.1666666
independent.varRed  -0.1666666          0.2451282

The logit of TRUE for a "Red" case is $ \beta_{0}+\beta_{1} \approx -0.039$. Is the variance for this estimate exactly

vcov(my.fit)[1,1] + vcov(my.fit)[2,2] + 2 * vcov(my.fit)[1,2]
[1] 0.07846154 ?

Or is this only an approximation, and a more accurate measure is to be found by sampling, e.g.

library(MASS)
var(rowSums(mvrnorm(n = 1E7, mu = coef(my.fit), Sigma = vcov(my.fit))))
[1] 0.07842985 ?

In this simple example, the sampling method does not seem to provide more accurate estimates of the variance (using 1E7 samples).

Here it is stated that "There is a correspondance between the covariance matrix of the fit parameters and Δχ2 confidence regions only for the case of Gaussian uncertainties on the input measurements.". Is that a reason against relying on the closed form above, or is there perhaps another reason for the advice to sample instead of deriving the variance analytically in cases like this?

EDIT: (In response to the answer given by StasK). The advice I originally got was to simulate from the full model, not from the vcov(), so here is the code to simulate from the full model:

library(arm)
sim.i <- sim(my.fit, 100000)
logit.for.TRUE.red <- sim.i@coef[,1] + sim.i@coef[,2]
var(logit.for.TRUE.red)
[1] 0.07781206

Answer 1

The standard error is indeed an approximation, but not in the sense that you are thinking about. It is an approximation in the statistical sense that the standard errors in logistic regression are only good in large samples (and biased down somewhat in small samples). There is a number of other ways to obtain the standard errors (e.g., by resampling methods such as the bootstrap or jackknife), but their justification is likewise asymptotic (large samples), and the standard errors obtained using these different methods should be close to one another.

Adding simulation on top of that estimate is relatively pointless. Bayesians out there might argue that you get better approximation for the distribution of the estimates if you sample from the posterior, but I don't think that's your question, and that's not what you are doing. But other than that, you won't in any way be better off simulating if all you use is vcov(). (And if you are simulating from the multivariate normal, the mean and the variance are independent, so if you are only interested in the variance, it does not matter what mean you use; since you are interested in variance, the vcov() is the only relevant part.)

Simulating from multivariate normal per se isn't interesting: nearly every distributional result of nearly any meaningful function of the multivariate normal has been derived in an analytic form. A potentially useful application of simulation could be to simulate the binary outcomes from your model (parametric bootstrap), fit the logistic model to them, and compare the mean of reported standard errors to the standard deviation of the estimates, to see how biased the latter are, and to correct for that bias. That is a somewhat paranoid simulation though; I don't think anybody does that in applied statistics publications that I tend to see, even though that is (a) technically feasible and (b) somewhat statistically meaningful. It's just too much work with rather little improvement in the standard errors on the back end.