Why do my p-values differ between logistic regression output, chi-squared test, and the confidence interval for the OR?

Question

I have built a logistic regression where the outcome variable is being cured after receiving treatment (Cure vs. No Cure). All patients in this study received treatment. I am interested in seeing if having diabetes is associated with this outcome.

In R my logistic regression output looks as follows:

Call:
glm(formula = Cure ~ Diabetes, family = binomial(link = "logit"), data = All_patients)
...
Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   1.2735     0.1306   9.749   <2e-16 ***
Diabetes     -0.5597     0.2813  -1.990   0.0466 *  
...
    Null deviance: 456.55  on 415  degrees of freedom
Residual deviance: 452.75  on 414  degrees of freedom
  (2 observations deleted due to missingness)
AIC: 456.75

However, the confidence interval for the odds ratio includes 1:

                   OR     2.5 %   97.5 %
(Intercept) 3.5733333 2.7822031 4.646366
Diabetes    0.5713619 0.3316513 1.003167

When I do a chi-squared test on these data I get the following:

data:  check
X-squared = 3.4397, df = 1, p-value = 0.06365

If you'd like to calculate it on your own the distribution of diabetes in the cured and uncured groups are as follows:

Diabetic cure rate:      49 /  73 (67%)
Non-diabetic cure rate: 268 / 343 (78%)

My question is: Why don't the p-values and the confidence interval including 1 agree?

Answer 1

With generalized linear models, there are three different types of statistical tests that can be run. These are: Wald tests, likelihood ratio tests, and score tests. The excellent UCLA statistics help site has a discussion of them here. The following figure (copied from their site) helps to illustrate them:

1. The Wald test assumes that the likelihood is normally distributed, and on that basis, uses the degree of curvature to estimate the standard error. Then, the parameter estimate divided by the SE yields a $z$-score. This holds under large $N$, but isn't quite true with smaller $N$s. It is hard to say when your $N$ is large enough for this property to hold, so this test can be slightly risky.
2. Likelihood ratio tests look at the ratio of the likelihoods (or difference in log likelihoods) at its maximum and at the null. This is often considered the best test.
3. The score test is based on the slope of the likelihood at the null value. This is typically less powerful, but there are times when the full likelihood cannot be computed and so this is a nice fallback option.

The tests that come with summary.glm() are Wald tests. You don't say how you got your confidence intervals, but I assume you used confint(), which in turn calls profile(). More specifically, those confidence intervals are calculated by profiling the likelihood (which is a better approach than multiplying the SE by $1.96$). That is, they are analogous to the likelihood ratio test, not the Wald test. The $\chi^2$-test, in turn, is a score test.

As your $N$ becomes indefinitely large, the three different $p$'s should converge on the same value, but they can differ slightly when you don't have infinite data. It is worth noting that the (Wald) $p$-value in your initial output is just barely significant and there is little real difference between just over and just under $\alpha=.05$ (quote). That line isn't 'magic'. Given that the two more reliable tests are just over $.05$, I would say that your data are not quite 'significant' by conventional criteria.

Below I profile the coefficients on the scale of the linear predictor and run the likelihood ratio test explicitly (via anova.glm()). I get the same results as you:

library(MASS)
x = matrix(c(343-268,268,73-49,49), nrow=2, byrow=T);  x
#      [,1] [,2]
# [1,]   75  268
# [2,]   24   49
D = factor(c("N","Diabetes"), levels=c("N","Diabetes"))
m = glm(x~D, family=binomial)
summary(m)
# ...
# Coefficients:
#             Estimate Std. Error z value Pr(>|z|)    
# (Intercept)  -1.2735     0.1306  -9.749   <2e-16 ***
# DDiabetes     0.5597     0.2813   1.990   0.0466 *  
# ...
confint(m)
# Waiting for profiling to be done...
#                    2.5 %    97.5 %
# (Intercept) -1.536085360 -1.023243
# DDiabetes   -0.003161693  1.103671
anova(m, test="LRT")
# ...
#      Df Deviance Resid. Df Resid. Dev Pr(>Chi)  
# NULL                     1     3.7997           
# D     1   3.7997         0     0.0000  0.05126 .
chisq.test(x)
#         Pearson's Chi-squared test with Yates' continuity correction
# 
# X-squared = 3.4397, df = 1, p-value = 0.06365

---

As @JWilliman pointed out in a comment (now deleted), in R, you can also get a score-based p-value using anova.glm(model, test="Rao"). In the example below, note that the p-value isn't quite the same as in the chi-squared test above, because by default, R's chisq.test() applies a continuity correction. If we change that setting, the p-values match:

anova(m, test="Rao")
# ...
#      Df Deviance Resid. Df Resid. Dev   Rao Pr(>Chi)  
# NULL                     1     3.7997                 
# D     1   3.7997         0     0.0000 4.024  0.04486 *
chisq.test(x, correct=FALSE)
#   Pearson's Chi-squared test
# 
# data:  x
# X-squared = 4.024, df = 1, p-value = 0.04486