स्टैटिसमॉडल के साथ एक लॉजिस्टिक रिग्रेशन के लिए पॉइंटवाइज मानक त्रुटियां
स्रोत
आर में अनुप्रयोगों के साथ सांख्यिकीय सीखने का एक परिचय , यहां पाया गया: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$कमाई की संभावना को संदर्भित नहीं करता है> 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अनुमानित संभावना के आसपास एक अंतराल देगा।