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

Question

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

Answer 1

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.

Answer 2

Following up: profile confidence intervals are more reliable (choosing the appropriate cutoff for the likelihood does involve an asymptotic (large sample) assumption, but this is a much weaker assumption than the quadratic-likelihood-surface assumption underlying the Wald confidence intervals). As far as I know, there is no argument for the Wald statistics over the profile confidence intervals except that the Wald statistics are much quicker to compute and may be "good enough" in many circumstances (but sometimes way off: look up the Hauck-Donner effect).

Answer 3

I believe if you look into the help file for confint() you will find that the confidence interval being constructed is a "profile" interval instead of a Wald confidence interval (your formula from HL).