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,

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

Odpowiedzi: