逻辑回归统计模型概率预测的置信区间
Confidence interval of probability prediction from logistic regression statsmodels
我正在尝试从 统计学习简介 中重新创建一个图，但我在弄清楚如何计算概率预测的置信区间时遇到了问题。具体来说，我正在尝试重新创建此图的右侧面板（figure 7.1) which is predicting the probability that wage>250 based on a degree 4 polynomial of age with associated 95% confidence intervals. The wage data is 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)
但是我不知道预测概率的置信区间是如何计算的。我曾多次考虑过引导数据以获得每个年龄段的概率分布，但我知道有一种更简单的方法超出了我的掌握范围。
我有估计系数协方差矩阵和与每个估计系数相关的标准误差。鉴于此信息，我将如何计算上图右侧面板中所示的置信区间？
谢谢！
您可以使用 delta method 找到预测概率的近似方差。即，
var(proba) = np.dot(np.dot(gradient.T, cov), gradient)
其中gradient是模型系数预测概率的导数向量，cov是系数的协方差矩阵。
Delta 方法被证明对所有最大似然估计都渐进地起作用。但是，如果你的训练样本很小，渐近方法可能效果不佳，你应该考虑 bootstrapping.
这是一个将 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()
它画出了如下漂亮的图：
对于您的示例，代码为
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()
然后会出现下图
看起来很像里面有大象的蟒蛇。
您可以将其与 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()
delta 方法的结果和 bootstrap 看起来几乎一样。
然而，这本书的作者选择了第三条路。他们使用
proba = np.exp(np.dot(x, 参数)) / (1 + np.exp(np.dot(x, 参数)))
并对线性部分计算置信区间，然后用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()
所以他们得到发散区间：
这些方法会产生如此不同的结果，因为它们假设不同的事物（预测概率和对数几率）呈正态分布。即，delta 方法假设预测概率是正常的，而在书中，log-odds 是正常的。实际上，其中none个在有限样本中是正态的，在无限样本中都收敛于正态，但它们的方差同时收敛于零。最大似然估计对重新参数化不敏感，但它们的估计分布是，这就是问题所在。
这是一种有指导意义的有效方法，可以在统计模型 Logit 的基础上计算拟合 ('mean_se') 和单个观察值 ('obs_se') 的标准误差 ('se') ().fit() object ('fit')，与 ISLR 书中的方法和 David Dale 的答案中的最后一个方法相同：
fit_mean = fit.model.exog.dot(fit.params)
fit_mean_se = ((fit.model.exog*fit.model.exog.dot(fit.cov_params())).sum(axis=1))**0.5
fit_obs_se = ( ((fit.model.endog-fit_mean).std(ddof=fit.params.shape[0]))**2 + \
fit_mean_se**2 )**0.5
A figure similar to the one in the book ISLR
阴影区域表示拟合和单个观察值的 95% 置信区间。
欢迎提出改进意见。
我正在尝试从 统计学习简介 中重新创建一个图，但我在弄清楚如何计算概率预测的置信区间时遇到了问题。具体来说，我正在尝试重新创建此图的右侧面板（figure 7.1) which is predicting the probability that wage>250 based on a degree 4 polynomial of age with associated 95% confidence intervals. The wage data is 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)
但是我不知道预测概率的置信区间是如何计算的。我曾多次考虑过引导数据以获得每个年龄段的概率分布，但我知道有一种更简单的方法超出了我的掌握范围。
我有估计系数协方差矩阵和与每个估计系数相关的标准误差。鉴于此信息，我将如何计算上图右侧面板中所示的置信区间？
谢谢！
您可以使用 delta method 找到预测概率的近似方差。即，
var(proba) = np.dot(np.dot(gradient.T, cov), gradient)
其中gradient是模型系数预测概率的导数向量，cov是系数的协方差矩阵。
Delta 方法被证明对所有最大似然估计都渐进地起作用。但是，如果你的训练样本很小，渐近方法可能效果不佳，你应该考虑 bootstrapping.
这是一个将 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()
它画出了如下漂亮的图：
对于您的示例，代码为
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()
然后会出现下图
看起来很像里面有大象的蟒蛇。
您可以将其与 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()
delta 方法的结果和 bootstrap 看起来几乎一样。
然而，这本书的作者选择了第三条路。他们使用
proba = np.exp(np.dot(x, 参数)) / (1 + np.exp(np.dot(x, 参数)))
并对线性部分计算置信区间，然后用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()
所以他们得到发散区间：
这些方法会产生如此不同的结果，因为它们假设不同的事物（预测概率和对数几率）呈正态分布。即，delta 方法假设预测概率是正常的，而在书中，log-odds 是正常的。实际上，其中none个在有限样本中是正态的，在无限样本中都收敛于正态，但它们的方差同时收敛于零。最大似然估计对重新参数化不敏感，但它们的估计分布是，这就是问题所在。
这是一种有指导意义的有效方法，可以在统计模型 Logit 的基础上计算拟合 ('mean_se') 和单个观察值 ('obs_se') 的标准误差 ('se') ().fit() object ('fit')，与 ISLR 书中的方法和 David Dale 的答案中的最后一个方法相同：
fit_mean = fit.model.exog.dot(fit.params)
fit_mean_se = ((fit.model.exog*fit.model.exog.dot(fit.cov_params())).sum(axis=1))**0.5
fit_obs_se = ( ((fit.model.endog-fit_mean).std(ddof=fit.params.shape[0]))**2 + \
fit_mean_se**2 )**0.5
A figure similar to the one in the book ISLR
阴影区域表示拟合和单个观察值的 95% 置信区间。
欢迎提出改进意见。