# Why is there a difference between manually calculating a logistic regression 95% confidence interval, and using the confint() function in R?

Dear everyone - I've noticed something strange that I can't explain, can you? In summary: the manual approach to calculating a confidence interval in a logistic regression model, and the R function confint() give different results.

I've been going through Hosmer & Lemeshow's Applied logistic regression (2nd edition). In the 3rd chapter there is an example of calculating the odds ratio and 95% confidence interval. Using R, I can easily reproduce the model:

Call:
glm(formula = dataset$CHD ~ as.factor(dataset$dich.age), family = "binomial")

Deviance Residuals:
Min      1Q  Median      3Q     Max
-1.734  -0.847  -0.847   0.709   1.549

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)                   -0.8408     0.2551  -3.296  0.00098 ***
as.factor(dataset$dich.age)1 2.0935 0.5285 3.961 7.46e-05 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 136.66 on 99 degrees of freedom Residual deviance: 117.96 on 98 degrees of freedom AIC: 121.96 Number of Fisher Scoring iterations: 4  However, when I calculate the confidence intervals of the parameters, I get a different interval to the one given in the text: > exp(confint(model)) Waiting for profiling to be done... 2.5 % 97.5 % (Intercept) 0.2566283 0.7013384 as.factor(dataset$dich.age)1 3.0293727 24.7013080


Hosmer & Lemeshow suggest the following formula:

$$e^{[\hat\beta_1\pm z_{1-\alpha/2}\times\hat{\text{SE}}(\hat\beta_1)]}$$

and they calculate the confidence interval for as.factor(dataset$dich.age)1 to be (2.9, 22.9). This seems straightforward to do in R: # upper CI for beta exp(summary(model)$coefficients[2,1]+1.96*summary(model)$coefficients[2,2]) # lower CI for beta exp(summary(model)$coefficients[2,1]-1.96*summary(model)$coefficients[2,2])  gives the same answer as the book. However, any thoughts on why confint() seems to give different results? I've seen lots of examples of people using confint(). - add comment ## 3 Answers After having fetched the data from the accompagnying website, here is how I would do: chdage <- read.table("chdage.dat", header=F, col.names=c("id","age","chd")) chdage$aged <- ifelse(chdage\$age>=55, 1, 0)
mod.lr <- glm(chd ~ aged, data=chdage, family=binomial)
summary(mod.lr)


The 95% CIs based on profile likelihood is obtained with

require(MASS)
exp(confint(mod.lr))


This often is the default if package MASS is automatically loaded. In this case, I get

                2.5 %     97.5 %
(Intercept) 0.2566283  0.7013384
aged        3.0293727 24.7013080


Now, if I want to compare with 95% Wald CIs (based on asymptotic normality) like the one you computed by hand, I would use confint.default() instead; this yields:

                2.5 %     97.5 %
(Intercept) 0.2616579  0.7111663
aged        2.8795652 22.8614705


Wald CIs are good in most situation, although profile likelihood-based may be useful with complex sampling strategies. If you want to grasp the idea of how they work, here is a brief overview of the main principles: Confidence intervals by the profile likelihood method, with applications in veterinary epidemiology. You can also take a look at Venables and Ripley's MASS book, §8.4, pp. 220-221.

-
Thanks! very helpful! –  Andrew Dec 9 '10 at 21:43