통계 모델에 맞는 로지스틱 회귀 분석에 대한 점별 표준 오차
출처
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)$ 0에 가깝고 책에 제공된 예에서이 값 주변의 신뢰 구간이 넓지 않고 확실히 음수 또는 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$250K 이상의 수익을 올릴 확률이 아니라 그 확률 의 로짓 을 의미합니다. 이것은 로지스틱 회귀가 로짓을 연결 함수로 사용하는 선형 모델이라는 말과 동일합니다.
따라서 확률과 로짓간에 변환 할 함수를 정의해야합니다 (아마도 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추정 된 확률 주위의 간격을 줄 것이다.