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().

After having fetched the data from the accompanying website, here is how I would do it:

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 are obtained with

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


This often is the default if the MASS package 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 wanted 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 situations, 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.
