İstatistik modellerine uygun lojistik regresyon için noktasal standart hatalar
Kaynak
R'deki Uygulamalar ile İstatistik Öğrenmeye Giriş , burada bulunan:https://faculty.marshall.usc.edu/gareth-james/ISL/ISLR%20Seventh%20Printing.pdf
Görev
267/8. Sayfadaki "Ücret" veri kümesindeki bir polinom lojistik regresyon örneğini kopyalamaya çalışıyorum.
Teori taslağı
Kitaba göre tahminler yapıldıktan sonra güven aralıkları bu şekilde hesaplanabiliyor. Formun bir modeli için$$\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,$$ Birlikte $5\times 5$ kovaryans matrisi $C$ ve vektör $l_0^T=(1, x_0, x_0^2, x_0^3, x_0^4)$noktasal standart hata, kareköktür $\text{Var}[\hat{f}(x_0)]=l_0^TCl_0$. Yani her biri için$x_0$ veri setimizde bir tahmin planımız var $\hat{f}(x_0)$ ve üst ve alt güven aralıklarının bir grafiği $\hat{f}(x_0)\pm(2\times \text{Var}[\hat{f}(x_0)])$.
Lojistik regresyon için, aynı ilke uygulanabilir, ancak doğrudan yukarıdaki formülden gelen tahminlerin aksine, güven koşullu olasılık logit fonksiyonu civarındadır.
Veri ve yaklaşım / yeniden kullanılabilir kod
Her şeyden önce, bu, lojistik regresyon modelini oluşturmak ve sonuçları çizmek için kullanılan koddur. Bu bit iyi ve kitapta yazanları başarıyla yeniden oluşturdum:
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")
Şimdi güven aralıklarını çizmeye çalışıyorum. Ben yok düşünüyorum statsmodels doğrudan bunu yapmanın bir yöntemi (yanılıyorsam beni düzeltin lütfen) bulunmaktadır.
Benim sorunum
Buradaki sorun, noktasal standart hataların ve güven aralıklarının hesaplanmasıdır. Lojistik regresyon modeli için yanıt değerlerinin olması gerektiğini biliyoruz$y\in [0, 1]$, çünkü şartlı bir olasılıktır.
Sorun şu ki, herkes için $x_0$, değeri $$\sqrt{l_0^TCl_0}$$nispeten büyük olacak. Bunu ilk yaş değerini kullanarak gösterebilirim,$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)
Bunun çıktısı pointwise_se = 6.14.
Yukarıdaki arsadan, tahmininin $\text{Pr}(\text{Wage} > 250 | x=18)$ sıfıra yakın ve kitapta verilen örnekten bu değerin etrafındaki güven aralığının geniş olmadığını ve kesinlikle negatif veya 1'den büyük olmadığını görebiliyorum.
Noktasal standart hatadan bir güven aralığı alacak olsaydım $6.14$olay örgüsü aptalca olurdu ve kitapta bunun bir kopyası olmaz.
Benim sorum
Noktasal standart hata hesaplamamda neyi yanlış yapıyorum?
Yanıtlar
Basit doğrusal regresyon yerine lojistik regresyon yaptığınız için, denklem $\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$kazanç> 250K olasılık bakın, bunlarla gelmez logit o olasılık. Bu, lojistik regresyonun, logit'i bağlantı işlevi olarak kullanan doğrusal bir model olduğunu söylemekle aynıdır.
Bu nedenle, olasılıklar ve günlükler arasında dönüştürmek için işlevler tanımlamanız gerekir (belki bunlar zaten Numpy'de veya başka bir yerde uygulanmıştır, ancak yazmak için yeterince basittir):
def logit(p):
return np.log(p/(1-p))
def invlogit(x):
# inverse function of logit
return 1/(1+np.exp(-x))
Şimdi, hesapladığınız noktasal SE'yi nokta tahminlerinin logitine uygulamalı ve sonra tekrar olasılıklara dönüştürmeliyiz:
upper_limit = invlogit(logit(y_pred)+1.96*std_err)
lower_limit = invlogit(logit(y_pred)-1.96*std_err)
std_errStandart hataları olan bir dizi nerede$\hat f(x)$doğru hesapladığın. Ardından, upper_limitve lower_limittahmini olasılık etrafında bir aralık verecektir.