Errores estándar puntuales para un ajuste de regresión logística con modelos de estadísticas

Aug 17 2020

Fuente

Una introducción al aprendizaje estadístico con aplicaciones en R , que se encuentra aquí:https://faculty.marshall.usc.edu/gareth-james/ISL/ISLR%20Seventh%20Printing.pdf

Tarea

Estoy tratando de replicar el ejemplo de una regresión logística polinomial en el conjunto de datos "Salario" en la página 267/8.

Esquema de la teoría

Según el libro, una vez realizadas las predicciones, los intervalos de confianza se pueden calcular así. Para un modelo de la forma$$\hat{f}(x_0)=\hat{\beta_0}+\hat{\beta_1}x_0+\hat{\beta_2}x_0^2+\hat{\beta_3}x_0^3+\hat{\beta_4}x_0^4,$$ con un $5\times 5$ Matriz de covarianza $C$ y vector $l_0^T=(1, x_0, x_0^2, x_0^3, x_0^4)$, el error estándar puntual es la raíz cuadrada de $\text{Var}[\hat{f}(x_0)]=l_0^TCl_0$. Entonces para cada$x_0$ en nuestro conjunto de datos tenemos una gráfica de predicciones $\hat{f}(x_0)$ y una gráfica de los intervalos de confianza superior e inferior $\hat{f}(x_0)\pm(2\times \text{Var}[\hat{f}(x_0)])$.

Para una regresión logística, se puede aplicar el mismo principio, pero la confianza está alrededor de la función logit de probabilidad condicional, a diferencia de las predicciones que provienen directamente de la fórmula anterior.

Datos y enfoque / código reutilizable

En primer lugar, este es el código para generar el modelo de regresión logística y graficar los resultados. Este bit está bien y he reproducido con éxito lo que está en el libro:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from statsmodels.datasets import get_rdataset
from statsmodels.discrete import discrete_model

from sklearn.preprocessing import PolynomialFeatures
polynomial_feat = PolynomialFeatures(degree=4)

# Get dataset from the R package
data = get_rdataset("Wage", package="ISLR") 
df = data.data.reset_index() 

# Split data into wage (response, y) and age (predictor, X_orig)
y = df.wage
X_orig = df.filter(['age'], axis=1)
# Get the polynomial features from the predictor variable
X = polynomial_feat.fit_transform(X_orig) 

# Set up the test ages for a smooth results plot
X_test = np.linspace(18, 80, 1000)
X_test = X_test[:,np.newaxis] 
X_test_poly = polynomial_feat.fit_transform(X_test) 

# Create a dummy response variable, 1 if wage > 250k and 0 otherwise
y_dummy = pd.DataFrame({'wage': y[:]}) 
y_dummy['wage_split'] = np.where(y_dummy['wage'] > 250, 1, 0) 
y_dummy = y_dummy.drop(['wage'], axis=1)

# Fit a logistic regression model with statsmodels
logit_model = discrete_model.Logit(y_dummy, X).fit() 
# Get predictions, i.e. Pr(Wage > 250 | Age)
y_preds = logit_model.predict(X_test_poly)

# Plot the results
plt.figure(figsize=(8, 8)) 
plt.plot(X_test, y_preds, 'b-') 
plt.ylim(top=0.2) 
plt.xlabel("Age")
plt.ylabel("P(Wage > 250 | Age)")
plt.title("Probability of Earning > 250k with Logistic Regression")

Así que ahora intento trazar los intervalos de confianza. No creo que haya un método para hacer esto directamente en statsmodels (corríjame si me equivoco).

Mi problema

Mi problema aquí es el cálculo de los errores estándar puntuales y los intervalos de confianza. Sabemos que los valores de respuesta para el modelo de regresión logística deben ser$y\in [0, 1]$, ya que es una probabilidad condicional.

El problema es que para cada $x_0$, El valor de $$\sqrt{l_0^TCl_0}$$va a ser relativamente grande. Puedo demostrar esto usando el primer valor de edad,$x_0=18$:

# Get the covariance matrix from the model class
C = logit_model.normalized_cov_params
x = 18.
L_T = np.array([1, x, x**2, x**3, x**4])

# Compute the pointwise standard error, as outlined above
L_T = np.matrix(L_T)
L = np.transpose(L_T)
C = np.matrix(C)

var_f = np.matmul(np.matmul(L_T, C), L)
var_f = np.asarray(var_f)[0][0]
pointwise_se = np.sqrt(var_f) 
print(pointwise_se)

El resultado de esto es pointwise_se = 6.14.

En la trama anterior, puedo ver que la predicción de $\text{Pr}(\text{Wage} > 250 | x=18)$ está cerca de cero, y del ejemplo proporcionado en el libro puedo ver que el intervalo de confianza alrededor de este valor no es amplio y definitivamente no es negativo o mayor que 1.

Si tuviera que obtener un intervalo de confianza a partir de un error estándar puntual de $6.14$, la trama sería una tontería y no una réplica de eso en el libro.

Mi pregunta

¿Qué estoy haciendo mal en mi cálculo del error estándar puntual?

Respuestas

4 PedroSebe Aug 17 2020 at 23:35

Como está haciendo regresión logística y no regresión lineal simple, la ecuación $\hat f(x_0)=\hat\beta_0+\hat\beta_1x_0+\hat\beta_2x_0^2+\hat\beta_3x_0^3+\hat\beta_4x_0^4$no se refiere a la probabilidad de ganar> 250K, sino al logit de esa probabilidad. Esto es lo mismo que decir que la regresión logística es un modelo lineal que usa logit como función de enlace.

Entonces, debe definir funciones para convertir entre probabilidades y logits (tal vez ya estén implementadas en Numpy o algo así, pero son lo suficientemente simples de escribir):

def logit(p):
    return np.log(p/(1-p))

def invlogit(x):
    # inverse function of logit
    return 1/(1+np.exp(-x))

Ahora, tenemos que aplicar el SE puntual que calculó al logit de las estimaciones puntuales, y luego convertirlo de nuevo a probabilidades:

upper_limit = invlogit(logit(y_pred)+1.96*std_err)
lower_limit = invlogit(logit(y_pred)-1.96*std_err)

¿Dónde std_errhay una matriz con los errores estándar de$\hat f(x)$que calculó correctamente. Luego, upper_limity lower_limitdará un intervalo alrededor de la probabilidad estimada.