Problema 1: advertencia de convergencia, la optimización falló
Usando el pitón modelos de estadísticas paquete, puede ejecutar un LR básico y recibir la advertencia “Advertencia de convergencia: la optimización de máxima probabilidad no pudo converger”.
import statsmodels.formula.api as smfmodel = smf.logit("target_satisfaction ~ C(age,Treatment('0-13')) + \
C(educ,Treatment('Some High School')) + \
C(emp,Treatment('Not Employed')) + \
C(gender,Treatment('Male')) + \
C(hispanic,Treatment('N')) + \
C(dept,Treatment('Medical Oncology')) + \
C(sign_status,Treatment('Activated')) + \
C(mode,Treatment('Internet')) + \
C(inf,Treatment('Missing'))", data=df_a)
results = model.fit(method='bfgs')
print(results.summary())
Los coeficientes y errores estándar pueden incluso poblarse, como es normal, incluidos los de Regresión logística en el sklearn.linear_model paquete.
Dados los resultados, es posible que tenga la tentación de detenerse aquí. Pero no lo hagas. Quiere asegurarse de que su modelo converja para producir la mejor función de costo (la más baja) y el mejor ajuste del modelo. [1].
Puedes resolver esto en modelos de estadísticas o aprender cambiando el solucionador/método o aumentando el maxiter parámetro. En aprenderel parámetro de sintonización C También se puede reducir para aplicar una mayor regularización L2 y se puede probar de forma iterativa en una escala logarítmica (100, 10, 1, 0,1, 0,01, etc.) [2].
En mi caso particular, aumenté a maxiter= 300000 y el modelo convergió. La razón por la que funcionó es porque le proporcioné al modelo (por ejemplo, solucionadores) más intentos de converger y encontrar los mejores parámetros. [3]. Y esos parámetros, como los coeficientes y los valores p, de hecho se vuelven más precisos.
Problema 2: se agregó una función, pero las salidas LR no se actualizaron
Éste es fácil de pasar por alto, pero fácil de diagnosticar. Si está trabajando con muchas funciones o no las detectó en la limpieza de datos, puede incluir accidentalmente una característica categórica en su modelo LR que sea casi constante o que tenga solo un nivel… cosas malas. La inclusión de dicha característica generará coeficientes y errores estándar sin ninguna advertencia.
Aquí está nuestro modelo y resultado sin la característica mala.
model = smf.logit("target_satisfaction ~ C(age,Treatment('0-13')) + \
C(educ,Treatment('Some High School')) + \
C(emp,Treatment('Not Employed')) + \
C(gender,Treatment('Male')) + \
C(hispanic,Treatment('N')) + \
C(dept,Treatment('Medical Oncology')) + \
C(sign_status,Treatment('Activated')) + \
C(mode,Treatment('Internet')) + \
C(inf,Treatment('Missing'))", data=df_a)results = model.fit(method='bfgs',maxiter=30000)
print(results.summary())
Ahora, agreguemos una característica de un nivel. tipocon su categoría de referencia establecida en “Desaparecido”.
model = smf.logit("target_satisfaction ~ C(age,Treatment('0-13')) + \
C(educ,Treatment('Some High School')) + \
C(emp,Treatment('Not Employed')) + \
C(gender,Treatment('Male')) + \
C(hispanic,Treatment('N')) + \
C(dept,Treatment('Medical Oncology')) + \
C(sign_status,Treatment('Activated')) + \
C(mode,Treatment('Internet')) + \
C(inf,Treatment('Missing')) + \
C(type,Treatment('Missing'))", data=df_a) results = model.fit(method='bfgs', maxiter=300000)
print(results.summary())
Y aquí vemos el impacto de esa característica superflua… que no fue impacto alguno. El valor de R cuadrado permanece igual, resaltado arriba y abajo.
La regresión logística no puede realizar estimaciones significativas sobre una característica con un nivel o valores constantes y puede descartarla del modelo. Esto no sucede sin ningún factor determinante para el modelo en sí, pero aún así, la mejor práctica es incluir solo las características en el modelo que están teniendo un impacto positivo.
Un simple valor_cuenta() de nuestra característica tipo revela que tiene un nivel, lo que indica que desea eliminar la característica de su modelo.
Problema 3: “Error al invertir el hessiano” (y errores estándar de NaN)
Este es un problema importante y sucede con frecuencia. Creemos este nuevo error incluyendo cuatro características más en nuestro modelo LR: VA, VB, VC y VI. Todos son categóricos, cada uno con 3 niveles (0, 1 y “Desaparecido” como referencia).
model = smf.logit("target_satisfaction ~ C(age,Treatment('0-13')) + \
C(educ,Treatment('Some High School')) + \
C(emp,Treatment('Not Employed')) + \
C(gender,Treatment('Male')) + \
C(hispanic,Treatment('N')) + \
C(dept,Treatment('Medical Oncology')) + \
C(sign_status,Treatment('Activated')) + \
C(mode,Treatment('Internet')) + \
C(inf,Treatment('Missing')) + \
C(type,Treatment('Missing')) + \
C(VA,Treatment('Missing')) + \
C(VB,Treatment('Missing')) + \
C(VC,Treatment('Missing')) + \
C(VI,Treatment('Missing'))", data=df_a) results = model.fit(method='bfgs', maxiter=300000)
print(results.summary())
Y aquí está nuestro nuevo error:
Al explorar un poco los resultados, vemos los coeficientes representados, pero ninguno de los errores estándar o valores p. Sklearn proporciona coeficientes también.
Pero, ¿por qué los valores de NaN? ¿Y por qué es necesario invertir el hessiano?
Muchos algoritmos de optimización utilizan el inverso de Hesse (o una estimación del mismo) para maximizar la función de verosimilitud. Entonces, cuando la inversión falla, significa que el hessiano (técnicamente, una segunda derivada de la probabilidad logarítmica) no es un definido positivo. [4]. Visualmente, algunas características tienen curvaturas enormes mientras que otras tienen curvaturas más pequeñas, y el resultado es que el modelo no puede producir errores estándar para los parámetros.
Más concretamente, cuando un hessiano no es invertible, ningún cambio computacional puede hacerlo inverso porque simplemente no existe. [5].
¿Cuáles son las causas de esto? Hay tres más comunes:
- Hay más variables características que observaciones.
- Las características categóricas tienen niveles de frecuencia muy bajos.
- Existe multicolinealidad entre características.
Dado que nuestro conjunto de datos tiene muchas filas (~25 000) en relación con las características (14), podemos ignorar con seguridad la primera posibilidad (aunque existen soluciones para este problema exacto). [6]). La tercera posibilidad puede estar en juego, y podemos comprobarlo con el factor de inflación de varianza (VIF). La segunda posibilidad es un poco más fácil de diagnosticar, así que comencemos por ahí.
Tenga en cuenta que las características VA, VB, VC y VI se componen principalmente de 1 y 0.
En particular, la característica VI tiene una frecuencia (relativa) muy baja para la categoría “Desaparecido”, en realidad del orden del 0,03% (9/24.874).
Digamos que confirmamos con nuestro contexto empresarial que podemos colapsar 1 y “Falta” juntos. O, al menos, cualquier posible consecuencia de refactorizar los datos de esta manera sería mucho menor que aceptar un modelo con errores conocidos (y sin errores estándar).
Entonces creamos VI_alt, que tiene 2 niveles y 0 puede servir como referencia.
Intercambiando VI_alt por VI en el modelo,
model = smf.logit("target_satisfaction ~ C(age,Treatment('0-13')) + \
C(educ,Treatment('Some High School')) + \
C(emp,Treatment('Not Employed')) + \
C(gender,Treatment('Male')) + \
C(hispanic,Treatment('N')) + \
C(dept,Treatment('Medical Oncology')) + \
C(sign_status,Treatment('Activated')) + \
C(mode,Treatment('Internet')) + \
C(inf,Treatment('Missing')) + \
C(type,Treatment('Missing')) + \
C(VA,Treatment('Missing')) + \
C(VB,Treatment('Missing')) + \
C(VC,Treatment('Missing')) + \
C(VI_alt,Treatment(0))", data=df_a) results = model.fit(method='bfgs', maxiter=300000)
print(results.summary())
Hay una ligera mejora en el modelo general porque converge sin errores, ahora aparecen los coeficientes y los errores estándar. Una vez más, sigue siendo un modelo mal ajustado, pero ahora funciona modelo, que es nuestro objetivo aquí.
Nuestra tercera y última posibilidad de inversión fallida de Hesse es la multicolinealidad. Ahora, la multicolinealidad es un tema candente en el aprendizaje automático, creo que porque (1) afecta a muchos modelos populares como LR, regresión lineal, KNN y Naive Bayes (2) la heurística VIF para eliminar características resulta ser más un arte que ciencia y (3) en última instancia, los profesionales no están de acuerdo sobre si eliminar o no tales características en primer lugar si se hace a expensas de inyectar un sesgo de selección o perder información clave del modelo. [5–9].
No bajaré por esa madriguera del conejo. Es profundo.
Pero mostraré cómo calcular los VIF y tal vez podamos sacar nuestras propias conclusiones.
En primer lugar, el VIF realmente pregunta: “¿Qué tan bien se explica una de las características en conjunto con todas mis otras características?” La estimación del VIF de una característica se “inflará” cuando exista una dependencia lineal con otras características.
Dadas todas las características de nuestro conjunto de datos, calculemos los VIF a continuación.
from statsmodels.stats.outliers_influence import variance_inflation_factorvariables = results.model.exog
vif = [variance_inflation_factor(variables, i) for i in range(variables.shape[1])]
vifs = pd.DataFrame({'variables':results.model.exog_names,'vif':[ '%.2f' % elem for elem in vif ]})
vifs.sort_index(ascending=False).head(14)
Tenga en cuenta que en esta lista parcial, las características VA, VB, VD muestran VIF hacia el infinito, lo que supera con creces el umbral de “regla general” de 5. Pero debemos tener cuidado con heurísticas como esta, ya que los umbrales VIF tienen dos advertencias:
- El límite de 5 es relativo a otras características; por ejemplo, si la gran mayoría de los VIF de características caen por debajo de 7 y un subconjunto más pequeño de VIF de características está por encima de 7, entonces 7 podría ser un umbral más razonable.
- Las características categóricas con un número menor de casos en la categoría de referencia en relación con otros niveles producirán VIF altos incluso si esas características no están correlacionadas con otras variables en el modelo. [8].
En nuestro caso, las características VA, VB, VC son todas altamente colineales. Las tablas cruzadas confirman esto, y una matriz de correlación de Pearson también lo haría si fueran variables continuas.
El consenso general para resolver esto es eliminar sistemáticamente una función a la vez, revisar todos los VIF, observar posibles mejoras y continuar hasta que todos los VIF caigan por debajo del umbral seleccionado. Se debe prestar especial atención a la pérdida de cualquier posible poder explicativo de las características, tanto en relación con el contexto empresarial como con la variable objetivo. Las pruebas estadísticas confirmatorias como chi-cuadrado y visualización pueden ayudar a decidir qué característica seleccionar entre dos opciones posibles.
Dejemos de lado VB y observemos los cambios de VIF:
model = smf.logit("target_satisfaction ~ C(age,Treatment('0-13')) + \
C(educ,Treatment('Some High School')) + \
C(emp,Treatment('Not Employed')) + \
C(gender,Treatment('Male')) + \
C(hispanic,Treatment('N')) + \
C(dept,Treatment('Medical Oncology')) + \
C(sign_status,Treatment('Activated')) + \
C(mode,Treatment('Internet')) + \
C(inf,Treatment('Missing')) + \
C(type,Treatment('Missing')) + \
C(VA,Treatment('Missing')) + \
C(VC,Treatment('Missing')) + \
C(VI_alt,Treatment(0))", data=df_a) results = model.fit(method='bfgs', maxiter=300000)
print(results.summary())
VA y VC todavía tienen VIF en el infinito. Peor aún, el modelo sigue generando errores estándar de NaN (no se muestran). Dejemos VC.
model = smf.logit("target_satisfaction ~ C(age,Treatment('0-13')) + \
C(educ,Treatment('Some High School')) + \
C(emp,Treatment('Not Employed')) + \
C(gender,Treatment('Male')) + \
C(hispanic,Treatment('N')) + \
C(dept,Treatment('Medical Oncology')) + \
C(sign_status,Treatment('Activated')) + \
C(mode,Treatment('Internet')) + \
C(inf,Treatment('Missing')) + \
C(type,Treatment('Missing')) + \
C(VA,Treatment('Missing')) + \
C(VI_alt,Treatment(0))", data=df_a) results = model.fit(method='bfgs', maxiter=300000)
print(results.summary())
Finalmente, el modelo produce errores estándar, por lo que aunque los VIF de la característica VA siguen siendo > 5, no es un problema debido a la segunda advertencia anterior: la categoría de referencia tiene una pequeña cantidad de casos en relación con los otros niveles.
Crédito adicional: supongamos que sabe absolutamente que VA, VB y VC son fundamentales para explicar la variable objetivo y los necesita en el modelo. Dadas las características adicionales, supongamos que la optimización funciona en un espacio complejo y que el “fallo en la inversión de Hesse” puede evitarse eligiendo nuevos solucionadores y puntos de partida (start_params en sklearn). También sería una opción entrenar un nuevo modelo que no asuma límites lineales.
Problema 4: Error de “cuasi-separación posiblemente completa” (y coeficientes y/o errores estándar excesivamente grandes)
Es posible que nos entusiasme ver coeficientes grandes y puntuaciones de alta precisión. Pero a menudo estas estimaciones son inverosímilmente grandes y se deben a otro problema común: la separación perfecta.
La separación perfecta (o casi perfecta) ocurre cuando una o más características están fuertemente asociadas con la variable objetivo. De hecho, pueden ser casi idénticos.
Puedo fabricar este error tomando la variable de destino satisfacción_objetivo y creando una nueva característica llamada nuevo_objetivo_satisfaccion eso es 95% similar a él.
df_a['new_target_satisfaction'] =
pd.concat([pd.DataFrame(np.where(df_a.target_satisfaction.iloc[:1000]==1,0,df_a.target_satisfaction[:1000]),columns=["target_satisfaction"]),
pd.DataFrame(df_a.target_satisfaction.iloc[1000:],columns=['target_satisfaction'])],axis=0)
Y pon nuevo_objetivo_satisfaccion en el modelo.
model = smf.logit("target_satisfaction ~ C(age,Treatment('0-13')) + \
C(educ,Treatment('Some High School')) + \
C(emp,Treatment('Not Employed')) + \
C(gender,Treatment('Male')) + \
C(hispanic,Treatment('N')) + \
C(dept,Treatment('Medical Oncology')) + \
C(sign_status,Treatment('Activated')) + \
C(mode,Treatment('Internet')) + \
C(inf,Treatment('Missing')) + \
C(type,Treatment('Missing')) + \
C(VA,Treatment('Missing')) + \
C(VI_alt,Treatment(0)) + \
C(new_target_satisfaction,Treatment(0))", data=df_a) results = model.fit(method='lbfgs', maxiter=1000000)
print(results.summary())
Se representan los coeficientes y los errores estándar, pero recibimos esta advertencia de cuasi separación:
La función tiene un coeficiente ridículamente alto y un índice de probabilidades de alrededor de 70.000.000, lo cual es inverosímil.
En el fondo, esto sucede porque las características que se separan “demasiado bien” crean una pendiente inflada, por lo tanto, un coeficiente y un error estándar grandes. Posiblemente, el modelo LR tampoco converja [10].
Esos dos puntos rojos, una vez eliminados los casos mal clasificados, en realidad habrían impedido la separación perfecta, ayudando al modelo LR a converger y a que las estimaciones del error estándar fueran más realistas.
Lo que hay que recordar sobre la separación perfecta en aprender es que puede generar un modelo que parece tener una precisión casi perfecta, en nuestro caso del 98%, pero que en realidad tiene una característica nuevo_objetivo_satisfaccion eso es casi un duplicado del objetivo satisfacción_objetivo.
categorical_features = ['educ','emp','gender','hispanic','dept','sign_status','mode','inf','type',
'VA','VI_alt','new_target_satisfaction']df1_y = df_a['target_satisfaction']
df1_x = df_a[['educ','emp','gender','hispanic','dept','sign_status','mode','inf','type',
'VA','VI_alt','new_target_satisfaction']]
# create a pipeline for all categorical features
categorical_transformer = Pipeline(steps=[
('ohe', OneHotEncoder(handle_unknown='ignore'))])
# create the overall column transformer
col_transformer = ColumnTransformer(transformers=[
('ohe', OneHotEncoder(handle_unknown='ignore'), categorical_features)],
# ('num', numeric_transformer, numeric_features)],
remainder='passthrough')
lr = Pipeline(steps = [('preprocessor', col_transformer),
('classifier', LogisticRegression(solver='lbfgs',max_iter=200000))])
#['liblinear', 'newton-cg', 'lbfgs', 'sag', 'saga']
X_train, X_test, y_train, y_test = train_test_split(df1_x, df1_y, test_size=0.2, random_state=101)
lr_model = lr.fit(X_train, y_train)
y_pred = lr_model.predict(X_test)
# to leverage threshold resetting
# THRESHOLD = 0.5
# y_pred = np.where(lr_model.predict_proba(X_test)[:,1] > THRESHOLD, 1, 0)
print(classification_report(y_test, y_pred))
La solución más común sería simplemente eliminar la función. Pero hay un número creciente de soluciones alternativas:
- Aplicar la corrección de Firth, que maximiza una función de probabilidad “penalizada”. Actualmente, modelos de estadísticas no tiene esta funcionalidad disponible pero R sí [11]
- La regresión penalizada también puede funcionar, específicamente probando combinaciones de solucionadores, penalizaciones y tasas de aprendizaje. [2]. Mantuve nuevo_objetivo_satisfaccion en el modelo y probé varias combinaciones, pero en este caso particular hubo poca diferencia
- Algunos practicantes intercambiarán manualmente algunos observaciones seleccionadas al azar en la característica problemática para que esté menos perfectamente separada del objetivo, como volver a agregar esos círculos rojos en la imagen de arriba [8, 10]. Ejecutar una tabla cruzada de esa característica con el objetivo puede ayudarle a determinar qué porcentaje de casos intercambiar. Mientras haces esto quizás te preguntes ¿Por qué? ¿Por qué estoy refactorizando una característica como esta sólo para que el modelo LR pueda aceptarla? Para ayudarle a dormir mejor, algunas investigaciones sostienen que la separación perfecta es solo un síntoma del modelo, no de nuestros datos. [8, 11]
- Finalmente, algunos practicantes contrarios en realidad no ven nada malo en coeficientes tan altos. [8]. Un odds-ratio muy alto simplemente sugiere que hay una fuerte asociación y predice casi perfectamente. Tenga en cuenta los hallazgos y déjelo así. La base de este argumento es que los coeficientes altos son una consecuencia desafortunada de la prueba de Wald y del índice de verosimilitud que requieren una evaluación de la información con una hipótesis alternativa.
Conclusión
La regresión logística es definitivamente versátil y poderosa si puede superar los desafíos de los conjuntos de datos del mundo real. Espero que esta descripción general proporcione una buena base para posibles soluciones. ¿Qué consejo te pareció más interesante? ¿Cuáles son algunos otros?