Dlaczego regresja logistyczna staje się niestabilna, gdy klasy są dobrze rozdzielone?

Dlaczego regresja logistyczna staje się niestabilna, gdy klasy są dobrze rozdzielone? Co oznaczają dobrze oddzielone klasy? Byłbym bardzo wdzięczny, gdyby ktoś mógł wyjaśnić na przykładzie.

Nieprawdą jest, że regresja logistyczna sama w sobie staje się niestabilna, gdy dochodzi do separacji. Separacja oznacza, że ​​istnieją pewne zmienne, które są bardzo dobrymi predyktorami, co jest dobre, lub separacja może być artefaktem zbyt małej liczby obserwacji / zbyt wielu zmiennych. W takim przypadku rozwiązaniem może być uzyskanie większej ilości danych. Ale sama separacja jest jedynie objawem, a nie problemem samym w sobie.

Są więc naprawdę różne przypadki do leczenia. Po pierwsze, jaki jest cel analizy? Jeśli końcowym wynikiem analizy jest jakaś klasyfikacja przypadków, separacja wcale nie stanowi problemu, to naprawdę oznacza, że ​​istnieją bardzo dobre zmienne dające bardzo dobrą klasyfikację. Ale jeśli celem jest oszacowanie ryzyka, potrzebujemy oszacowań parametrów, a przy rozdzieleniu zwykłe oszacowania mle (maksymalne prawdopodobieństwo) nie istnieją. Może więc musimy zmienić metodę szacowania. W literaturze jest kilka propozycji, wrócę do tego.

Następnie są (jak wspomniano powyżej) dwie różne możliwe przyczyny separacji. W pełnej populacji może wystąpić separacja lub separacja może wynikać z niewielu zaobserwowanych przypadków / zbyt wielu zmiennych.

Z separacją rozpada się procedura szacowania maksymalnego prawdopodobieństwa. Oszacowania parametru mle (lub przynajmniej niektóre z nich) stają się nieskończone. Powiedziałem w pierwszej wersji tej odpowiedzi, że można to łatwo rozwiązać, być może za pomocą ładowania początkowego, ale to nie działa, ponieważ w każdym ponownym próbkowaniu bootowania będzie separacja, przynajmniej przy zwykłej procedurze ładowania początkowego. Ale regresja logistyczna jest nadal prawidłowym modelem, ale potrzebujemy innej procedury szacowania. Niektóre propozycje to:

1. regularyzacja, np. kalenica lub lasso, może być połączona z bootstrapem.
2. dokładna warunkowa regresja logistyczna
3. testy permutacyjne, patrz https://www.ncbi.nlm.nih.gov/pubmed/15515134
4. Procedura estymacji zredukowanej uprzedzeniami Firthsa, patrz Model regresji logistycznej nie jest zbieżny
5. na pewno inni ...

Jeśli używasz R, na CRAN znajduje się pakiet SafeBinaryRegression, który pomaga w diagnozowaniu problemów z separacją, za pomocą matematycznych metod optymalizacji, aby sprawdzić, czy istnieje separacja czy quasi-separacja! Poniżej podam symulowany przykład z wykorzystaniem tego pakietu i elrmpakietu przybliżonej warunkowej regresji logistycznej.

Po pierwsze, prosty przykład z safeBinaryRegressionpakietem. Ten pakiet po prostu przedefiniowuje glmfunkcję, przeciążając ją testem separacji, stosując metody programowania liniowego. Jeśli wykryje separację, wychodzi z warunkiem błędu, deklarując, że mle nie istnieje. W przeciwnym razie po prostu uruchamia zwykłą glmfunkcję z stats. Przykład pochodzi ze stron pomocy:

library(safeBinaryRegression)   # Some testing of that package,
                                # based on its examples
# complete separation:
x  <-  c(-2, -1, 1, 2)
y  <-  c(0, 0, 1, 1)
glm(y ~ x, family=binomial)
glm(y ~ x,  family=binomial,  separation="test")
stats::glm(y~ x, family=binomial)
# Quasicomplete separation:
x  <-  c(-2, 0, 0, 2)
y  <-  c(0, 0, 1, 1)
glm(y ~ x, family=binomial)
glm(y ~ x,  family=binomial,  separation="test")
stats::glm(y~ x, family=binomial)

Dane wyjściowe z jego uruchomienia:

> # complete separation:
> x  <-  c(-2, -1, 1, 2)
> y  <-  c(0, 0, 1, 1)
> glm(y ~ x, family=binomial)
Error in glm(y ~ x, family = binomial) : 
  The following terms are causing separation among the sample points: (Intercept), x
> glm(y ~ x,  family=binomial,  separation="test")
Error in glm(y ~ x, family = binomial, separation = "test") : 
  Separation exists among the sample points.
    This model cannot be fit by maximum likelihood.
> stats::glm(y~ x, family=binomial)

Call:  stats::glm(formula = y ~ x, family = binomial)

Coefficients:
(Intercept)            x  
 -9.031e-08    2.314e+01  

Degrees of Freedom: 3 Total (i.e. Null);  2 Residual
Null Deviance:      5.545 
Residual Deviance: 3.567e-10    AIC: 4
Warning message:
glm.fit: fitted probabilities numerically 0 or 1 occurred 
> # Quasicomplete separation:
> x  <-  c(-2, 0, 0, 2)
> y  <-  c(0, 0, 1, 1)
> glm(y ~ x, family=binomial)
Error in glm(y ~ x, family = binomial) : 
  The following terms are causing separation among the sample points: x
> glm(y ~ x,  family=binomial,  separation="test")
Error in glm(y ~ x, family = binomial, separation = "test") : 
  Separation exists among the sample points.
    This model cannot be fit by maximum likelihood.
> stats::glm(y~ x, family=binomial)

Call:  stats::glm(formula = y ~ x, family = binomial)

Coefficients:
(Intercept)            x  
  5.009e-17    9.783e+00  

Degrees of Freedom: 3 Total (i.e. Null);  2 Residual
Null Deviance:      5.545 
Residual Deviance: 2.773    AIC: 6.773

Teraz symulujemy z modelu, który może być ściśle aproksymowany przez model logistyczny, z tym wyjątkiem, że powyżej pewnego poziomu odcięcia prawdopodobieństwo zdarzenia wynosi dokładnie 1,0. Pomyśl o problemie z testem biologicznym, ale ponad punktem odcięcia trucizna zawsze zabija:

pl  <-  function(a, b, x) 1/(1+exp(-a-b*x))
a  <-  0
b  <-  1.5
x_cutoff  <-  uniroot(function(x) pl(0,1.5,x)-0.98,lower=1,upper=3.5)$root
### circa 2.6
pltrue  <-  function(a, b, x) ifelse(x < x_cutoff, pl(a, b, x), 1.0)

x <- -3:3

### Let us simulate many times from this model,  and try to estimate it
### with safeBinaryRegression::glm  That way we can estimate the probability
### of separation from this model

set.seed(31415926)  ### May I have a large container of coffee 
replications  <-  1000
p  <-  pltrue(a, b, x)
err  <-  0
good  <- 0

for (i in 1:replications) {
    y  <- rbinom(length(x), 1, p)
    res  <-  try(glm(y~x, family=binomial), silent=TRUE)
    if (inherits(res,"try-error")) err <-  err+1 else good <- good+1
}
P_separation  <-  err/replications
P_separation

Podczas uruchamiania tego kodu prawdopodobieństwo separacji szacujemy na 0,759. Uruchom kod sam, jest szybki!

Następnie rozszerzamy ten kod, aby wypróbować różne procedury szacowania, mle i przybliżoną warunkową regresję logistyczną z elrm. Uruchomienie tej symulacji na moim komputerze zajmuje około 40 minut.

library(elrm)  # from CRAN
set.seed(31415926)  ### May I have a large container of coffee
replications  <-  1000
GOOD  <-  numeric(length=replications) ### will be set to one when MLE exists!
COEFS <- matrix(as.numeric(NA), replications, 2)
COEFS.elrm <- matrix(as.numeric(NA), replications, 2) # But we'll only use second col for x
p  <-  pltrue(a, b, x)
err  <-  0
good  <- 0

for (i in 1:replications) {
    y  <- rbinom(length(x), 1, p)
    res  <-  try(glm(y~x, family=binomial), silent=TRUE)
    if (inherits(res,"try-error")) err <-  err+1 else{ good <- good+1
                                                     GOOD[i] <- 1 }
    # Using stats::glm
    mod  <-  stats::glm(y~x, family=binomial)
    COEFS[i, ]  <-  coef(mod)
    # Using elrm:
    DATASET  <-  data.frame(x=x, y=y, n=1)
    mod.elrm  <-  elrm(y/n ~ x,  interest= ~ x -1, r=4, iter=10000, burnIn=1000,
                       dataset=DATASET)
    COEFS.elrm[i, 2 ]  <-  mod.erlm$coeffs       
}
### Now we can compare coefficient estimates of x,
###  when there are separation,  and when not:

non  <-  which(GOOD==1)
cof.mle.non  <-  COEFS[non, 2, drop=TRUE]
cof.mle.sep  <-  COEFS[-non, 2, drop=TRUE]
cof.elrm.non  <-  COEFS.elrm[non, 2, drop=TRUE]
cof.elrm.sep  <-  COEFS.elrm[-non, 2, drop=TRUE]

Teraz chcemy wykreślić wyniki, ale wcześniej zauważmy, że WSZYSTKIE szacunki warunkowe są równe! To jest naprawdę dziwne i powinno wymagać wyjaśnienia ... Wspólna wartość to 0,9523975. Ale przynajmniej uzyskaliśmy oszacowania skończone, z przedziałami ufności, które zawierają prawdziwą wartość (nie pokazano tutaj). Pokażę więc histogram szacunków mle tylko w przypadkach bez separacji:

hist(cof.mle.non, prob=TRUE)

[

Godne uwagi jest to, że wszystkie szacunki są mniejsze niż wartość rzeczywista 1,5. Może to mieć związek z faktem, że symulowaliśmy na podstawie zmodyfikowanego modelu, wymaga zbadania.

Tutaj są dobre odpowiedzi od @ sean501 i @kjetilbhalvorsen. Poprosiłeś o przykład. Rozważ poniższy rysunek. Można natknąć się sytuacji, w której proces generujący dane są tak przedstawione w panelu A . Jeśli tak, to jest całkiem możliwe, że dane które rzeczywiście wyglądają jak te gromadzą się w panelu B . Teraz, gdy używasz danych do budowy modelu statystycznego, chodzi o odzyskanie prawdziwego procesu generowania danych lub przynajmniej uzyskanie przybliżenia, które jest dość zbliżone. Pytanie zatem brzmi: czy dopasowanie regresji logistycznej do danych w B da model zbliżony do niebieskiej linii w A ? Jeśli spojrzysz na panel C, widać, że szara linia lepiej przybliża dane niż prawdziwa funkcja, więc przy szukaniu najlepszego dopasowania regresja logistyczna „woli” zwrócić szarą linię zamiast niebieskiej. Na tym jednak się nie kończy. Patrząc na panel D, czarna linia aproksymuje dane lepiej niż szara - w rzeczywistości jest to najlepsze dopasowanie, jakie może wystąpić. To jest linia, którą realizuje model regresji logistycznej. Odpowiada to punktowi przechwytywania ujemnej nieskończoności i nachyleniu nieskończoności. Jest to oczywiście bardzo dalekie od prawdy, którą masz nadzieję odzyskać. Całkowite rozdzielenie może również powodować problemy z obliczaniem wartości p dla zmiennych, które są standardowo dostarczane z wyjściem regresji logistycznej (wyjaśnienie jest nieco inne i bardziej skomplikowane). Co więcej, próba połączenia dopasowania tutaj z innymi próbami, na przykład z metaanalizą, sprawi, że inne ustalenia będą mniej dokładne.

Oznacza to, że istnieje hiperpłaszczyzna taka, że ​​po jednej stronie znajdują się wszystkie punkty dodatnie, a po drugiej wszystkie ujemne. Rozwiązanie największego prawdopodobieństwa to zatem płaska 1 z jednej strony i płaska 0 z drugiej strony, która jest „osiągana” za pomocą funkcji logistycznej poprzez posiadanie współczynników w nieskończoności.

To, co nazywacie „separacją” (a nie „oddzieleniem”), obejmuje dwie różne sytuacje, które w końcu powodują ten sam problem - którego nie nazwałbym jednak „niestabilnością”, jak to robicie.

Ilustracja: Przetrwanie na Titanicu

• $DV \in (0, 1)$$SV$

  $SV$$DV$$0$$1$

• $SV$$DV$

  Tak by było, gdyby wszyscy pasażerowie pierwszej klasy na Titanicu przeżyli wrak, a żaden z pasażerów drugiej klasy nie przeżył.

• $SV$$DV = 0$$DV = 1$

  Tak byłoby w przypadku niektórych pasażerów pierwszej klasy w Internecie$SV$$DV = 1$$DV = 0$

  Odwrotnie, jeśli tylko niektórzy pasażerowie drugiej klasy na Titanic zginęli we wraku, to klasa pasażerów$SV$$DV = 0$$DV = 1$

$DV$$SV$$SV$

Dlaczego regresja logistyczna jest w tych przypadkach „niestabilna”?

Jest to dobrze wyjaśnione w Rainey 2016 i Zorn 2005 .

• $DV$$1$$SV = 1$$DV$$0$$SV = 0$

  $SV = 1$

  $0$$1$$SV$$DV$

• $0$$1$$DV$$SV = 0$$SV = 1$

W obu przypadkach funkcja wiarygodności twojego modelu nie będzie w stanie znaleźć oszacowania maksymalnego prawdopodobieństwa: znajdzie przybliżenie tej wartości tylko poprzez zbliżenie się do niej asymptotycznie.

To, co nazywacie „niestabilnością”, to fakt, że w przypadkach całkowitej lub quasi-pełnej separacji nie ma skończonego prawdopodobieństwa, że ​​model logistyczny osiągnie. Nie użyłbym tego terminu: funkcja prawdopodobieństwa jest w rzeczywistości dość „stabilna” (monotoniczna) w przypisywaniu wartości współczynników do nieskończoności.

Uwaga: mój przykład jest fikcyjny. Przetrwanie na Titanicu nie sprowadzało się tylko do przynależności do klasy pasażerskiej. Patrz Hall (1986) .