Bagaimana cara mengatasi pemisahan yang sempurna dalam regresi logistik?

163

Jika Anda memiliki variabel yang dengan sempurna memisahkan nol dan yang ada di variabel target, R akan menghasilkan pesan peringatan "pemisahan sempurna atau kuasi sempurna":

Warning message:
glm.fit: fitted probabilities numerically 0 or 1 occurred 

Kami masih mendapatkan model tetapi estimasi koefisien meningkat.

Bagaimana Anda menangani hal ini dalam praktik?

pengguna333
sumber
4
pertanyaan
user603
1
pertanyaan dan demo terkait dengan regularisasi di sini
Haitao Du

Jawaban:

100

Solusi untuk ini adalah dengan menggunakan bentuk regresi yang dihukum. Bahkan, ini adalah alasan asli beberapa bentuk regresi yang dihukum dikembangkan (meskipun mereka ternyata memiliki sifat menarik lainnya.

Instal dan muat paket glmnet dalam R dan Anda siap untuk pergi. Salah satu aspek glmnet yang kurang ramah pengguna adalah Anda hanya bisa memberinya matriks, bukan formula seperti yang biasa kita lakukan. Namun, Anda dapat melihat model.matrix dan sejenisnya untuk membuat matriks ini dari data.frame dan rumus ...

Sekarang, ketika Anda berharap bahwa pemisahan sempurna ini bukan hanya produk sampingan dari sampel Anda, tetapi bisa benar dalam populasi, Anda secara khusus tidak ingin menangani ini: gunakan variabel pemisah ini hanya sebagai satu-satunya prediktor untuk hasil Anda, bukan menggunakan model apa pun.

Nick Sabbe
sumber
20
Anda juga dapat menggunakan antarmuka rumus untuk glmnet melalui paket caret.
Zach
"Sekarang, ketika kamu mengharapkan ..." Pertanyaan tentang ini. Saya memiliki studi kasus / kontrol yang melihat hubungan dengan microbiome. Kami juga memiliki perawatan yang hampir hanya ditemukan di antara kasus-kasus. Namun, kami pikir perawatannya mungkin juga mempengaruhi microbiome. Apakah ini contoh peringatan Anda? Secara hipotesis kami dapat menemukan lebih banyak kasus yang tidak menggunakan perawatan jika kami mencoba, tetapi kami memiliki apa yang kami miliki.
abalter
142

Anda memiliki beberapa opsi:

  1. Hapus beberapa bias.

    (A) Dengan menghukum kemungkinan sesuai saran @ Nick. Paket logistf dalam R atau FIRTHopsi dalam SAS PROC LOGISTICmengimplementasikan metode yang diusulkan dalam Firth (1993), "Bias pengurangan estimasi kemungkinan maksimum", Biometrika , 80 , 1; yang menghilangkan bias tingkat pertama dari perkiraan kemungkinan maksimum. ( Di sini @Gavin merekomendasikan brglmpaket, yang saya tidak kenal, tapi saya kumpulkan mengimplementasikan pendekatan yang serupa untuk fungsi tautan non-kanonik misalnya probit.)

    (B) Dengan menggunakan estimasi median-tidak bias dalam regresi logistik bersyarat yang tepat. Paket elrm atau logistiX dalam R, atau EXACTpernyataan dalam SAS PROC LOGISTIC.

  2. Kecualikan kasus di mana kategori atau nilai prediktor yang menyebabkan pemisahan terjadi. Ini mungkin berada di luar jangkauan Anda; atau layak untuk diselidiki lebih lanjut, terfokus. (Paket R safeBinaryRegression berguna untuk menemukannya.)

  3. Pasang kembali model. Biasanya ini adalah sesuatu yang akan Anda lakukan sebelumnya jika Anda memikirkannya, karena terlalu rumit untuk ukuran sampel Anda.

    (a) Hapus prediktor dari model. Dicey, dengan alasan yang diberikan oleh @Simon: "Anda menghapus prediktor yang paling menjelaskan respons".

    (B) Dengan menciutkan kategori prediktor / binning nilai-nilai prediktor. Hanya jika ini masuk akal.

    (c) Mengekspresikan kembali prediktor sebagai dua (atau lebih) faktor silang tanpa interaksi. Hanya jika ini masuk akal.

  4. 5212

  5. Tidak melakukan apapun. (Tapi hitung interval kepercayaan berdasarkan kemungkinan profil, karena estimasi Wald tentang kesalahan standar akan sangat salah.) Pilihan yang sering diabaikan. Jika tujuan dari model ini hanya untuk menggambarkan apa yang telah Anda pelajari tentang hubungan antara prediktor & respons, tidak ada salahnya mengutip interval kepercayaan untuk rasio odds, katakanlah, 2,3 ke atas. (Memang bisa terlihat mencurigakan mengutip interval kepercayaan berdasarkan estimasi yang tidak bias yang mengecualikan rasio odds yang paling didukung oleh data.) Masalah muncul ketika Anda mencoba untuk memprediksi menggunakan estimasi titik, & prediksi di mana pemisahan terjadi membanjiri yang lain.

  6. Gunakan model regresi logistik tersembunyi, seperti yang dijelaskan dalam Rousseeuw & Christmann (2003), "Ketangguhan terhadap pemisahan dan pencilan dalam regresi logistik", Statistik Komputasi & Analisis Data , 43 , 3, dan diimplementasikan dalam paket R hlr . (@ user603 menyarankan ini. ) Saya belum membaca makalah, tetapi mereka mengatakan secara abstrak "model yang sedikit lebih umum diusulkan di mana respons yang diamati sangat terkait tetapi tidak sama dengan respons sejati yang tidak dapat diamati", yang menyarankan untuk saya mungkin bukan ide yang baik untuk menggunakan metode ini kecuali itu terdengar masuk akal.

  7. "Ubah beberapa pengamatan yang dipilih secara acak dari 1 menjadi 0 atau 0 ke 1 di antara variabel yang menunjukkan pemisahan total": @ komentar RobertF . Saran ini tampaknya muncul dari menganggap pemisahan sebagai masalah per se daripada sebagai gejala kurangnya informasi dalam data yang mungkin membuat Anda lebih memilih metode lain untuk estimasi kemungkinan maksimum, atau untuk membatasi kesimpulan kepada mereka yang dapat Anda buat dengan presisi yang masuk akal — pendekatan yang memiliki kelebihannya sendiri & bukan sekadar "perbaikan" untuk pemisahan. (Selain dari ad hoc tanpa malu-malu , itu tidak menyenangkan bagi sebagian besar analis yang menanyakan pertanyaan yang sama dari data yang sama, membuat asumsi yang sama, harus memberikan jawaban yang berbeda karena hasil lemparan koin atau apa pun.)

Scortchi
sumber
1
@Scortchi Ada opsi (sesat) lainnya. Bagaimana dengan mengubah beberapa pengamatan yang dipilih secara acak dari 1 menjadi 0 atau 0 menjadi 1 di antara variabel yang menunjukkan pemisahan total?
RobertF
@RobertF: Terima kasih! Saya belum memikirkan yang ini - jika Anda memiliki referensi mengenai kinerjanya, saya akan berterima kasih. Sudahkah Anda menemukan orang yang menggunakannya dalam praktik?
Scortchi
@ Scortchi - Tidak, ada referensi untuk peneliti menambahkan data buatan untuk menghilangkan pemisahan lengkap, tapi saya belum menemukan artikel tentang modifikasi selektif data. Saya tidak tahu seberapa efektif metode ini.
RobertF
1
@tatami: Tidak semua (banyak?) program memperingatkan tentang pemisahan per se, yang bisa sulit dikenali ketika itu pada kombinasi linear dari beberapa variabel, tetapi tentang kegagalan konvergensi & / atau nilai-nilai yang dipasang dekat dengan nol atau satu - saya akan selalu periksa ini.
Scortchi
2
@Scortchi: ringkasan yang sangat bagus dalam jawaban Anda. Secara pribadi saya menyukai pendekatan Bayesian tetapi perlu menyebutkan analisis indah dari fenomena umum dari sudut pandang yang sering di projecteuclid.org/euclid.ejs/1239716414 . Penulis menawarkan beberapa interval kepercayaan satu sisi yang dapat digunakan bahkan di hadapan pemisahan lengkap dalam regresi logistik.
Cyan
55

Ini adalah perluasan jawaban Scortchi dan Manoel, tetapi karena Anda tampaknya menggunakan RI, saya pikir saya akan menyediakan beberapa kode. :)

Saya percaya solusi termudah dan paling langsung untuk masalah Anda adalah dengan menggunakan analisis Bayesian dengan asumsi sebelumnya yang tidak informatif seperti yang diusulkan oleh Gelman et al (2008). Seperti Scortchi menyebutkan, Gelman merekomendasikan untuk menempatkan Cauchy sebelumnya dengan median 0,0 dan skala 2,5 pada masing-masing koefisien (dinormalisasi untuk memiliki rata-rata 0,0 dan SD 0,5). Ini akan mengatur koefisien dan menariknya sedikit ke nol. Dalam hal ini persis apa yang Anda inginkan. Karena memiliki ekor yang sangat lebar, Cauchy masih memungkinkan untuk koefisien besar (sebagai lawan dari Normal berekor pendek), dari Gelman:

masukkan deskripsi gambar di sini

Bagaimana cara menjalankan analisis ini? Gunakan bayesglmfungsi dalam paket lengan yang mengimplementasikan analisis ini!

library(arm)

set.seed(123456)
# Faking some data where x1 is unrelated to y
# while x2 perfectly separates y.
d <- data.frame(y  =  c(0,0,0,0, 0, 1,1,1,1,1),
                x1 = rnorm(10),
                x2 = sort(rnorm(10)))

fit <- glm(y ~ x1 + x2, data=d, family="binomial")

## Warning message:
## glm.fit: fitted probabilities numerically 0 or 1 occurred 

summary(fit)
## Call:
## glm(formula = y ~ x1 + x2, family = "binomial", data = d)
##
## Deviance Residuals: 
##       Min          1Q      Median          3Q         Max  
## -1.114e-05  -2.110e-08   0.000e+00   2.110e-08   1.325e-05  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept)    -18.528  75938.934       0        1
## x1              -4.837  76469.100       0        1
## x2              81.689 165617.221       0        1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1.3863e+01  on 9  degrees of freedom
## Residual deviance: 3.3646e-10  on 7  degrees of freedom
## AIC: 6
## 
## Number of Fisher Scoring iterations: 25

Tidak berfungsi dengan baik ... Sekarang versi Bayesian:

fit <- bayesglm(y ~ x1 + x2, data=d, family="binomial")
display(fit)
## bayesglm(formula = y ~ x1 + x2, family = "binomial", data = d)
##             coef.est coef.se
## (Intercept) -1.10     1.37  
## x1          -0.05     0.79  
## x2           3.75     1.85  
## ---
## n = 10, k = 3
## residual deviance = 2.2, null deviance = 3.3 (difference = 1.1)

Sangat sederhana, bukan?

Referensi

Gelman et al (2008), "Sebuah distribusi sebelumnya standar informatif lemah untuk logistik & model regresi lainnya", Ann. Appl. Stat., 2, 4 http://projecteuclid.org/euclid.aoas/1231424214

Rasmus Bååth
sumber
6
Tidak. Terlalu sederhana. Bisakah Anda menjelaskan apa yang baru saja Anda lakukan? Apa yang bayesglmdigunakan sebelumnya ? Jika estimasi ML setara dengan Bayesian dengan flat sebelumnya, bagaimana prior non-informatif membantu di sini?
Tugas
5
Menambahkan beberapa info lagi! Sebelumnya tidak jelas tetapi tidak datar. Ini memiliki beberapa pengaruh karena mengatur perkiraan dan menariknya sedikit ke arah 0,0 yang saya yakin Anda inginkan dalam kasus ini.
Rasmus Bååth
> m = bayesglm (kecocokan ~., keluarga = binomial (tautan = 'logit'), data = df) Pesan peringatan: probabilitas pas dipasang secara numerik 0 atau 1 terjadi Tidak baik!
Chris
Sebagai pemula, cobalah regularisasi yang sedikit lebih kuat dengan meningkatkan prior.dfdefault ke 1.0dan / atau mengurangi prior.scaledefault ke 2.5, mungkin mulai mencoba:m=bayesglm(match ~. , family = binomial(link = 'logit'), data = df, prior.df=5)
Rasmus Bååth
1
Apa sebenarnya yang kita lakukan ketika kita meningkatkan prior.df dalam model. Apakah ada batasan seberapa tinggi kita ingin pergi? Pemahaman saya adalah bahwa itu membatasi model untuk memungkinkan konvergensi dengan estimasi kesalahan yang akurat?
hamilthj
7

Salah satu penjelasan paling menyeluruh tentang masalah "pemisahan semu" dalam kemungkinan maksimum adalah makalah Paul Allison. Dia menulis tentang perangkat lunak SAS, tetapi masalah yang dia tangani dapat digeneralisasi untuk perangkat lunak apa pun:

  • Pemisahan sempurna terjadi setiap kali fungsi linear x dapat menghasilkan prediksi y yang sempurna

  • Pemisahan semu -lengkap terjadi ketika (a) ada beberapa koefisien vektor b sedemikian sehingga bxi ≥ 0 setiap kali yi = 1 , dan bxi ≤ 0 * setiap kali ** yi = 0 dan persamaan ini berlaku untuk setidaknya satu kasus di setiap kategori variabel tak bebas. Dengan kata lain dalam kasus paling sederhana, untuk setiap variabel independen dikotomis dalam regresi logistik, jika ada nol dalam tabel 2 × 2 yang dibentuk oleh variabel itu dan variabel dependen, estimasi ML untuk koefisien regresi tidak ada.

Allison membahas banyak solusi yang telah disebutkan termasuk penghapusan variabel masalah, kategori runtuh, tidak melakukan apa-apa, meningkatkan regresi logistik yang tepat , estimasi Bayesian dan estimasi kemungkinan maksimum hukuman.

http://www2.sas.com/proceedings/forum2008/360-2008.pdf

Mike Hunter
sumber
3

warning

Dengan data yang dihasilkan di sepanjang baris

x <- seq(-3, 3, by=0.1)
y <- x > 0
summary(glm(y ~ x, family=binomial))

Peringatan itu dibuat:

Warning messages:
1: glm.fit: algorithm did not converge 
2: glm.fit: fitted probabilities numerically 0 or 1 occurred 

yang sangat jelas mencerminkan ketergantungan yang dibangun ke dalam data ini.

Di R tes Wald ditemukan dengan summary.glmatau dengan waldtestdalam lmtestpaket. Tes rasio kemungkinan dilakukan dengan anovaatau dengan lrtestdalam lmtestpaket. Dalam kedua kasus tersebut, matriks informasi dihargai secara tidak terbatas, dan tidak ada kesimpulan yang tersedia. Sebaliknya, R memang menghasilkan output, tetapi Anda tidak bisa mempercayainya. Kesimpulan yang dihasilkan R dalam kasus-kasus ini memiliki nilai-p sangat dekat dengan satu. Ini karena kehilangan presisi dalam OR adalah urutan besarnya lebih kecil daripada hilangnya presisi dalam matriks varians-kovarians.

Beberapa solusi yang diuraikan di sini:

Gunakan penduga satu langkah,

Ada banyak teori yang mendukung bias, efisiensi, dan kemampuan generalisasi satu penduga satu langkah yang rendah. Mudah untuk menentukan penduga satu langkah dalam R dan hasilnya biasanya sangat menguntungkan untuk prediksi dan inferensi. Dan model ini tidak akan pernah menyimpang, karena iterator (Newton-Raphson) tidak memiliki kesempatan untuk melakukannya!

fit.1s <- glm(y ~ x, family=binomial, control=glm.control(maxit=1))
summary(fit.1s)

Memberi:

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.03987    0.29569  -0.135    0.893    
x            1.19604    0.16794   7.122 1.07e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1

Sehingga Anda bisa melihat prediksi yang mencerminkan arah tren. Dan kesimpulannya sangat menunjukkan kecenderungan yang kami yakini benar.

masukkan deskripsi gambar di sini

melakukan tes skor,

The Score (atau Rao) statistik berbeda dari rasio kemungkinan dan wald statistik. Itu tidak memerlukan evaluasi varian di bawah hipotesis alternatif. Kami cocok dengan model di bawah nol:

mm <- model.matrix( ~ x)
fit0 <- glm(y ~ 1, family=binomial)
pred0 <- predict(fit0, type='response')
inf.null <- t(mm) %*% diag(binomial()$variance(mu=pred0)) %*% mm
sc.null <- t(mm) %*% c(y - pred0)
score.stat <- t(sc.null) %*% solve(inf.null) %*% sc.null ## compare to chisq
pchisq(score.stat, 1, lower.tail=F)

χ2

> pchisq(scstat, df=1, lower.tail=F)
             [,1]
[1,] 1.343494e-11

Dalam kedua kasus Anda memiliki inferensi untuk OR tanpa batas.

, dan gunakan estimasi median yang tidak bias untuk interval kepercayaan.

Anda dapat menghasilkan median rata-rata, CI 95% non-singular untuk rasio odds tak terbatas dengan menggunakan estimasi median bias. Paket epitoolsdi R dapat melakukan ini. Dan saya memberikan contoh penerapan penaksir ini di sini: Interval kepercayaan untuk pengambilan sampel Bernoulli

AdamO
sumber
2
Ini bagus, tetapi saya punya beberapa quibbles, tentu saja: (1) Tes kemungkinan-rasio tidak menggunakan matriks informasi; hanya tes Wald yang berhasil, & yang gagal secara dahsyat di hadapan perpisahan. (2) Saya sama sekali tidak terbiasa dengan penduga satu langkah, tetapi perkiraan kemiringan di sini tampaknya sangat tidak masuk akal. (3) Interval kepercayaan tidak rata-rata. Apa yang Anda tautkan pada bagian itu adalah interval kepercayaan pertengahan-p. (4) Anda dapat memperoleh interval kepercayaan dengan membalikkan LR atau skor tes. ...
Scortchi
... (5) Anda dapat melakukan tes skor dalam R dengan memberikan argumen test="Rao"ke anovafungsi. (Ya, dua yang terakhir adalah not, bukan quibbles.)
Scortchi
@scortchi senang mengetahui anova memiliki tes skor default! Mungkin implementasi dengan tangan berguna. CIs tidak median bias, tetapi CIs untuk estimator rata-rata median memberikan inferensi yang konsisten untuk parameter batas. Mid p adalah penaksir seperti itu. P dapat ditransformasikan ke rasio odds b / c itu tidak berubah untuk transformasi satu-ke-satu. Apakah tes LR konsisten untuk parameter batas?
AdamO
Hanya hipotesis nol yang tidak boleh berisi parameter pada batas yang berlaku untuk teorema Wilks, meskipun skor & LR diperkirakan dalam sampel terbatas.
Scortchi
2

Hati-hati dengan pesan peringatan ini dari R. Lihatlah posting blog ini oleh Andrew Gelman, dan Anda akan melihat bahwa itu tidak selalu merupakan masalah perpisahan yang sempurna, tetapi kadang-kadang bug dengan glm. Tampaknya jika nilai awal terlalu jauh dari perkiraan kemungkinan maksimum, nilai itu akan meledak. Jadi, periksa dulu dengan perangkat lunak lain, seperti Stata.

Jika Anda benar-benar memiliki masalah ini, Anda dapat mencoba menggunakan pemodelan Bayesian, dengan prior informatif.

Tetapi dalam praktiknya saya hanya menyingkirkan prediktor yang menyebabkan masalah, karena saya tidak tahu bagaimana memilih yang informatif sebelumnya. Tapi saya kira ada makalah dari Gelman tentang menggunakan informasi sebelumnya ketika Anda memiliki masalah dengan masalah pemisahan yang sempurna ini. Hanya google saja. Mungkin Anda harus mencobanya.

Manoel Galdino
sumber
8
Masalah dengan menghapus prediktor adalah bahwa Anda menghapus prediktor yang paling baik menjelaskan responsnya, yang biasanya adalah apa yang ingin Anda lakukan! Saya berpendapat bahwa ini hanya masuk akal jika Anda sudah melebihi model Anda, misalnya dengan memasang terlalu banyak interaksi yang rumit.
Simon Byrne
4
Bukan bug, tetapi masalah dengan perkiraan awal yang terlalu jauh dari MLE, yang tidak akan muncul jika Anda tidak mencoba memilihnya sendiri.
Scortchi
Saya mengerti ini, tapi saya pikir ini adalah Bug dalam algoritme.
Manoel Galdino
5
Yah saya tidak ingin berdalih tentang definisi 'bug'. Tetapi perilaku itu tidak terduga dan tidak dapat diperbaiki dalam basis R - Anda tidak perlu "memeriksa dengan perangkat lunak lain". Jika Anda ingin berurusan secara otomatis dengan banyak masalah non-konvergensi, glm2paket mengimplementasikan pemeriksaan bahwa kemungkinan sebenarnya meningkat pada setiap langkah pemberian skor, & membagi dua ukuran langkah jika tidak.
Scortchi
3
Ada (pada CRAN) paket R safeBinaryRegression yang dirancang untuk mendiagnosis dan memperbaiki masalah seperti itu, menggunakan metode optimasi untuk memeriksa apakah ada pemisahan atau quasiseparation. Cobalah!
kjetil b halvorsen
2

Saya tidak yakin bahwa saya setuju dengan pernyataan dalam pertanyaan Anda.

Saya pikir pesan peringatan berarti, untuk beberapa tingkat X yang diamati dalam data Anda, probabilitas yang dipasang adalah angka 0 atau 1. Dengan kata lain, pada resolusi, itu ditampilkan sebagai 0 atau 1.

Anda dapat menjalankan predict(yourmodel,yourdata,type='response')dan Anda akan menemukan 0's atau / dan 1's di sana sebagai probabilitas diprediksi.

Sebagai hasilnya, saya pikir tidak apa-apa untuk hanya menggunakan hasilnya.

StayLearning
sumber
-1

Saya mengerti ini adalah posting lama, namun saya masih akan melanjutkan dengan menjawab ini karena saya telah berjuang berhari-hari dengan itu dan dapat membantu orang lain.

Pemisahan total terjadi ketika variabel yang Anda pilih agar sesuai dengan model dapat dengan sangat akurat membedakan antara 0 dan 1 atau ya dan tidak. Seluruh pendekatan kami dalam ilmu data didasarkan pada estimasi probabilitas tetapi gagal dalam kasus ini.

Langkah-langkah perbaikan: -

  1. Gunakan bayesglm () alih-alih glm (), ketika dalam kasus varians antara variabel rendah

  2. Pada saat menggunakan (maksit = ”beberapa nilai numerik”) bersama dengan bayesglm () dapat membantu

3. Periksa ketiga dan paling penting untuk variabel yang Anda pilih untuk pemasangan model, harus ada variabel yang multi collinearity dengan variabel Y (outout) sangat tinggi, buang variabel itu dari model Anda.

Seperti dalam kasus saya, saya memiliki data churn telekomunikasi untuk memprediksi churn untuk data validasi. Saya memiliki variabel dalam data pelatihan saya yang bisa sangat membedakan antara ya dan tidak. Setelah menjatuhkannya saya bisa mendapatkan model yang benar. Selanjutnya Anda dapat menggunakan stepwise (fit) untuk membuat model Anda lebih akurat.

yash
sumber
2
Saya tidak melihat bahwa jawaban ini menambah banyak diskusi. Pendekatan Bayesian sepenuhnya tercakup dalam jawaban sebelumnya, penghapusan prediktor "bermasalah" juga telah disebutkan (dan tidak disarankan). Seleksi stepwise variabel jarang merupakan ide bagus sejauh yang saya tahu.
einar