¿Cómo tratar la separación perfecta en la regresión logística?

Question

Si tiene una variable que separa perfectamente los ceros y los unos en la variable de destino, R dará el siguiente mensaje de advertencia de "separación perfecta o casi perfecta":

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

Seguimos obteniendo el modelo, pero las estimaciones de los coeficientes están infladas.

¿Cómo se aborda esto en la práctica?

Answer 1

Tienes varias opciones:

1. Quitar parte de la parcialidad.

   (a) Penalizando la probabilidad según la sugerencia de @Nick. Paquete logísticaf en R o el FIRTH de SAS PROC LOGISTIC aplicar el método propuesto en Firth (1993), "Bias reduction of maximum likelihood estimates", Biometrika , 80 ,1.; que elimina el sesgo de primer orden de las estimaciones de máxima verosimilitud. ( Aquí @Gavin recomienda el brglm con el que no estoy familiarizado, pero deduzco que implementa un enfoque similar para las funciones de enlace no canónicas, por ejemplo, probit).

   (b) Utilizando estimaciones insesgadas por la mediana en la regresión logística condicional exacta. Paquete elrm o logistiX en R, o el EXACT de SAS PROC LOGISTIC .

2. Excluye casos donde se produce la categoría o valor predictor que provoca la separación. Es posible que esto quede fuera de su alcance, o que merezca una investigación más detallada. (El paquete R safeBinaryRegression es útil para encontrarlos).

3. Vuelve a fundir el modelo. Normalmente esto es algo que habrías hecho de antemano si hubieras pensado en ello, porque es demasiado complejo para el tamaño de tu muestra.

   (a) Retire el predictor del modelo. Dicey, por las razones dado por @Simon: "Estás eliminando el predictor que mejor explica la respuesta".

   (b) Colapsando las categorías de los predictores / agrupando los valores de los predictores. Sólo si esto tiene sentido.

   (c) Reexpresar el predictor como dos (o más) factores cruzados sin interacción. Sólo si esto tiene sentido.

4. Utilice un análisis bayesiano como el de @Manoel sugerencia . Aunque parece poco probable que quieras sólo debido a la separación, vale la pena considerarlo por sus otros méritos.El papel que recomienda es Gelman et al (2008), "A weakly informative default prior distribution for logistic & other regression models", Ann. Appl. Stat. , 2 , 4 el valor predeterminado en cuestión es una prioridad Cauchy independiente para cada coeficiente, con una media de cero y una escala de $\frac{5}{2}$ que se utilizará después de estandarizar todos los predictores continuos para que tengan una media de cero y una desviación estándar de $\frac{1}{2}$ . Si se pueden dilucidar los antecedentes fuertemente informativos, mucho mejor.

5. No hacer nada. (Pero calcular los intervalos de confianza basados en las probabilidades de los perfiles, ya que las estimaciones de Wald del error estándar serán muy erróneas). Una opción que a menudo se pasa por alto. Si el propósito del modelo es sólo describir lo que se ha aprendido sobre las relaciones entre los predictores y la respuesta, no hay que avergonzarse de citar un intervalo de confianza para una razón de momios de, digamos, 2,3 en adelante. (De hecho, podría parecer sospechoso citar intervalos de confianza basados en estimaciones insesgadas que excluyan los cocientes de probabilidades mejor respaldados por los datos). Los problemas surgen cuando se intenta predecir mediante estimaciones puntuales, y el predictor sobre el que se produce la separación empapa a los demás.

6. Utilice un modelo de regresión logística oculta, tal y como se describe en Rousseeuw y Christmann (2003), "Robustness against separation and outliers in logistic regression", Estadística computacional y análisis de datos , 43 , 3, e implementado en el paquete R hlr . (@usuario603 sugiere esto. ) No he leído el artículo, pero en el resumen dicen que "se propone un modelo ligeramente más general bajo el cual la respuesta observada está fuertemente relacionada pero no es igual a la respuesta verdadera no observable", lo que me sugiere que podría no ser una buena idea utilizar el método a menos que suene plausible.

7. "Cambiar algunas observaciones seleccionadas al azar de 1 a 0 o de 0 a 1 entre las variables que muestran una separación completa": @RobertF's comentario . Esta sugerencia parece surgir de considerar la separación como un problema por sí mismo y no como un síntoma de una escasez de información en los datos que podría llevarle a preferir otros métodos a la estimación de máxima verosimilitud, o a limitar las inferencias a las que pueda hacer con una precisión razonable, enfoques que tienen sus propios méritos y no son sólo "arreglos" para la separación. (Aparte de su carácter descaradamente ad hoc En este sentido, la mayoría de los analistas que se plantean la misma pregunta sobre los mismos datos y hacen las mismas suposiciones deben dar respuestas diferentes debido al resultado de un lanzamiento de una moneda o lo que sea).

Answer 2

Una solución a esto es utilizar una forma de regresión penalizada. De hecho, esta es la razón original por la que se desarrollaron algunas de las formas de regresión penalizada (aunque resultaron tener otras propiedades interesantes.

Instale y cargue el paquete glmnet en R y estará casi listo para empezar. Uno de los aspectos menos fáciles de usar de glmnet es que sólo puedes alimentarlo con matrices, no con fórmulas como estamos acostumbrados. Sin embargo, puedes consultar model.matrix y similares para construir esta matriz a partir de un data.frame y una fórmula...

Ahora bien, cuando usted espera que esta separación perfecta no sea sólo un subproducto de su muestra, sino que pueda ser cierta en la población, usted específicamente no quiere manejar esto: utilizar esta variable de separación simplemente como el único predictor de su resultado, sin emplear un modelo de ningún tipo.

Answer 3

Esto es una ampliación de las respuestas de Scortchi y Manoel, pero como parece que usas R he pensado en aportar algo de código. :)

Creo que la solución más fácil y directa a su problema es utilizar un análisis bayesiano con supuestos a priori no informativos, como proponen Gelman et al (2008). Como menciona Scortchi, Gelman recomienda poner una previa de Cauchy con una mediana de 0,0 y una escala de 2,5 en cada coeficiente (normalizada para tener una media de 0,0 y una DE de 0,5). Esto regularizará los coeficientes y los acercará ligeramente a cero. En este caso es exactamente lo que quieres. Debido a que tiene colas muy anchas, el Cauchy todavía permite coeficientes grandes (a diferencia de la Normal de cola corta), de Gelman:

¿Cómo realizar este análisis? Utilice la herramienta bayesglm función en paquete de brazos ¡que implementa este análisis!

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

No funciona tan bien... Ahora la versión bayesiana:

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)

Súper sencillo, ¿no?

Referencias

Gelman et al (2008), "A weakly informative default prior distribution for logistic & other regression models", Ann. Appl. Stat., 2, 4 http://projecteuclid.org/euclid.aoas/1231424214

Answer 4

Una de las explicaciones más completas de los problemas de "separación cuasi-completa" en la máxima probabilidad es el artículo de Paul Allison. Escribe sobre el software SAS, pero las cuestiones que aborda son generalizables a cualquier software:

• La separación completa se produce cuando una función lineal de x puede generar predicciones perfectas de y

• La separación cuasi-completa se produce cuando (a) existe algún vector de coeficientes b tal que bxi 0 siempre que yi = 1 , y bxi 0* siempre que **yi = 0 y esta igualdad se mantiene para at al menos un caso en cada categoría de la variable dependiente. En otras palabras, en el caso más sencillo, para cualquier variable independiente dicotómica en una regresión logística, si hay un cero en la tabla 2 × 2 formada por esa variable y la variable dependiente, la estimación ML para el coeficiente de coeficiente de regresión no existe.

Allison discute muchas de las soluciones ya mencionadas, incluyendo la eliminación de las variables del problema, el colapso de las categorías, no hacer nada, aprovechar exacto regresión logística, estimación bayesiana y estimación de máxima verosimilitud penalizada.

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

Answer 5

La pregunta original está mal planteada y muchas de las respuestas son problemáticas. El hecho de que una estimación de máxima verosimilitud sea $\infty$ cuando hay una separación perfecta sólo es un problema porque seguimos utilizando los estadísticos de Wald (es decir, utilizamos la matriz de información y los errores estándar) para la inferencia. Un $\infty$ $\beta$ da lugar a una probabilidad predicha de 1,0. Esto no tiene nada de malo, aunque es probable que los modelos bayesianos o la contracción en un modelo frecuentista den como resultado un modelo mejor calibrado. Sólo hay que utilizar la razón de verosimilitud $\chi^2$ y los intervalos de confianza de la probabilidad del perfil y obtendrá una inferencia válida sin cambiar el modelo. Véase, por ejemplo, este paquete de R: https://cran.r-project.org/web/packages/ProfileLikelihood/ProfileLikelihood.pdf .

Creo que deberíamos utilizar rutinariamente modelos bayesianos, pero reconozcamos que $\infty$ es un MLE válido.