Точечные стандартные ошибки для логистической регрессии, соответствующие моделям статистики

Aug 17 2020

Источник

Введение в статистическое обучение с приложениями на 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$, сюжет получился бы глупым, а не копией того, что в книге.

Мой вопрос

Что я делаю неправильно при вычислении стандартной точечной ошибки?

Ответы

4 PedroSebe Aug 17 2020 at 23:35

Поскольку вы выполняете логистическую регрессию, а не простую линейную регрессию, уравнение $\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даст интервал около предполагаемой вероятности.