Errori standard puntuali per una regressione logistica adatta a statsmodels
fonte
Un'introduzione all'apprendimento statistico con applicazioni in R , disponibile qui:https://faculty.marshall.usc.edu/gareth-james/ISL/ISLR%20Seventh%20Printing.pdf
Compito
Sto cercando di replicare l'esempio di una regressione logistica polinomiale sul set di dati "Salario" a pagina 267/8.
Cenni teorici
Secondo il libro, una volta fatte le previsioni, gli intervalli di confidenza possono essere calcolati in questo modo. Per un modello della 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$ matrice di covarianza $C$ e vettoriale $l_0^T=(1, x_0, x_0^2, x_0^3, x_0^4)$, l'errore standard puntuale è la radice quadrata di $\text{Var}[\hat{f}(x_0)]=l_0^TCl_0$. Quindi per ogni$x_0$ nel nostro set di dati abbiamo una trama di previsioni $\hat{f}(x_0)$ e un grafico degli intervalli di confidenza superiore e inferiore $\hat{f}(x_0)\pm(2\times \text{Var}[\hat{f}(x_0)])$.
Per una regressione logistica, può essere applicato lo stesso principio, ma la confidenza è intorno alla funzione logit di probabilità condizionale, al contrario delle previsioni che derivano direttamente dalla formula precedente.
Dati e codice di approccio / riutilizzabile
Prima di tutto, questo è il codice per generare il modello di regressione logistica e tracciare i risultati. Questo bit va bene e ho riprodotto con successo ciò che è nel 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")
Quindi ora provo a tracciare gli intervalli di confidenza. Non penso che ci sia un metodo per farlo direttamente in statsmodels (per favore correggimi se sbaglio).
Il mio problema
Il mio problema qui è nel calcolo degli errori standard puntuali e degli intervalli di confidenza. Sappiamo che i valori di risposta per il modello di regressione logistica devono essere$y\in [0, 1]$, poiché è una probabilità condizionata.
Il problema è quello per ogni $x_0$, il valore di $$\sqrt{l_0^TCl_0}$$sarà relativamente grande. Posso dimostrarlo utilizzando il primo valore di età,$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)
L'output di questo è pointwise_se = 6.14.
Dalla trama sopra, posso vedere che la previsione di $\text{Pr}(\text{Wage} > 250 | x=18)$ è vicino a zero e dall'esempio fornito nel libro posso vedere che l'intervallo di confidenza attorno a questo valore non è ampio e sicuramente non diventa negativo o maggiore di 1.
Se dovessi ottenere un intervallo di confidenza da un errore standard puntuale di $6.14$, la trama sarebbe sciocca e non una replica di quella nel libro.
La mia domanda
Cosa sto facendo di sbagliato nel calcolo dell'errore standard puntuale?
Risposte
Poiché stai eseguendo una regressione logistica e non una semplice regressione lineare, l'equazione $\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$non si riferisce alla probabilità di guadagnare> 250K, ma al logit di quella probabilità. Ciò equivale a dire che la regressione logistica è un modello lineare che utilizza logit come funzione di collegamento.
Quindi, devi definire le funzioni per convertire tra probabilità e logit (forse sono già implementate in Numpy o qualcosa del genere, ma sono abbastanza semplici da digitare):
def logit(p):
return np.log(p/(1-p))
def invlogit(x):
# inverse function of logit
return 1/(1+np.exp(-x))
Ora, dobbiamo applicare l'SE puntuale calcolato al logit delle stime puntuali e quindi riconvertirlo in probabilità:
upper_limit = invlogit(logit(y_pred)+1.96*std_err)
lower_limit = invlogit(logit(y_pred)-1.96*std_err)
Dov'è std_errun array con gli errori standard di$\hat f(x)$che hai calcolato correttamente. Quindi, upper_limite lower_limitfornirà un intervallo intorno alla probabilità stimata.