How are the standard errors computed for the fitted values from a logistic regression?

Question

When you predict a fitted value from a logistic regression model, how are standard errors computed? I mean for the fitted values, not for the coefficients (which involves Fishers information matrix).

I only found out how to get the numbers with R (e.g., here on r-help, or here on Stack Overflow), but I cannot find the formula.

pred <- predict(y.glm, newdata= something, se.fit=TRUE)

If you could provide online source (preferably on a university website), that would be fantastic.

Answer 1

The prediction is just a linear combination of the estimated coefficients. The coefficients are asymptotically normal so a linear combination of those coefficients will be asymptotically normal as well. So if we can obtain the covariance matrix for the parameter estimates we can obtain the standard error for a linear combination of those estimates easily. If I denote the covariance matrix as $\Sigma$ and and write the coefficients for my linear combination in a vector as $C$ then the standard error is just $\sqrt{C' \Sigma C}$

# Making fake data and fitting the model and getting a prediction
set.seed(500)
dat <- data.frame(x = runif(20), y = rbinom(20, 1, .5))
o <- glm(y ~ x, data = dat)
pred <- predict(o, newdata = data.frame(x=1.5), se.fit = TRUE)

# To obtain a prediction for x=1.5 I'm really
# asking for yhat = b0 + 1.5*b1 so my
# C = c(1, 1.5)
# and vcov applied to the glm object gives me
# the covariance matrix for the estimates
C <- c(1, 1.5)
std.er <- sqrt(t(C) %*% vcov(o) %*% C)

> pred$se.fit
[1] 0.4246289
> std.er
          [,1]
[1,] 0.4246289

We see that the 'by hand' method I show gives the same standard error as reported via predict