Kesalahan standar menunjuk untuk kesesuaian regresi logistik dengan statsmodels
Sumber
Pengantar Pembelajaran Statistik dengan Aplikasi di R , ditemukan di sini:https://faculty.marshall.usc.edu/gareth-james/ISL/ISLR%20Seventh%20Printing.pdf
Tugas
Saya mencoba mereplikasi contoh regresi logistik polinomial pada kumpulan data "Upah" di halaman 267/8.
Garis besar teori
Menurut buku tersebut, setelah prediksi dibuat, interval kepercayaan dapat dihitung seperti itu. Untuk model formulir$$\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,$$ dengan $5\times 5$ matriks kovarians $C$ dan vektor $l_0^T=(1, x_0, x_0^2, x_0^3, x_0^4)$, kesalahan standar pointwise adalah akar kuadrat dari $\text{Var}[\hat{f}(x_0)]=l_0^TCl_0$. Jadi untuk setiap$x_0$ dalam kumpulan data kami, kami memiliki plot prediksi $\hat{f}(x_0)$ dan plot interval kepercayaan atas dan bawah $\hat{f}(x_0)\pm(2\times \text{Var}[\hat{f}(x_0)])$.
Untuk regresi logistik, prinsip yang sama dapat diterapkan, tetapi keyakinannya ada di sekitar fungsi logit probabilitas bersyarat, sebagai lawan dari prediksi yang berasal langsung dari rumus di atas.
Data dan pendekatan / kode yang dapat digunakan kembali
Pertama-tama, ini adalah kode untuk menghasilkan model regresi logistik dan memplot hasilnya. Sedikit ini bagus dan saya berhasil mereproduksi apa yang ada di dalam buku:
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")

Jadi sekarang saya mencoba memplot interval kepercayaan. Saya tidak berpikir ada metode untuk melakukan hal ini secara langsung di statsmodels (perbaiki saya jika saya salah).
Masalah saya
Masalah saya di sini adalah dalam perhitungan kesalahan standar pointwise dan interval kepercayaan. Kita tahu bahwa nilai respon untuk model regresi logistik haruslah$y\in [0, 1]$, karena ini adalah probabilitas bersyarat.
Masalahnya adalah untuk setiap $x_0$, nilai dari $$\sqrt{l_0^TCl_0}$$akan menjadi relatif besar. Saya dapat mendemonstrasikan ini dengan menggunakan nilai usia pertama,$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)
Outputnya adalah pointwise_se = 6.14
.
Dari plot di atas, saya dapat melihat prediksi dari $\text{Pr}(\text{Wage} > 250 | x=18)$ mendekati nol, dan dari contoh yang diberikan di buku saya dapat melihat bahwa interval kepercayaan di sekitar nilai ini tidak lebar, dan pasti tidak negatif atau lebih besar dari 1.
Jika saya mendapatkan interval kepercayaan dari kesalahan standar pointwise dari $6.14$, plotnya akan konyol, dan bukan replikasi di buku.
Pertanyaan saya
Apa yang saya lakukan salah dalam perhitungan saya dari kesalahan standar pointwise?
Jawaban
Karena Anda melakukan regresi logistik dan bukan regresi linier sederhana, persamaannya $\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$tidak mengacu pada probabilitas menghasilkan> 250K, tetapi pada logit dari probabilitas itu. Ini sama dengan mengatakan bahwa regresi logistik adalah model linier yang menggunakan logit sebagai fungsi tautan.
Jadi, Anda harus menentukan fungsi untuk mengonversi antara probabilitas dan logit (mungkin sudah diterapkan di Numpy atau semacamnya, tetapi cukup sederhana untuk diketik):
def logit(p):
return np.log(p/(1-p))
def invlogit(x):
# inverse function of logit
return 1/(1+np.exp(-x))
Sekarang, kita harus menerapkan SE pointwise yang Anda hitung ke logit estimasi titik, dan kemudian mengubahnya kembali ke probabilitas:
upper_limit = invlogit(logit(y_pred)+1.96*std_err)
lower_limit = invlogit(logit(y_pred)-1.96*std_err)
Di mana std_err
array dengan kesalahan standar$\hat f(x)$yang Anda hitung dengan benar. Kemudian, upper_limit
dan lower_limit
akan memberikan interval di sekitar perkiraan probabilitas.