Calculation confidence interval for estimated logit in logistic regression in R

Question

I have read from here and understand how to calculate the estimated logit from a fitted logistic regression model, but how to work on the confidence interval? As it involved a variance-covariance matrix and I think it is better to have a program to do the calculation, rather then doing it by myself.

Thanks.

Edit 01

I have added a script here:

chdage.dummy <- data.frame(chd=c(rep(1,50),rep(0,50)),
                           race=c(rep("white",5),rep("black",20),rep("hispanic",15),rep("other",10),
                                  rep("white",20),rep("black",10),rep("hispanic",10),rep("other",10)),
                           stringsAsFactors=FALSE)
chdage.dummy[,"race"] <- factor(chdage.dummy[,"race"],levels=c("white","black","hispanic","other"))
chdage.lr.02 <- glm(chd~race,data=chdage.dummy,family="binomial")
predict(chdage.lr.02,newdata=data.frame(race="white"))

predict function can give me an estimate, but I can't use confint outside predict, so what can I do?

Answer 1

In version 3.4-0 of the R rms package (available now in CRAN for Linux and will probably be there by 19Jan12 for Mac and Windows) there are multiple ways to get confidence intervals. The following is from the help file for the plot.Predict function in rms.

fit <- lrm(y ~ blood.pressure + sex * (age + rcs(cholesterol,4)),
             x=TRUE, y=TRUE)

# For males at the median blood pressure and cholesterol, plot 3 types
# of confidence intervals for the probability on one plot, for varying age
ages <- seq(20, 80, length=100)
p1 <- Predict(fit, age=ages, sex='male', fun=plogis)  # standard pointwise
p2 <- Predict(fit, age=ages, sex='male', fun=plogis,
              conf.type='simultaneous')               # simultaneous
p3 <- Predict(fit, age=c(60,65,70), sex='male', fun=plogis,
              conf.type='simultaneous')               # simultaneous 3 pts
# The previous only adjusts for a multiplicity of 3 points instead of 100
f <- update(fit, x=TRUE, y=TRUE)
g <- bootcov(f, B=500, coef.reps=TRUE)
p4 <- Predict(g, age=ages, sex='male', fun=plogis)    # bootstrap percentile
p <- rbind(Pointwise=p1, 'Simultaneous 100 ages'=p2,
           'Simultaneous     3 ages'=p3, 'Bootstrap nonparametric'=p4)
xYplot(Cbind(yhat, lower, upper) ~ age, groups=.set.,
       data=p, type='l', method='bands', label.curve=list(keys='lines'))

Answer 2

new.data <- data.frame(race="white")
predict.fit.CI <- function(glmobj, newdata, level=0.05) {
  fit <- predict(glmobj, newdata=newdata, se.fit=TRUE)
  return(data.frame(fit=fit$fit, lower=fit$fit-(pnorm(1-(level/2))*fit$se.fit), upper=fit$fit+(pnorm(1-(level/2))*fit$se.fit)))
}
predict.fit.CI(chdage.lr.02, newdata)