Calculation confidence interval for estimated logit in logistic regression in R

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?

• What about confint? Jan 10 '12 at 3:18
• There's a se.fit option in predict(). A confidence intervals can be constructed around each fitted logit is given by taking the fitting logit value +- 1.96 * the corresponding estimated standard error. Alternatively, use matrix algebra, and construct confidence intervals for the fitted logits, by treating them as linear combinations of the estimated coefficients, returned by e.g. coef(chdage.lr.o2). (The former method is likely easier) Jan 10 '12 at 7:14
• @guest can you give us an example about se.fit? Thanks. Jan 10 '12 at 7:17

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'))
• @Harrell: Why the option of simultaneous prediction intervals is not available in the most recent version at CRAN? for the line: p3 <- Predict(fit, age=c(60,65,70), sex='male', fun=plogis,conf.type='simultaneous') I got the following error: Error in match.arg(conf.type) : 'arg' should be one of “mean”, “individual” Jan 18 '12 at 1:15
• Please read my note carefully. You apparently are not using Linux. See my note for Mac and Windows availability. Jan 18 '12 at 3:26
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)