Точечные стандартные ошибки для логистической регрессии, соответствующие моделям статистики
Источник
Введение в статистическое обучение с приложениями на R можно найти здесь:https://faculty.marshall.usc.edu/gareth-james/ISL/ISLR%20Seventh%20Printing.pdf
Задача
Я пытаюсь воспроизвести пример полиномиальной логистической регрессии в наборе данных «Заработная плата» на странице 267/8.
Набросок теории
Согласно книге, как только предсказания сделаны, доверительные интервалы могут быть рассчитаны таким образом. Для модели вида$$\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,$$ с $5\times 5$ ковариационная матрица $C$ и вектор $l_0^T=(1, x_0, x_0^2, x_0^3, x_0^4)$поточечная стандартная ошибка - это квадратный корень из $\text{Var}[\hat{f}(x_0)]=l_0^TCl_0$. Так что для каждого$x_0$ в нашем наборе данных есть график прогнозов $\hat{f}(x_0)$ и график верхнего и нижнего доверительных интервалов $\hat{f}(x_0)\pm(2\times \text{Var}[\hat{f}(x_0)])$.
Для логистической регрессии может быть применен тот же принцип, но достоверность зависит от функции логита условной вероятности, в отличие от прогнозов, которые исходят непосредственно из формулы выше.
Данные и подход / повторно используемый код
Прежде всего, это код для создания модели логистической регрессии и построения графика результатов. Это нормально, и я успешно воспроизвел то, что написано в книге:
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")

Итак, теперь я пытаюсь построить доверительные интервалы. Я не думаю, что есть способ сделать это прямо в статистических моделях (пожалуйста, поправьте меня, если я ошибаюсь).
Моя проблема
Моя проблема здесь в вычислении точечных стандартных ошибок и доверительных интервалов. Мы знаем, что значения отклика для модели логистической регрессии должны быть$y\in [0, 1]$, поскольку это условная вероятность.
Проблема в том, что для каждого $x_0$, значение $$\sqrt{l_0^TCl_0}$$будет относительно большим. Я могу продемонстрировать это, используя первое значение возраста,$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)
Результатом этого является pointwise_se = 6.14
.
Из графика выше я вижу, что предсказание $\text{Pr}(\text{Wage} > 250 | x=18)$ близко к нулю, и из примера, приведенного в книге, я вижу, что доверительный интервал вокруг этого значения невелик и определенно не становится отрицательным или больше 1.
Если бы я получил доверительный интервал от точечной стандартной ошибки $6.14$, сюжет получился бы глупым, а не копией того, что в книге.
Мой вопрос
Что я делаю неправильно при вычислении стандартной точечной ошибки?
Ответы
Поскольку вы выполняете логистическую регрессию, а не простую линейную регрессию, уравнение $\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$не относятся к вероятности зарабатывания> 250К, но к логит этой вероятности. Это то же самое, что сказать, что логистическая регрессия - это линейная модель, которая использует логит как функцию связи.
Итак, вам нужно определить функции для преобразования между вероятностями и логитами (возможно, они уже реализованы в Numpy или чем-то еще, но их достаточно просто ввести):
def logit(p):
return np.log(p/(1-p))
def invlogit(x):
# inverse function of logit
return 1/(1+np.exp(-x))
Теперь мы должны применить вычисленную вами точечную SE к логиту точечных оценок, а затем преобразовать обратно в вероятности:
upper_limit = invlogit(logit(y_pred)+1.96*std_err)
lower_limit = invlogit(logit(y_pred)-1.96*std_err)
Где std_err
массив со стандартными ошибками$\hat f(x)$что вы правильно рассчитали. Затем, upper_limit
и lower_limit
даст интервал около предполагаемой вероятности.