python - 如何绘制 statsmodels 拟合的置信区间？

Question

我想在我为数据的三次样条绘制的图上显示置信区间，但我不知道应该如何完成。从理论上讲，我知道当我们接近边缘时 CI 应该偏离拟合线，但我想出的唯一解决方案是这个 janky 加法，它没有显示正确的 CI。
这是代码:

import pandas as pd
from patsy import dmatrix
import statsmodels.api as sm
import matplotlib.pyplot as plt
import numpy as np

plt.figure(figsize=(7,5))
df = pd.read_csv('http://web.stanford.edu/~oleg2/hse/wage/wage.csv').sort_values(by=['age'])
ind_df = df[['wage', 'age']].copy()

def get_bse(bse, k, m, ma, x):
  prev, ans = m, []
  k.append(ma+1)
  for i, k_ in enumerate(k):
    ans += [bse[i]]*np.sum( ((x >= prev) & (x < k_)) ); prev = k_
  return np.array(ans)


plt.scatter(df.age, df.wage, color='none', edgecolor='silver', s=10)
plt.xlabel('Age', fontsize=15)
plt.ylabel('Wage', fontsize=15)
plt.ylim((0,333))

d = 4
knots = [df.age.quantile(0.25), df.age.quantile(0.5), df.age.quantile(0.75)]

my_spline_transformation = f"bs(train, knots={knots}, degree={d}, include_intercept=True)"

transformed = dmatrix( my_spline_transformation, {"train": df.age}, return_type='dataframe' )

ft = sm.GLS(df.wage, transformed).fit()

lft = sm.Logit( (df.age > 250), transformed )
y_grid1 = lft.predict(transformed.transpose())
y_grid = ft.predict(transformed)
plt.plot(df.age, y_grid, color='crimson', linewidth=2)
plt.plot(df.age, y_grid + get_bse(ft.bse, knots, df.age.min(), df.age.max(), df.age), color='crimson', linewidth=2, linestyle='--')
plt.plot(df.age, y_grid - get_bse(ft.bse, knots, df.age.min(), df.age.max(), df.age), color='crimson', linewidth=2, linestyle='--')

plt.show()

请注意，该图应该是具有四个自由度的自然三次样条函数，但我不确定我的解决方案是否正确。实现它的正确方法是什么？

Answer 1

如果您使用:

predictions = ft.get_prediction()
df_predictions = predictions.summary_frame()
df_predictions.index = df.age.values

您将拥有一个带有 CI 拟合结果的 DataFrame。然后，将绘图更改为:

y_grid1 = lft.predict(transformed.transpose())
y_grid = ft.predict(transformed)
plt.plot(df.age, y_grid, color='crimson', linewidth=2)
#plt.plot(df.age, y_grid + get_bse(ft.bse, knots, df.age.min(), df.age.max(), df.age), color='crimson', linewidth=2, linestyle='--')
#plt.plot(df.age, y_grid - get_bse(ft.bse, knots, df.age.min(), df.age.max(), df.age), color='crimson', linewidth=2, linestyle='--')
predictions = ft.get_prediction()
df_predictions = predictions.summary_frame()
df_predictions.index = df.age.values
plt.plot(df_predictions['mean'], color='crimson')
plt.fill_between(df_predictions.index, df_predictions.mean_ci_lower, df_predictions.mean_ci_upper, alpha=.1, color='crimson')
plt.fill_between(df_predictions.index, df_predictions.obs_ci_lower, df_predictions.obs_ci_upper, alpha=.1, color='crimson')
plt.show()

产量:

其中深红色区域是平均 CI，浅红色区域是数据集中观测值的 CI。