Erreurs standard ponctuelles pour un ajustement de régression logistique avec statsmodels
La source
Une introduction à l'apprentissage statistique avec des applications dans R , disponible ici:https://faculty.marshall.usc.edu/gareth-james/ISL/ISLR%20Seventh%20Printing.pdf
Tâche
J'essaye de reproduire l'exemple d'une régression logistique polynomiale sur le jeu de données «Salaire» à la page 267/8.
Aperçu de la théorie
Selon le livre, une fois que les prédictions ont été faites, les intervalles de confiance peuvent être calculés comme ça. Pour un modèle de la forme$$\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,$$ avec un $5\times 5$ matrice de covariance $C$ et vecteur $l_0^T=(1, x_0, x_0^2, x_0^3, x_0^4)$, l'erreur standard ponctuelle est la racine carrée de $\text{Var}[\hat{f}(x_0)]=l_0^TCl_0$. Donc pour chaque$x_0$ dans notre jeu de données, nous avons un graphique de prédictions $\hat{f}(x_0)$ et un graphique des intervalles de confiance supérieur et inférieur $\hat{f}(x_0)\pm(2\times \text{Var}[\hat{f}(x_0)])$.
Pour une régression logistique, le même principe peut être appliqué, mais la confiance se situe autour de la fonction logit de probabilité conditionnelle, par opposition aux prédictions qui proviennent directement de la formule ci-dessus.
Données et approche / code réutilisable
Tout d'abord, c'est le code pour générer le modèle de régression logistique et tracer les résultats. Ce bit est bien et j'ai réussi à reproduire ce qui est dans le livre:
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")
Alors maintenant, j'essaye de tracer les intervalles de confiance. Je ne pense pas qu'il existe une méthode pour faire cela directement dans statsmodels (veuillez me corriger si je me trompe).
Mon problème
Mon problème ici est dans le calcul des erreurs standard ponctuelles et des intervalles de confiance. Nous savons que les valeurs de réponse pour le modèle de régression logistique doivent être$y\in [0, 1]$, puisqu'il s'agit d'une probabilité conditionnelle.
Le problème est que pour chaque $x_0$, la valeur de $$\sqrt{l_0^TCl_0}$$va être relativement important. Je peux le démontrer en utilisant la première valeur d'âge,$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)
Le résultat de ceci est pointwise_se = 6.14
.
D'après l'intrigue ci-dessus, je peux voir que la prédiction de $\text{Pr}(\text{Wage} > 250 | x=18)$ est proche de zéro, et à partir de l'exemple fourni dans le livre, je peux voir que l'intervalle de confiance autour de cette valeur n'est pas large et ne devient certainement pas négatif ou supérieur à 1.
Si je devais obtenir un intervalle de confiance à partir d'une erreur standard ponctuelle de $6.14$, l'intrigue serait idiote, et non une réplique de celle du livre.
Ma question
Qu'est-ce que je fais mal dans mon calcul de l'erreur standard ponctuelle?
Réponses
Puisque vous effectuez une régression logistique et non une simple régression linéaire, l'équation $\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$ne fait pas référence à la probabilité de gagner> 250K, mais au logit de cette probabilité. Cela revient à dire que la régression logistique est un modèle linéaire qui utilise logit comme fonction de lien.
Donc, vous devez définir des fonctions à convertir entre probabilités et logits (peut-être qu'elles sont déjà implémentées dans Numpy ou quelque chose du genre, mais elles sont assez simples à taper):
def logit(p):
return np.log(p/(1-p))
def invlogit(x):
# inverse function of logit
return 1/(1+np.exp(-x))
Maintenant, nous devons appliquer le SE ponctuel que vous avez calculé au logit des estimations ponctuelles, puis reconvertir en probabilités:
upper_limit = invlogit(logit(y_pred)+1.96*std_err)
lower_limit = invlogit(logit(y_pred)-1.96*std_err)
Où std_err
est un tableau avec les erreurs standard de$\hat f(x)$que vous avez correctement calculé. Ensuite, upper_limit
et lower_limit
donnera un intervalle autour de la probabilité estimée.