Mengapa ada perbedaan antara menghitung interval kepercayaan 95% regresi logistik secara manual, dan menggunakan fungsi confint () di R?

34

Dear everyone - Saya telah memperhatikan sesuatu yang aneh yang tidak dapat saya jelaskan, bukan? Singkatnya: pendekatan manual untuk menghitung interval kepercayaan dalam model regresi logistik, dan fungsi R confint()memberikan hasil yang berbeda.

Saya telah melalui regresi logistik Terapan Hosmer & Lemeshow (edisi ke-2). Dalam bab ke-3 ada contoh menghitung rasio odds dan interval kepercayaan 95%. Menggunakan R, saya dapat dengan mudah mereproduksi 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

Namun, ketika saya menghitung interval kepercayaan parameter, saya mendapatkan interval berbeda dengan yang diberikan dalam teks:

> 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 menyarankan rumus berikut:

e^{[{\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)]}$

dan mereka menghitung interval kepercayaan untuk as.factor(dataset$dich.age)1menjadi (2.9, 22.9).

Ini tampaknya mudah dilakukan di 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])

memberikan jawaban yang sama seperti buku itu.

Namun, ada pemikiran mengapa confint()nampaknya memberikan hasil yang berbeda? Saya telah melihat banyak contoh orang menggunakan 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 Andrew
sumber

1

Maukah Anda menambahkan referensi literatur yang tepat untuk Hosmer & Lemeshow? Saya sudah lama mencari saran di kepompong dan buku-buku mereka, tetapi belum menemukannya.

DavidR

36

Setelah mengambil data dari situs web yang menyertainya , berikut adalah bagaimana saya akan melakukannya:

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)

95% CI berdasarkan kemungkinan profil diperoleh dengan

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

Ini sering merupakan default jika MASSpaket dimuat secara otomatis. Dalam hal ini, saya mengerti

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

Sekarang, jika saya ingin membandingkan dengan 95% Wald CIs (berdasarkan normalitas asimptotik) seperti yang Anda hitung dengan tangan, saya akan menggunakan confint.default()sebagai gantinya; ini menghasilkan

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

Wald CIs baik dalam kebanyakan situasi, meskipun profil berbasis kemungkinan dapat berguna dengan strategi pengambilan sampel yang kompleks. Jika Anda ingin memahami cara kerjanya, berikut adalah ikhtisar singkat tentang prinsip-prinsip utama: Interval kepercayaan dengan metode kemungkinan profil, dengan aplikasi dalam epidemiologi veteriner . Anda juga dapat melihat buku MASS Venables dan Ripley, §8.4, hlm. 220-221.

chl
sumber

25

Tindak lanjut: interval kepercayaan profil lebih dapat diandalkan (memilih cutoff yang sesuai untuk kemungkinan melibatkan asumsi asimptotik (sampel besar), tetapi ini adalah asumsi yang jauh lebih lemah daripada asumsi kuadratik-kemungkinan-permukaan yang mendasari interval kepercayaan Wald). Sejauh yang saya tahu, tidak ada argumen untuk statistik Wald atas interval kepercayaan profil kecuali bahwa statistik Wald lebih cepat untuk dihitung dan mungkin "cukup baik" dalam banyak keadaan (tapi kadang-kadang jauh: mencari Hauck- Efek donner).

Ben Bolker
sumber

2

Terima kasih untuk ini, dan untuk saran saya mencari efek Hauck-Donner. Efeknya tidak mendapatkan banyak perawatan di buku pelajaran tetapi tampaknya cukup penting!

Andrew

18

Saya percaya jika Anda melihat file bantuan untuk confint () Anda akan menemukan bahwa interval kepercayaan yang sedang dibangun adalah interval "profil" alih-alih interval kepercayaan Wald (formula Anda dari HL).

B_Miner
sumber

5

Ahh. Itu menjawab pertanyaan. Namun, itu mengarah ke yang berikutnya - mana yang lebih disukai?

Andrew

Mengapa ada perbedaan antara menghitung interval kepercayaan 95% regresi logistik secara manual, dan menggunakan fungsi confint () di R?

Jawaban: