sklearn - scikit learn: cómo verificar la importancia de los coeficientes

sklearn adaboost

Question

Traté de hacer un LR con SKLearn para un conjunto de datos bastante grande con ~ 600 dummy y pocas variables de intervalo (y 300 K líneas en mi conjunto de datos) y la matriz de confusión resultante parece sospechosa. Quería verificar el significado de los coeficientes devueltos y ANOVA, pero no puedo encontrar cómo acceder a él. ¿Es posible? ¿Y cuál es la mejor estrategia para los datos que contienen muchas variables ficticias? ¡Muchas gracias!

Answer 1

Scikit-learn deliberadamente no es compatible con la inferencia estadística. Si desea pruebas de significación de coeficientes listos para usar (y mucho más), puede usar el estimador Logit de Statsmodels . Este paquete imita los modelos de interfaz glm en R, por lo que podría encontrarlo familiar.

Si aún desea seguir con LogisticRegression de scikit-learn, puede usar una aproximación asimétrica a la distribución de las estimaciones de likelihiood máximas. Precisamente, para un vector de estimaciones de máxima verosimilitud theta , su matriz de varianza-covarianza puede estimarse como inverse(H) , donde H es la matriz de Hesse de probabilidad logarítmica en theta . Esto es exactamente lo que hace la función siguiente:

import numpy as np
from scipy.stats import norm
from sklearn.linear_model import LogisticRegression

def logit_pvalue(model, x):
    """ Calculate z-scores for scikit-learn LogisticRegression.
    parameters:
        model: fitted sklearn.linear_model.LogisticRegression with intercept and large C
        x:     matrix on which the model was fit
    This function uses asymtptics for maximum likelihood estimates.
    """
    p = model.predict_proba(x)
    n = len(p)
    m = len(model.coef_[0]) + 1
    coefs = np.concatenate([model.intercept_, model.coef_[0]])
    x_full = np.matrix(np.insert(np.array(x), 0, 1, axis = 1))
    ans = np.zeros((m, m))
    for i in range(n):
        ans = ans + np.dot(np.transpose(x_full[i, :]), x_full[i, :]) * p[i,1] * p[i, 0]
    vcov = np.linalg.inv(np.matrix(ans))
    se = np.sqrt(np.diag(vcov))
    t =  coefs/se  
    p = (1 - norm.cdf(abs(t))) * 2
    return p

# test p-values
x = np.arange(10)[:, np.newaxis]
y = np.array([0,0,0,1,0,0,1,1,1,1])
model = LogisticRegression(C=1e30).fit(x, y)
print(logit_pvalue(model, x))

# compare with statsmodels
import statsmodels.api as sm
sm_model = sm.Logit(y, sm.add_constant(x)).fit(disp=0)
print(sm_model.pvalues)
sm_model.summary()

Las salidas de print() son idénticas, y resultan ser coeficientes p-values.

[ 0.11413093  0.08779978]
[ 0.11413093  0.08779979]

sm_model.summary() también imprime un resumen HTML muy bien formateado.