Jak radzić sobie z idealną separacją w regresji logistycznej?

Jeśli masz zmienną, która doskonale oddziela zera i jedynki w zmiennej docelowej, R wyświetli następujący komunikat ostrzegawczy „idealna lub quasi idealna separacja”:

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

Nadal otrzymujemy model, ale szacunki współczynników są zawyżone.

Jak sobie z tym radzisz w praktyce?

Rozwiązaniem tego jest zastosowanie formy regresji karnej. W rzeczywistości jest to pierwotny powód, dla którego opracowano niektóre formy regresji karanej (choć okazały się mieć inne interesujące właściwości).

Zainstaluj i załaduj pakiet glmnet w R i jesteś w większości gotowy do pracy. Jednym z mniej przyjaznych dla użytkownika aspektów glmnet jest to, że możesz karmić go tylko matrycami, a nie formułami, do których jesteśmy przyzwyczajeni. Możesz jednak spojrzeć na model.matrix i tym podobne, aby skonstruować tę matrycę z data.frame i formuły ...

Teraz, gdy oczekujesz, że ten idealny rozdział nie jest tylko produktem ubocznym twojej próbki, ale może być prawdziwy w populacji, w szczególności nie chcesz sobie z tym poradzić: użyj tej zmiennej oddzielającej po prostu jako jedynego predyktora dla twojego wyniku, a nie stosowanie dowolnego modelu.

Masz kilka opcji:

1. Usuń niektóre uprzedzenia.

   (a) Przez karanie prawdopodobieństwa zgodnie z sugestią @ Nicka. Logistf pakietu w R lub FIRTHopcja w SAS PROC LOGISTICimplementują metodę zaproponowaną w Firth (1993), „Redukcja błędu szacunków maksymalnego prawdopodobieństwa”, Biometrika , 80 , 1 .; co usuwa tendencyjność pierwszego rzędu z szacunków maksymalnego prawdopodobieństwa. ( Tutaj @Gavin zaleca brglmpakiet, którego nie znam, ale wydaje mi się, że implementuje podobne podejście do niekanonicznych funkcji łącza, np. Probit.)

   (b) Poprzez zastosowanie medianowo-obiektywnych szacunków w dokładnej warunkowej regresji logistycznej. Pakiet elrm lub LogistiX w R lub EXACToświadczenie w SAS PROC LOGISTIC.

2. Wyklucz przypadki, w których występuje kategoria predykcyjna lub wartość powodująca separację. Mogą one również znajdować się poza twoim zakresem; lub warte dalszego, ukierunkowanego dochodzenia. (Pakiet R safeBinaryRegression jest przydatny do ich znalezienia.)

3. Ponownie rzuć model. Zazwyczaj jest to coś, co zrobiłbyś wcześniej, gdybyś o tym pomyślał, ponieważ jest to zbyt skomplikowane jak na twoją próbkę.

   (a) Usuń predyktor z modelu. Dicey, z powodów podanych przez @ Simon: „Usuwasz predyktor, który najlepiej wyjaśnia odpowiedź”.

   (b) Zwijając kategorie predyktorów / dzieląc wartości predyktorów. Tylko jeśli ma to sens.

   (c) Ponowne wyrażanie predyktora jako dwóch (lub więcej) skrzyżowanych czynników bez interakcji. Tylko jeśli ma to sens.

4. $\frac{5}{2}$$\frac{1}{2}$

5. Nic nie robić. (Ale oblicz przedziały ufności w oparciu o prawdopodobieństwa profilu, ponieważ szacunki Walda dotyczące błędu standardowego będą bardzo błędne.) Często przeoczona opcja. Jeśli celem tego modelu jest po prostu opisanie tego, czego dowiedziałeś się o związkach między predyktorami i reakcją, nie ma wstydu w cytowaniu przedziału ufności dla ilorazu szans, powiedzmy, 2,3 w górę. (Rzeczywiście, może wydawać się podejrzane cytowanie przedziałów ufności w oparciu o obiektywne szacunki, które wykluczają iloraz szans najlepiej poparty danymi.) Problemy pojawiają się, gdy próbujesz przewidzieć przy użyciu oszacowań punktowych, a predyktor wystąpienia separacji zalewa inne.

6. Użyj ukrytego modelu regresji logistycznej, jak opisano w Rousseeuw i Christmann (2003), „Odporność na separację i wartości odstające w regresji logistycznej”, Statystyka obliczeniowa i analiza danych , 43 , 3, i zaimplementowane w pakiecie R hlr . (@ user603 to sugeruje. ) Nie czytałem tego artykułu, ale w streszczeniu mówią, że „zaproponowano nieco bardziej ogólny model, w którym obserwowana reakcja jest silnie powiązana, ale nie równa nieobserwowalnej prawdziwej odpowiedzi”, co sugeruje ja może nie być dobrym pomysłem użycie tej metody, chyba że brzmi to realistycznie.

7. „Zmiana kilka losowo wybranych uwag od 1 do 0 lub 0 do 1 spośród zmiennych wykazujących całkowite rozdzielenie”: @ robertf w komentarzu . Ta propozycja wydaje się wynikać z dotyczące separacji jako problem per se , a nie jako objaw niedostatku informacji w danych, które mogłyby doprowadzić do preferują inne metody szacowania największej wiarygodności, lub ograniczenia wnioskowania do tych można dokonać z rozsądna precyzja - podejścia, które mają swoje zalety i nie są jedynie „poprawkami” do separacji. (Poza tym, że jest to bezwstydnie ad hoc , dla większości analityków zadawanie tego samego pytania o te same dane, przyjmowanie tych samych założeń, powinno dawać różne odpowiedzi z powodu rzutu monetą lub cokolwiek innego).

Jest to rozszerzenie odpowiedzi Scortchi i Manoela, ale ponieważ wydaje się, że używasz RI, pomyślałem, że dostarczę trochę kodu. :)

Uważam, że najłatwiejszym i najprostszym rozwiązaniem problemu jest zastosowanie analizy bayesowskiej z wcześniejszymi nieinformacyjnymi założeniami, jak zaproponowali Gelman i wsp. (2008). Jak wspomina Scortchi, Gelman zaleca umieszczenie Cauchy'ego przed medianą 0,0 i skalą 2,5 na każdym współczynniku (znormalizowanym, aby mieć średnią 0,0 i SD 0,5). Spowoduje to wyrównanie współczynników i nieznaczne przesunięcie ich do zera. W tym przypadku jest dokładnie to, czego chcesz. Ze względu na bardzo szerokie ogony Cauchy nadal dopuszcza duże współczynniki (w przeciwieństwie do krótkiego ogona Normalnego), od Gelmana:

Jak uruchomić tę analizę? Użyj bayesglmfunkcji w pakiecie uzbrojenia, która implementuje tę analizę!

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

Nie działa tak dobrze ... Teraz wersja 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)

Super-proste, nie?

Bibliografia

Gelman i wsp. (2008), „Słabo informacyjna domyślna wcześniejsza dystrybucja dla modeli logistycznych i innych modeli regresji”, Ann. Appl. Stat., 2, 4 http://projecteuclid.org/euclid.aoas/1231424214

Jednym z najdokładniejszych wyjaśnień kwestii „quasi-całkowitej separacji” przy najwyższym prawdopodobieństwie jest praca Paula Allisona. Pisze o oprogramowaniu SAS, ale problemy, które rozwiązuje, można uogólnić na każde oprogramowanie:

• Całkowite rozdzielenie występuje, ilekroć funkcja liniowa x może wygenerować doskonałe przewidywania y

• Quasi-całkowite rozdzielenie występuje, gdy (a) istnieje jakiś wektor współczynnika b taki, że bxi ≥ 0 za każdym razem, gdy yi = 1 , i bxi ≤ 0 * za każdym razem ** yi = 0 i ta równość dotyczy co najmniej jednego przypadku w każdej kategorii zmienna zależna. Innymi słowy, w najprostszym przypadku, dla dowolnej dychotomicznej zmiennej niezależnej w regresji logistycznej, jeśli w tabeli 2 × 2 utworzonej przez tę zmienną i zmienną zależną jest zero, oszacowanie ML dla współczynnika regresji nie istnieje.

Allison omawia wiele już wspomnianych rozwiązań, w tym usuwanie zmiennych problemowych, zwijanie kategorii, nie robienie nic, wykorzystanie dokładnej regresji logistycznej, estymację Bayesa i karane oszacowanie maksymalnego prawdopodobieństwa.

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

warning$\infty$

Z danymi generowanymi zgodnie z

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

Ostrzeżenie jest wysyłane:

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

co bardzo wyraźnie odzwierciedla zależność wbudowaną w te dane.

W R test Walda znajduje się z summary.glmlub waldtestw lmtestpakiecie. Test współczynnika prawdopodobieństwa jest wykonywany z anovalub lrtestw lmtestpakiecie. W obu przypadkach matryca informacji jest nieskończenie ceniona i nie jest możliwe wnioskowanie. R tworzy raczej wynik, ale nie można mu ufać. Wnioskowanie, które R zwykle wytwarza w tych przypadkach, ma wartości p bardzo zbliżone do jednego. Jest tak, ponieważ utrata precyzji w OR jest o rząd wielkości mniejsza niż utrata precyzji w macierzy wariancji-kowariancji.

Niektóre rozwiązania przedstawione tutaj:

Użyj estymatora jednoetapowego,

Istnieje wiele teorii potwierdzających niski błąd systematyczny, wydajność i uogólnienie estymatorów jednoetapowych. Łatwo jest określić estymator jednoetapowy w R, a wyniki są zazwyczaj bardzo korzystne dla przewidywania i wnioskowania. I ten model nigdy się nie rozejdzie, ponieważ iterator (Newton-Raphson) po prostu nie ma na to szans!

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

Daje:

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

Dzięki temu prognozy odzwierciedlają kierunek trendu. A wnioskowanie bardzo sugeruje trendy, które uważamy za prawdziwe.

wykonać test punktowy,

Wynik (lub Rao) statystyczny różni się od stosunku prawdopodobieństwo i wald statystyki. Nie wymaga oceny wariancji zgodnie z alternatywną hipotezą. Dopasowujemy model do zera:

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)

$\chi^2$

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

W obu przypadkach wnioskujesz o OR nieskończoności.

i użyj mediany obiektywnych szacunków dla przedziału ufności.

Za pomocą mediany szacunku bezstronnego można wygenerować medianę obiektywnego 95% CI dla nieskończonego ilorazu szans. Pakiet epitoolsw R może to zrobić. Podam tutaj przykład implementacji tego estymatora: Przedział ufności dla próbkowania Bernoulliego

Bądź ostrożny z tym ostrzeżeniem od R. Spójrz na ten post na blogu Andrew Gelmana, a zobaczysz, że nie zawsze jest to problem idealnej separacji, ale czasem błąd glm. Wydaje się, że jeśli wartości początkowe są zbyt dalekie od oszacowania maksymalnego prawdopodobieństwa, wybuchnie. Sprawdź najpierw inne oprogramowanie, takie jak Stata.

Jeśli naprawdę masz ten problem, możesz spróbować użyć modelowania bayesowskiego z pouczającymi priorytetami.

Ale w praktyce po prostu pozbywam się predyktorów powodujących problemy, ponieważ nie wiem, jak wybrać pouczającego przeora. Ale wydaje mi się, że jest artykuł autorstwa Gelmana o używaniu informacyjnego przeora, kiedy masz problem z idealną separacją. Po prostu wyszukaj to w Google. Może powinieneś spróbować.

Nie jestem pewien, czy zgadzam się ze stwierdzeniami zawartymi w pytaniu.

Myślę, że ten komunikat ostrzegawczy oznacza, że ​​dla niektórych zaobserwowanych poziomów X w twoich danych, dopasowane prawdopodobieństwo wynosi 0 lub 1. Innymi słowy, przy rozdzielczości pokazuje 0 lub 1.

Możesz biec, predict(yourmodel,yourdata,type='response')a znajdziesz tam 0 i / lub 1 jako przewidywane prawdopodobieństwa.

W rezultacie myślę, że można po prostu użyć wyników.

Rozumiem, że to stary post, jednak nadal będę odpowiadać na to pytanie, ponieważ miałem z nim wiele dni i może to pomóc innym.

Całkowite rozdzielenie ma miejsce, gdy wybrane przez Ciebie zmienne pasujące do modelu mogą bardzo dokładnie rozróżniać między 0 a 1 lub tak i nie. Całe nasze podejście do analizy danych opiera się na oszacowaniu prawdopodobieństwa, ale w tym przypadku się nie powiedzie.

Kroki naprawy:

1. Użyj bayesglm () zamiast glm (), w przypadku gdy różnica między zmiennymi jest niska

2. Czasami pomocne może być użycie (maxit = „jakaś wartość liczbowa”) wraz z bayesglm ()

3. Trzecia i najważniejsza kontrola wybranych zmiennych dla dopasowania modelu, musi istnieć zmienna, dla której wieloliniowość ze zmienną Y (out) jest bardzo wysoka, odrzuć tę zmienną z modelu.

Podobnie jak w moim przypadku miałem dane o rezygnacji z usług telekomunikacyjnych, aby przewidzieć rezygnację z danych sprawdzania poprawności. Miałem zmienną w moich danych treningowych, która mogła bardzo odróżnić tak i nie. Po upuszczeniu mogłem uzyskać odpowiedni model. Co więcej, możesz zastosować stopniowanie (dopasowanie), aby zwiększyć dokładność modelu.