Dlaczego moje wartości p różnią się między wynikami regresji logistycznej, testem chi-kwadrat i przedziałem ufności dla OR?

Zbudowałem regresję logistyczną, w której zmienna wynikowa jest leczona po otrzymaniu leczenia ( Curevs. No Cure). Wszyscy pacjenci w tym badaniu zostali poddani leczeniu. Interesuje mnie, czy cukrzyca jest związana z tym wynikiem.

W R mój wynik regresji logistycznej wygląda następująco:

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

Jednak przedział ufności dla ilorazu szans obejmuje 1 :

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

Kiedy wykonuję test chi-kwadrat na tych danych, otrzymuję:

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

Jeśli chcesz to obliczyć samodzielnie, rozkład cukrzycy w wyleczonych i nieutwardzonych grupach przedstawia się następująco:

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

Moje pytanie brzmi: dlaczego wartości p i przedział ufności, w tym 1, nie są zgodne?

Dzięki uogólnionym modelom liniowym można uruchomić trzy różne typy testów statystycznych. Są to: testy Walda, testy współczynnika prawdopodobieństwa i testy punktowe. Doskonała strona pomocy na temat statystyk UCLA omawia je tutaj . Poniższy rysunek (skopiowany z ich witryny) pomaga je zilustrować:

1. Test Walda zakłada, że ​​prawdopodobieństwo jest zwykle rozkładane i na tej podstawie wykorzystuje stopień krzywizny do oszacowania błędu standardowego. Następnie oszacowanie parametru podzielone przez SE daje wynik- . Odnosi się pod dużym , ale nie jest to do końca prawda z mniejszą s. Trudno powiedzieć, kiedy twoje jest wystarczająco duże, aby utrzymać tę właściwość, więc ten test może być nieco ryzykowny. $z$$N$$N$$N$
2. Testy ilorazu wiarygodności sprawdzają iloraz prawdopodobieństw (lub różnic w prawdopodobieństwach logarytmicznych) na jego wartości maksymalnej i zerowej. Jest to często uważane za najlepszy test.
3. Test punktowy opiera się na nachyleniu prawdopodobieństwa przy wartości zerowej. Zwykle jest to mniej wydajne, ale zdarzają się sytuacje, w których nie można obliczyć pełnego prawdopodobieństwa, więc jest to dobra opcja rezerwowa.

Testy, które są dołączone, summary.glm()to testy Walda. Nie mówisz, jak uzyskałeś przedziały ufności, ale zakładam, że wykorzystałeś confint(), co z kolei wymaga profile(). Mówiąc dokładniej, te przedziały ufności są obliczane przez profilowanie prawdopodobieństwa (co jest lepszym podejściem niż pomnożenie SE przez ). Oznacza to, że są one analogiczne do testu współczynnika wiarygodności, a nie testu Walda. Z kolei test jest testem punktacji. $1.96$$\chi^2$

Gdy twoje staje się nieskończenie duże, trzy różne powinny zbiegać się na tej samej wartości, ale mogą się nieznacznie różnić, gdy nie masz nieskończonych danych. Warto zauważyć, że wartość (Wald) w początkowym wyniku jest ledwie znacząca i istnieje niewielka prawdziwa różnica między nieco powyżej i poniżej ( cytat ). Ta linia nie jest „magią”. Biorąc pod uwagę, że dwa bardziej wiarygodne testy są nieco ponad , powiedziałbym, że twoje dane nie są „znaczące” według konwencjonalnych kryteriów. $N$$p$$p$$\alpha=.05$.05$.05$

Poniżej profiluję współczynniki na skali predyktora liniowego i przeprowadzam test prawdopodobieństwa jawnie (via anova.glm()). Otrzymuję takie same wyniki jak Ty:

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

Jak zauważył @JWilliman w komentarzu (teraz usunięty), w Rmożna również uzyskać wartość p opartą na wynikach, używając anova.glm(model, test="Rao"). W poniższym przykładzie należy pamiętać, że wartość p nie jest całkiem takie same jak w teście chi-kwadrat powyżej, ponieważ domyślnie Rjest chisq.test()zastosowanie korekty ciągłości. Jeśli zmienimy to ustawienie, wartości p będą zgodne:

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