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

163

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

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

Todavía obtenemos el modelo, pero los coeficientes estimados están inflados.

¿Cómo manejas esto en la práctica?

usuario333
fuente
44
pregunta
usuario603
1
pregunta relacionada y demostración sobre la regularización aquí
Haitao Du

Respuestas:

100

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 penalizadas (aunque resultaron tener otras propiedades interesantes).

Instale y cargue el paquete glmnet en R y estará listo para comenzar. Uno de los aspectos menos fáciles de usar de glmnet es que solo puede alimentarlo con matrices, no con fórmulas como estamos acostumbrados. Sin embargo, puede mirar model.matrix y similares para construir esta matriz a partir de un marco de datos y una fórmula ...

Ahora, cuando espera que esta separación perfecta no sea solo un subproducto de su muestra, sino que podría ser cierto en la población, específicamente no desea manejar esto: use esta variable de separación simplemente como el único predictor de su resultado, no empleando un modelo de cualquier tipo.

Nick Sabbe
fuente
20
También puede usar una interfaz de fórmula para glmnet a través del paquete caret.
Zach
"Ahora, cuando esperas ..." Pregunta sobre esto. Tengo un estudio de caso / control que analiza la relación con el microbioma. También tenemos un tratamiento que casi solo se encuentra entre los casos. Sin embargo, creemos que el tratamiento también podría afectar el microbioma. ¿Es este un ejemplo de su advertencia? Hipotéticamente, podríamos encontrar muchos más casos sin el tratamiento si lo intentáramos, pero tenemos lo que tenemos.
abalter
142

Tienes varias opciones:

  1. Eliminar algunos de los prejuicios.

    (a) Al penalizar la probabilidad según la sugerencia de @ Nick. El paquete logistf en R o la FIRTHopción en SAS PROC LOGISTICimplementan el método propuesto en Firth (1993), "Reducción de sesgo de las estimaciones de máxima verosimilitud", Biometrika , 80 , 1 .; que elimina el sesgo de primer orden de las estimaciones de máxima verosimilitud. ( Aquí @Gavin recomienda el brglmpaquete, con el que no estoy familiarizado, pero creo que implementa un enfoque similar para las funciones de enlace no canónicas, por ejemplo, probit).

    (b) Mediante el uso de estimaciones no sesgadas medianas en regresión logística condicional exacta. Paquete elrm o logistiX en R, o la EXACTdeclaración en SAS PROC LOGISTIC.

  2. Excluya los casos en que se produce la categoría o el valor del predictor que causa la separación. Estos pueden estar fuera de su alcance; o digno de más investigación centrada. (El paquete R safeBinaryRegression es útil para encontrarlos).

  3. Vuelve a lanzar el modelo. Por lo general, esto es algo que habría hecho de antemano si lo hubiera pensado, porque es demasiado complejo para el tamaño de su muestra.

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

    (b) Al contraer categorías de predictores / agrupación de los valores de predictores. Solo si esto tiene sentido.

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

  4. 5212

  5. Hacer nada. (Pero calcule los intervalos de confianza en función de las probabilidades de perfil, ya que las estimaciones de Wald del error estándar serán muy incorrectas). Una opción que a menudo se pasa por alto. Si el propósito del modelo es solo describir lo que ha aprendido sobre las relaciones entre los predictores y la respuesta, no hay vergüenza en citar un intervalo de confianza para una razón de probabilidades de, digamos, 2.3 hacia arriba. (De hecho, podría parecer sospechoso citar intervalos de confianza basados ​​en estimaciones imparciales que excluyen los odds ratios mejor respaldados por los datos). Los problemas surgen cuando intenta predecir el uso de estimaciones puntuales, y el predictor en el que se produce la separación inunda a los demás.

  6. Utilice un modelo de regresión logística oculto, como se describe en Rousseeuw y Christmann (2003), "Robustez contra la separación y valores atípicos en la regresión logística", Computational Statistics & Data Analysis , 43 , 3, e implementado en el paquete R hlr . (@ user603 sugiere esto ) . No he leído el documento, pero dicen en resumen "se propone un modelo un poco más general bajo el cual la respuesta observada está fuertemente relacionada pero no es igual a la respuesta verdadera no observable", que sugiere Para mí, no sería una buena idea usar el método a menos que parezca plausible.

  7. "Cambie algunas observaciones seleccionadas al azar de 1 a 0 o de 0 a 1 entre las variables que muestran una separación completa": comentario de @ RobertF . Esta sugerencia parece surgir de considerar la separación como un problema per se en lugar de como un síntoma de escasez de información en los datos que podría llevarlo a preferir otros métodos a la estimación de máxima verosimilitud, o limitar las inferencias a las que puede hacer con precisión razonable: enfoques que tienen sus propios méritos y no son solo "soluciones" para la separación. (Además de ser descaradamente ad hoc , es desagradable para la mayoría de los analistas que hacen la misma pregunta de los mismos datos, hacen los mismos supuestos, deberían dar diferentes respuestas debido al resultado de un lanzamiento de moneda o lo que sea).

Scortchi
fuente
1
@Scortchi Hay otra opción (herética). ¿Qué pasa con el cambio de algunas observaciones seleccionadas al azar de 1 a 0 o de 0 a 1 entre las variables que muestran una separación completa?
RobertF
@RobertF: ¡Gracias! No había pensado en este: si tienes alguna referencia sobre su rendimiento, te lo agradecería. ¿Te has encontrado con personas que lo usan en la práctica?
Scortchi
@Scortchi: No, hay referencias a investigadores que agregan datos artificiales para eliminar la separación completa, pero no he encontrado ningún artículo sobre la modificación selectiva de los datos. No tengo idea de cuán efectivo sería este método.
RobertF
1
@tatami: No todos los programas (¿muchos?) advierten sobre la separación per se, que puede ser difícil de detectar cuando se trata de una combinación lineal de varias variables, sino sobre la falla de convergencia y / o valores ajustados cercanos a la nada o uno. siempre revisa estos.
Scortchi
2
@Scortchi: muy buen resumen en su respuesta. Personalmente, estoy a favor del enfoque bayesiano, pero vale la pena mencionar el hermoso análisis del fenómeno general desde un punto de vista frecuente en projecteuclid.org/euclid.ejs/1239716414 . El autor ofrece algunos intervalos de confianza unilaterales que pueden usarse incluso en presencia de una separación completa en la regresión logística.
Cian
55

Esta es una expansión de las respuestas de Scortchi y Manoel, pero como parece que usas RI, pensé que te daría un código. :)

Creo que la solución más sencilla y directa a su problema es utilizar un análisis bayesiano con suposiciones previas no informativas como lo propuso Gelman et al (2008). Como menciona Scortchi, Gelman recomienda poner un Cauchy antes con una mediana de 0.0 y una escala de 2.5 en cada coeficiente (normalizado para tener una media de 0.0 y una DE de 0.5). Esto regularizará los coeficientes y los arrastrará ligeramente hacia cero. En este caso es exactamente lo que quieres. Debido a que tiene colas muy anchas, el Cauchy aún permite coeficientes grandes (en oposición al Normal de cola corta), de Gelman:

ingrese la descripción de la imagen aquí

¿Cómo ejecutar este análisis? ¡Use la bayesglmfunción en el paquete arm 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 simple, no?

Referencias

Gelman et al (2008), "Una distribución previa por defecto débilmente informativa para modelos logísticos y otros modelos de regresión", Ann. Appl. Stat., 2, 4 http://projecteuclid.org/euclid.aoas/1231424214

Rasmus Bååth
fuente
66
No. Demasiado simple. ¿Puedes explicar lo que acabas de hacer? ¿Cuál es el previo que bayesglmusa? Si la estimación de ML es equivalente a Bayesian con un previo plano, ¿cómo ayudan aquí los previos no informativos?
StasK
55
¡Agregué más información! Lo anterior es vago pero no plano. Tiene cierta influencia, ya que regulariza las estimaciones y las empuja ligeramente hacia 0.0, que es lo que creo que quiere en este caso.
Rasmus Bååth
> m = bayesglm (match ~., family = binomial (link = 'logit'), data = df) Mensaje de advertencia: probabilidades ajustadas numéricamente 0 o 1 ocurrieron ¡ No es bueno!
Chris
Como titular, intente una regularización un poco más fuerte aumentando prior.dfel valor predeterminado 1.0y / o disminuyendo prior.scaleel valor predeterminado 2.5, tal vez comience a intentar:m=bayesglm(match ~. , family = binomial(link = 'logit'), data = df, prior.df=5)
Rasmus Bååth
1
¿Qué estamos haciendo exactamente cuando aumentamos prior.df en el modelo? ¿Hay un límite de cuán alto queremos llegar? ¿Tengo entendido que limita el modelo para permitir la convergencia con estimaciones precisas de error?
hamilthj
7

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

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

  • La separación casi completa se produce cuando (a) existe algún vector de coeficiente b tal que bxi ≥ 0 siempre que yi = 1 , y bxi ≤ 0 * siempre ** yi = 0 y esta igualdad se cumple al menos para un caso en cada categoría de la variable dependiente. En otras palabras, en el caso más simple, 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 regresión no existe.

Allison analiza muchas de las soluciones ya mencionadas, incluida la eliminación de variables problemáticas, el colapso de categorías, no hacer nada, aprovechar la regresión logística exacta , la estimación bayesiana y la estimación de máxima probabilidad penalizada.

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

Mike Hunter
fuente
3

warning

Con datos generados a lo largo de las líneas de

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

Se hace la advertencia:

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

lo que obviamente refleja la dependencia que está integrada en estos datos.

En R, la prueba de Wald se encuentra con summary.glmo con waldtesten el lmtestpaquete. La prueba de razón de probabilidad se realiza con anovao con lrtestel lmtestpaquete. En ambos casos, la matriz de información tiene un valor infinito y no hay inferencia disponible. Más bien, R no producir una salida, pero no se puede confiar en ella. La inferencia que R produce típicamente en estos casos tiene valores p muy cercanos a uno. Esto se debe a que la pérdida de precisión en el OR es un orden de magnitud menor que la pérdida de precisión en la matriz de varianza-covarianza.

Algunas soluciones descritas aquí:

Use un estimador de un paso,

Existe mucha teoría que respalda el bajo sesgo, la eficiencia y la generalización de los estimadores de un paso. Es fácil especificar un estimador de un paso en R y los resultados suelen ser muy favorables para la predicción y la inferencia. ¡Y este modelo nunca divergerá, porque el iterador (Newton-Raphson) simplemente no tiene la oportunidad de hacerlo!

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

Da:

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

Para que pueda ver las predicciones reflejan la dirección de la tendencia. Y la inferencia es muy sugerente de las tendencias que creemos que son ciertas.

ingrese la descripción de la imagen aquí

realizar una prueba de puntaje,

La estadística de Puntuación (o Rao) difiere de la razón de probabilidad y las estadísticas de wald. No requiere una evaluación de la varianza bajo la hipótesis alternativa. Encajamos el modelo debajo del nulo:

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

En ambos casos tienes inferencia para un OR de infinito.

, y utilice estimaciones medias imparciales para un intervalo de confianza.

Puede generar un IC del 95% imparcial, no singular para la razón de probabilidades infinita utilizando la estimación imparcial mediana. El paquete epitoolsen R puede hacer esto. Y doy un ejemplo de implementación de este estimador aquí: intervalo de confianza para el muestreo de Bernoulli

AdamO
fuente
2
Esto es genial, pero tengo algunas objeciones, por supuesto: (1) La prueba de razón de probabilidad no utiliza la matriz de información; es solo la prueba de Wald que lo hace, y que falla catastróficamente en presencia de separación. (2) No estoy familiarizado con los estimadores de un solo paso, pero la estimación de la pendiente aquí parece absurdamente baja. (3) Un intervalo de confianza no es medianamente imparcial. Lo que vincula en esa sección es el intervalo de confianza a mitad de p. (4) Puede obtener intervalos de confianza invirtiendo el LR o las pruebas de puntaje. ...
Scortchi
... (5) Puede realizar la prueba de puntuación en R dando el argumento test="Rao"a la anovafunción. (Bueno, las dos últimas son notas, no objeciones.)
Scortchi
¡@scortchi es bueno saber que anova tiene pruebas de puntaje predeterminadas! Quizás una implementación manual sea útil. Los IC no son medianos insesgados, pero los CI para el estimador insesgado mediano proporcionan inferencia consistente para los parámetros de límite. El p medio es un estimador de este tipo. El p puede transformarse en una razón de probabilidad b / c, es invariable para transformaciones uno a uno. ¿La prueba LR es consistente para los parámetros de límite?
AdamO
Solo la hipótesis nula no debe contener parámetros en un límite para que se aplique el teorema de Wilks, aunque las pruebas de puntuación y LR son aproximadas en muestras finitas.
Scortchi
2

Tenga cuidado con este mensaje de advertencia de R. Eche un vistazo a esta publicación de blog de Andrew Gelman, y verá que no siempre es un problema de separación perfecta, sino que a veces es un error glm. Parece que si los valores iniciales están demasiado lejos de la estimación de máxima verosimilitud, explota. Entonces, verifique primero con otro software, como Stata.

Si realmente tiene este problema, puede intentar utilizar el modelado bayesiano, con antecedentes informativos.

Pero en la práctica, me deshago de los predictores que causan el problema, porque no sé cómo elegir un informe informativo. Pero supongo que hay un documento de Gelman sobre el uso de información previa cuando tienes este problema de separación perfecta. Solo googlealo. Tal vez deberías intentarlo.

Manoel Galdino
fuente
8
El problema con la eliminación de predictores es que está eliminando el predictor que mejor explica la respuesta, que generalmente es lo que pretende hacer. Yo diría que esto solo tiene sentido si ha sobreajustado su modelo, por ejemplo, ajustando demasiadas interacciones complicadas.
Simon Byrne
44
No es un error, sino un problema con las estimaciones iniciales que están demasiado lejos del MLE, lo que no surgirá si no intentas elegirlas tú mismo.
Scortchi
Entiendo esto, pero creo que esto es un error en el algoritmo.
Manoel Galdino
55
Bueno, no quiero discutir sobre la definición de 'error'. Pero el comportamiento no es insondable ni inamovible en la base R: no es necesario "consultar con otro software". Si desea lidiar automáticamente con muchos problemas de no convergencia, el glm2paquete implementa una verificación de que la probabilidad realmente aumenta en cada paso de puntuación y reduce a la mitad el tamaño del paso si no es así.
Scortchi
3
Existe (en CRAN) el paquete R safeBinaryRegression que está diseñado para diagnosticar y solucionar estos problemas, utilizando métodos de optimización para verificar con seguridad si hay separación o cuasiseparación. ¡Intentalo!
kjetil b halvorsen
2

No estoy seguro de estar de acuerdo con las declaraciones en su pregunta.

Creo que ese mensaje de advertencia significa que, para algunos de los niveles de X observados en sus datos, la probabilidad ajustada es numéricamente 0 o 1. En otras palabras, en la resolución, se muestra como 0 o 1.

Puedes correr predict(yourmodel,yourdata,type='response')y encontrarás 0's y / y 1's allí como probabilidades pronosticadas.

Como resultado, creo que está bien usar los resultados.

StayLearning
fuente
-1

Entiendo que esta es una publicación antigua, sin embargo, continuaré respondiendo esto, ya que he luchado días con ella y puede ayudar a otros.

La separación completa ocurre cuando las variables seleccionadas para ajustarse al modelo pueden diferenciar con mucha precisión entre 0 y 1 o sí y no. Todo nuestro enfoque de la ciencia de datos se basa en la estimación de probabilidad, pero falla en este caso.

Pasos de rectificación: -

  1. Use bayesglm () en lugar de glm (), en caso de que la varianza entre las variables sea baja

  2. A veces, usar (maxit = "algún valor numérico") junto con bayesglm () puede ayudar

3. La tercera y más importante verificación para las variables seleccionadas para el ajuste del modelo, debe haber una variable para la cual la colinealidad múltiple con la variable Y (outout) es muy alta, descarte esa variable de su modelo.

Como en mi caso, tenía datos de abandono de telecomunicaciones para predecir el abandono de los datos de validación. Tenía una variable en mis datos de entrenamiento que podía diferenciar mucho entre el sí y el no. Después de soltarlo, pude obtener el modelo correcto. Además, puede usar paso a paso (ajuste) para que su modelo sea más preciso.

yash
fuente
2
No veo que esta respuesta agregue mucho a la discusión. El enfoque bayesiano está completamente cubierto en respuestas anteriores, la eliminación de predictores "problemáticos" también ya se menciona (y se desaconseja). La selección de variables paso a paso rara vez es una gran idea hasta donde yo sé.
Einar