Dlaczego istnieje różnica pomiędzy ręcznym obliczeniem regresji logistycznej 95% przedziału ufności a użyciem funkcji confint () w R?

34

Drodzy wszyscy - zauważyłem coś dziwnego, czego nie potrafię wyjaśnić, prawda? Podsumowując: ręczne podejście do obliczania przedziału ufności w modelu regresji logistycznej oraz funkcja R confint()dają różne wyniki.

Przechodziłem przez regresję logistyczną stosowaną przez Hosmer & Lemeshow (2. edycja). W trzecim rozdziale znajduje się przykład obliczenia ilorazu szans i 95% przedziału ufności. Za pomocą R mogę łatwo odtworzyć 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

Kiedy jednak obliczam przedziały ufności parametrów, otrzymuję inny przedział niż podany w tekście:

> 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 i Lemeshow sugerują następującą formułę:

{mi}^{[{\hat{β}}_{1} \pm z_{1 - α / 2)} \times \hat{SE} ({\hat{β}}_{1})]}

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

i obliczają przedział ufności dla as.factor(dataset$dich.age)1(2,9; 22,9).

Wydaje się to proste w 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])

daje taką samą odpowiedź jak książka.

Jednak jakieś przemyślenia na temat tego, dlaczego confint()wydaje się dawać różne wyniki? Widziałem wiele przykładów ludzi confint().

r regression logistic confidence-interval profile-likelihood correlation mcmc error mixture measurement data-augmentation r logistic goodness-of-fit r time-series exponential descriptive-statistics average expected-value data-visualization anova teaching hypothesis-testing multivariate-analysis r r mixed-model clustering categorical-data unsupervised-learning r logistic anova binomial estimation variance expected-value r r anova mixed-model multiple-comparisons repeated-measures project-management r poisson-distribution control-chart project-management regression residuals r distributions data-visualization r unbiased-estimator kurtosis expected-value regression spss meta-analysis r censoring regression classification data-mining mixture

— Andrzej
źródło

1

Czy miałbyś coś przeciwko dodaniu dokładnego odniesienia do literatury Hosmer & Lemeshow? Od dłuższego czasu szukałem sugestii w ich puplications i książkach, ale jeszcze jej nie znalazłem.

— DavidR

36

Po pobraniu danych z towarzyszącej strony internetowej , oto jak bym to zrobił:

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)

Wartości 95% CI oparte na prawdopodobieństwie profilu są uzyskiwane z

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

Jest to często ustawienie domyślne, jeśli MASSpakiet jest automatycznie ładowany. W tym przypadku rozumiem

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

Teraz, gdybym chciał porównać z 95% CI CI (opartymi na asymptotycznej normalności), takim jak ten, który obliczyłeś ręcznie, użyłbym confint.default()zamiast tego; to daje

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

Wartości CI Wald są dobre w większości sytuacji, chociaż oparte na prawdopodobieństwie profilu mogą być przydatne w przypadku złożonych strategii próbkowania. Jeśli chcesz zrozumieć, jak działają, oto krótkie omówienie głównych zasad: Przedziały ufności według metody prawdopodobieństwa profilu, z zastosowaniem w epidemiologii weterynaryjnej . Możesz także zajrzeć do książki MASY Venables i Ripleya, § 8.4, s. 220–221.

— chl
źródło

25

Kontynuacja: przedziały ufności profilu są bardziej wiarygodne (wybór odpowiedniego punktu odcięcia dla prawdopodobieństwa wiąże się z założeniem asymptotycznym (dla dużej próby), ale jest to założenie znacznie słabsze niż założenie kwadratowej powierzchni prawdopodobieństwa leżące u podstaw przedziałów ufności Walda). O ile mi wiadomo, nie ma argumentów za statystykami Walda w przedziałach ufności profilu, z wyjątkiem tego, że statystyki Walda są znacznie szybsze do obliczenia i mogą być „wystarczająco dobre” w wielu okolicznościach (ale czasami bardzo daleko: spójrz na Haucka- Efekt Donnera).

— Ben Bolker
źródło

2

Dzięki za to i za sugestię, żebym sprawdził efekt Haucka-Donnera. W podręcznikach efekt nie jest zbytnio traktowany, ale wydaje się bardzo ważny!

— Andrew,

18

Wierzę, że jeśli przejrzysz plik pomocy dla confint (), zobaczysz, że budowany przedział ufności jest przedziałem „profilowym” zamiast przedziału ufności Wald (twoja formuła z HL).

— B_Miner
źródło

5

Ahh To odpowiada na pytanie. Prowadzi to jednak do następnego - który jest preferowany?

— Andrew,