首页 文章

Logistic回归statsmodels的概率预测置信区间

提问于
浏览 次
6

我正在尝试从“统计学习简介”重新创建一个情节，并且我试图重新创建该图的右侧面板（figure 7.1），该面板预测工资> 250的概率基于年龄为4的多项式相关的95％置信区间 . 如果有人关心，工资数据是here .

我可以使用以下代码预测并绘制预测概率

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
from sklearn.preprocessing import PolynomialFeatures

wage = pd.read_csv('../../data/Wage.csv', index_col=0)
wage['wage250'] = 0
wage.loc[wage['wage'] > 250, 'wage250'] = 1

poly = Polynomialfeatures(degree=4)
age = poly.fit_transform(wage['age'].values.reshape(-1, 1))

logit = sm.Logit(wage['wage250'], age).fit()

age_range_poly = poly.fit_transform(np.arange(18, 81).reshape(-1, 1))

y_proba = logit.predict(age_range_poly)

plt.plot(age_range_poly[:, 1], y_proba)

但我对如何计算预测概率的置信区间感到茫然 . 我已经考虑过多次引导数据以获得每个时代的概率分布，但我知道有一种更简单的方法，这是我无法掌握的 .

我有估计的系数协方差矩阵和与每个估计系数相关的标准误差 . 如果给出这些信息，我将如何计算上图中右侧面板所示的置信区间？

谢谢！

python logistic-regression statsmodels confidence-interval

1 回答

  • 11

    您可以使用delta method查找预测概率的近似方差 . 也就是说，

    var(proba) = np.dot(np.dot(gradient.T, cov), gradient)

    其中 gradient 是模型系数的预测概率导数的向量， cov 是系数的协方差矩阵 .

    事实证明，Delta方法可以渐近地用于所有最大似然估计 . 但是，如果您有一个小的训练样本，渐近方法可能无法正常工作，您应该考虑自举 .

    以下是将delta方法应用于逻辑回归的玩具示例：

    import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt

# generate data
np.random.seed(1)
x = np.arange(100)
y = (x * 0.5 + np.random.normal(size=100,scale=10)>30)
# estimate the model
X = sm.add_constant(x)
model = sm.Logit(y, X).fit()
proba = model.predict(X) # predicted probability

# estimate confidence interval for predicted probabilities
cov = model.cov_params()
gradient = (proba * (1 - proba) * X.T).T # matrix of gradients for each observation
std_errors = np.array([np.sqrt(np.dot(np.dot(g, cov), g)) for g in gradient])
c = 1.96 # multiplier for confidence interval
upper = np.maximum(0, np.minimum(1, proba + std_errors * c))
lower = np.maximum(0, np.minimum(1, proba - std_errors * c))

plt.plot(x, proba)
plt.plot(x, lower, color='g')
plt.plot(x, upper, color='g')
plt.show()

    它绘制了以下精美图片：
    enter image description here

    对于您的示例，代码将是

    proba = logit.predict(age_range_poly)
cov = logit.cov_params()
gradient = (proba * (1 - proba) * age_range_poly.T).T 
std_errors = np.array([np.sqrt(np.dot(np.dot(g, cov), g)) for g in gradient])
c = 1.96 
upper = np.maximum(0, np.minimum(1, proba + std_errors * c))
lower = np.maximum(0, np.minimum(1, proba - std_errors * c))

plt.plot(age_range_poly[:, 1], proba)
plt.plot(age_range_poly[:, 1], lower, color='g')
plt.plot(age_range_poly[:, 1], upper, color='g')
plt.show()

    它会给出如下图片

    enter image description here

    看起来非常像一个里面有大象的蟒蛇 .

    您可以将它与bootstrap估计值进行比较：

    preds = []
for i in range(1000):
    boot_idx = np.random.choice(len(age), replace=True, size=len(age))
    model = sm.Logit(wage['wage250'].iloc[boot_idx], age[boot_idx]).fit(disp=0)
    preds.append(model.predict(age_range_poly))
p = np.array(preds)
plt.plot(age_range_poly[:, 1], np.percentile(p, 97.5, axis=0))
plt.plot(age_range_poly[:, 1], np.percentile(p, 2.5, axis=0))
plt.show()

    enter image description here

    delta方法和bootstrap的结果看起来几乎相同 .

    然而，本书的作者走的是第三种方式 . 他们使用的事实

    proba = np.exp（np.dot（x，params））/（1 np.exp（np.dot（x，params）））

    并计算线性部分的置信区间，然后用logit函数进行变换

    xb = np.dot(age_range_poly, logit.params)
std_errors = np.array([np.sqrt(np.dot(np.dot(g, cov), g)) for g in age_range_poly])
upper_xb = xb + c * std_errors
lower_xb = xb - c * std_errors
upper = np.exp(upper_xb) / (1 + np.exp(upper_xb))
lower = np.exp(lower_xb) / (1 + np.exp(lower_xb))
plt.plot(age_range_poly[:, 1], upper)
plt.plot(age_range_poly[:, 1], lower)
plt.show()

    所以他们得到了不同的间隔：

    enter image description here

    这些方法产生了如此不同的结果，因为它们假定正常分布的不同事物（预测概率和对数概率） . 也就是说，delta方法假设预测的概率是正常的，并且在书中，对数概率是正常的 . 事实上，它们在有限样本中都不是正常的，但它们都收敛于无限样本中，但它们的方差同时收敛于零 . 最大似然估计对重新参数化不敏感，但它们的估计分布是，这就是问题所在 .

    回复于

相关问题