¿Por qué la regresión logística se vuelve inestable cuando las clases están bien separadas?

¿Por qué la regresión logística se vuelve inestable cuando las clases están bien separadas? ¿Qué significan las clases bien separadas? Realmente agradecería si alguien puede explicar con un ejemplo.

No es correcto que la regresión logística en sí misma se vuelva inestable cuando hay separación. La separación significa que hay algunas variables que son muy buenos predictores, lo cual es bueno, o la separación puede ser un artefacto de muy pocas observaciones / demasiadas variables. Si ese es el caso, la solución podría ser obtener más datos. Pero la separación en sí misma, entonces, es solo un síntoma y no un problema en sí mismo.

Así que hay casos realmente diferentes para ser tratados. Primero, ¿cuál es el objetivo del análisis? Si el resultado final del análisis es alguna clasificación de casos, la separación no es un problema en absoluto, realmente significa que hay muy buenas variables que dan una clasificación muy buena. Pero si el objetivo es la estimación del riesgo, necesitamos las estimaciones de los parámetros, y con la separación no existen las estimaciones habituales de mle (máxima probabilidad). Entonces debemos cambiar el método de estimación, tal vez. Hay varias propuestas en la literatura, volveré sobre eso.

Luego hay (como se dijo anteriormente) dos posibles causas diferentes de separación. Puede haber separación en la población total, o la separación puede ser causada por pocos casos observados / demasiadas variables.

Lo que se rompe con la separación es el procedimiento de estimación de máxima verosimilitud. Las estimaciones de los parámetros mle (o al menos algunas de ellas) se vuelven infinitas. Dije en la primera versión de esta respuesta que eso se puede resolver fácilmente, tal vez con bootstrapping, pero eso no funciona, ya que habrá separación en cada remuestreo de bootstrap, al menos con el procedimiento habitual de bootstrapping. Pero la regresión logística sigue siendo un modelo válido, pero necesitamos algún otro procedimiento de estimación. Algunas propuestas han sido:

1. regularización, como cresta o lazo, tal vez combinada con bootstrap.
2. regresión logística condicional exacta
3. pruebas de permutación, consulte https://www.ncbi.nlm.nih.gov/pubmed/15515134
4. Procedimiento de estimación de sesgo reducido de Firths, ver El modelo de regresión logística no converge
5. seguramente otros ...

Si usa R, hay un paquete en CRAN, SafeBinaryRegressionque ayuda a diagnosticar problemas con la separación, ¡usando métodos matemáticos de optimización para verificar con seguridad si hay separación o cuasiseparación! A continuación, daré un ejemplo simulado usando este paquete y el elrmpaquete para la regresión logística condicional aproximada.

Primero, un ejemplo simple con el safeBinaryRegressionpaquete. Este paquete simplemente redefine la glmfunción, sobrecargándola con una prueba de separación, utilizando métodos de programación lineal. Si detecta separación, sale con una condición de error, declarando que el archivo no existe. De lo contrario, solo ejecuta la glmfunción ordinaria desde stats. El ejemplo es de sus páginas de ayuda:

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)

La salida de ejecutarlo:

> # 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

Ahora simulamos a partir de un modelo que puede ser aproximado por un modelo logístico, excepto que por encima de cierto límite la probabilidad de evento es exactamente 1.0. Piensa en un problema de bioensayo, pero por encima del límite, el veneno siempre mata:

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

Al ejecutar este código, estimamos la probabilidad de separación como 0.759. Ejecute el código usted mismo, ¡es rápido!

Luego extendemos este código para probar diferentes procedimientos de estimación, mle y regresión logística condicional aproximada de elrm. Ejecutar esta simulación toma alrededor de 40 minutos en mi computadora.

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]

Ahora queremos graficar los resultados, pero antes de eso, ¡tenga en cuenta que TODAS las estimaciones condicionales son iguales! Eso es realmente extraño y debería necesitar una explicación ... El valor común es 0.9523975. Pero al menos obtuvimos estimaciones finitas, con intervalos de confianza que contienen el valor verdadero (no se muestra aquí). Por lo tanto, solo mostraré un histograma de las estimaciones mle en los casos sin separación:

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

[

Lo notable es que todas las estimaciones son menores que el valor verdadero 1.5. Eso puede tener que ver con el hecho de que simulamos a partir de un modelo modificado, necesita investigación.

Aquí hay buenas respuestas de @ sean501 y @kjetilbhalvorsen. Pediste un ejemplo. Considere la figura a continuación. Usted puede venir a través de alguna situación en la que el proceso de generación de datos es como la que se muestra en el panel A . Si es así, es muy posible que los datos que se reúnen en realidad se parecen a los del panel B . Ahora, cuando usa datos para construir un modelo estadístico, la idea es recuperar el verdadero proceso de generación de datos o al menos llegar a una aproximación que sea razonablemente cercana. Por lo tanto, la pregunta es: ¿ajustar una regresión logística a los datos en B producirá un modelo que se aproxima a la línea azul en A ? Si nos fijamos en el panel C, puede ver que la línea gris se aproxima mejor a los datos que la función verdadera, por lo que al buscar el mejor ajuste, la regresión logística 'preferirá' devolver la línea gris en lugar de la azul. Sin embargo, no se detiene allí. Mirando el panel D, la línea negra aproxima los datos mejor que la gris; de hecho, es el mejor ajuste que podría ocurrir. Esa es la línea que sigue el modelo de regresión logística. Corresponde a una intersección de infinito negativo y una pendiente de infinito. Eso, por supuesto, está muy lejos de la verdad que espera recuperar. La separación completa también puede causar problemas con el cálculo de los valores p para las variables que vienen de manera estándar con la salida de regresión logística (la explicación allí es ligeramente diferente y más complicada). Además, tratar de combinar el ajuste aquí con otros intentos, por ejemplo con un metanálisis, hará que los otros hallazgos sean menos precisos.

Significa que hay un hiperplano tal que en un lado hay todos los puntos positivos y en el otro todos los negativos. La solución de máxima verosimilitud es entonces el plano 1 en un lado y el plano 0 en el otro lado, que se 'logra' con la función logística al tener los coeficientes en el infinito.

Lo que usted llama "separación" (no 'separación') cubre dos situaciones diferentes que terminan causando el mismo problema, que no llamaría, sin embargo, un problema de "inestabilidad" como usted.

Una ilustración: supervivencia en el Titanic

• Deje que sea ​​una variable dependiente binaria, y una variable separadora e independiente.$DV \in (0, 1)$$SV$

  Supongamos que es la clase de los pasajeros en el Titanic , y que indica si sobrevivieron a los restos, con indicando la muerte y indicando la supervivencia.$SV$$DV$$0$$1$

• La separación completa es la situación en la que predice todos los valores de .$SV$$DV$

  Ese sería el caso si todos los pasajeros de primera clase en el Titanic hubieran sobrevivido a los restos, y ninguno de los pasajeros de segunda clase hubiera sobrevivido.

• Separación-Quasi completa es la situación en la que predice bien donde todos los casos , o todos los casos en , pero no ambos.$SV$$DV = 0$$DV = 1$

  Ese sería el caso si algunos pasajeros de primera clase en el Titanic hubieran sobrevivido a los restos, y ninguno de los pasajeros de segunda clase hubiera sobrevivido. En ese caso, la clase de pasajeros predice todos los casos donde , pero no todos los casos donde .$SV$$DV = 1$$DV = 0$

  A la inversa, si solo algunos pasajeros de segunda clase en el Titanic hubieran muerto en los restos, entonces la clase de pasajeros predice todos los casos donde , pero no todos los casos donde , que incluye tanto pasajeros de primera clase como de segunda clase. .$SV$$DV = 0$$DV = 1$

Lo que usted llama "clases bien separadas" es la situación en la que una variable de resultado binaria (p. Ej., Supervivencia en el Titanic ) puede asignarse completa o casi completamente a un predictor (p. Ej., Pertenencia a la clase de pasajeros; no necesita ser binario como está en mi ejemplo).$DV$$SV$$SV$

¿Por qué la regresión logística es "inestable" en estos casos?

Esto se explica bien en Rainey 2016 y Zorn 2005 .

• Bajo separación completa , su modelo logístico buscará una curva logística que asigne, por ejemplo, todas las probabilidades de a cuando , y todas las probabilidades de a cuando .$DV$$1$$SV = 1$$DV$$0$$SV = 0$

  Esto corresponde a la situación antes mencionada donde solo y todos los pasajeros de primera clase del Titanic sobreviven, con indicando la membresía de pasajeros de primera clase.$SV = 1$

  Esto es problemático porque la curva logística se encuentra estrictamente entre y , lo que significa que, para modelar los datos observados, la maximización empujará algunos de sus términos hacia el infinito, para, si lo desea, hacer que "infinitamente" predictivo de .$0$$1$$SV$$DV$

• El mismo problema surge bajo una separación casi completa , ya que la curva logística todavía necesitará asignar solo valores de o a en uno de dos casos, o .1 D V S V = 0 S V = $0$$1$$DV$$SV = 0$$SV = 1$

En ambos casos, la función de probabilidad de su modelo será incapaz de encontrar una estimación de máxima probabilidad: solo encontrará una aproximación de ese valor al acercarse asintóticamente.

Lo que usted llama "inestabilidad" es el hecho de que, en casos de separación completa o casi completa, no existe una probabilidad finita de que el modelo logístico alcance. Sin embargo, no usaría ese término: la función de probabilidad es, de hecho, ser bastante "estable" (monotónica) en su asignación de valores de coeficientes hacia el infinito.

Nota: mi ejemplo es ficticio. La supervivencia en el Titanic no se redujo solo a la membresía de la clase de pasajeros. Ver Hall (1986) .