使用 StatsModels 绘制二阶多项式的分位数回归

Posted 2023-03-12

【中文标题】使用 StatsModels 绘制二阶多项式的分位数回归

【英文标题】：Using StatsModels to plot quantile regression for 2nd order polynomial

【发布时间】：2016-05-13 02:55:42

【问题描述】：

我正在按照 StatsModels 示例 here 绘制分位数回归线。只需对我的数据稍作修改，该示例就可以很好地生成此图（请注意，我已修改代码以仅绘制 0.05、0.25、0.5、0.75 和 0.95 分位数）：

但是，我想为二阶多项式拟合（而不是线性）绘制 OLS 拟合和相应的分位数。例如，这里是相同数据的二阶 OLS 行：

如何修改链接示例中的代码以生成非线性分位数？

这是我从链接示例修改的相关代码以生成第一个图：

d = 'temp': x, 'dens': y
df = pd.DataFrame(data=d)

# Least Absolute Deviation
# 
# The LAD model is a special case of quantile regression where q=0.5

mod = smf.quantreg('dens ~ temp', df)
res = mod.fit(q=.5)
print(res.summary())

# Prepare data for plotting
# 
# For convenience, we place the quantile regression results in a Pandas DataFrame, and the OLS results in a dictionary.

quantiles = [.05, .25, .50, .75, .95]
def fit_model(q):
    res = mod.fit(q=q)
    return [q, res.params['Intercept'], res.params['temp']] + res.conf_int().ix['temp'].tolist()

models = [fit_model(x) for x in quantiles]
models = pd.DataFrame(models, columns=['q', 'a', 'b','lb','ub'])

ols = smf.ols('dens ~ temp', df).fit()
ols_ci = ols.conf_int().ix['temp'].tolist()
ols = dict(a = ols.params['Intercept'],
           b = ols.params['temp'],
           lb = ols_ci[0],
           ub = ols_ci[1])

print(models)
print(ols)

x = np.arange(df.temp.min(), df.temp.max(), 50)
get_y = lambda a, b: a + b * x

for i in range(models.shape[0]):
    y = get_y(models.a[i], models.b[i])
    plt.plot(x, y, linestyle='dotted', color='grey')

y = get_y(ols['a'], ols['b'])
plt.plot(x, y, color='red', label='OLS')

plt.scatter(df.temp, df.dens, alpha=.2)
plt.xlim((-10, 40))
plt.ylim((0, 0.4))
plt.legend()
plt.xlabel('temp')
plt.ylabel('dens')
plt.show()

【参考方案1】：

经过一天的研究，想出了一个解决方案，所以发布了我自己的答案。非常感谢 StatsModels 的 Josef Perktold 提供的帮助。

以下是相关代码和情节：

d = 'temp': x, 'dens': y
df = pd.DataFrame(data=d)

x1 = pd.DataFrame('temp': np.linspace(df.temp.min(), df.temp.max(), 200))

poly_2 = smf.ols(formula='dens ~ 1 + temp + I(temp ** 2.0)', data=df).fit()
plt.plot(x, y, 'o', alpha=0.2)
plt.plot(x1.temp, poly_2.predict(x1), 'r-', 
         label='2nd order poly fit, $R^2$=%.2f' % poly_2.rsquared, 
         alpha=0.9)
plt.xlim((-10, 50))
plt.ylim((0, 0.25))
plt.xlabel('mean air temp')
plt.ylabel('density')
plt.legend(loc="upper left")


# with quantile regression

# Least Absolute Deviation
# The LAD model is a special case of quantile regression where q=0.5

mod = smf.quantreg('dens ~ temp + I(temp ** 2.0)', df)
res = mod.fit(q=.5)
print(res.summary())

# Quantile regression for 5 quantiles

quantiles = [.05, .25, .50, .75, .95]

# get all result instances in a list
res_all = [mod.fit(q=q) for q in quantiles]

res_ols = smf.ols('dens ~ temp + I(temp ** 2.0)', df).fit()


plt.figure()

# create x for prediction
x_p = np.linspace(df.temp.min(), df.temp.max(), 50)
df_p = pd.DataFrame('temp': x_p)

for qm, res in zip(quantiles, res_all):
    # get prediction for the model and plot
    # here we use a dict which works the same way as the df in ols
    plt.plot(x_p, res.predict('temp': x_p), linestyle='--', lw=1, 
             color='k', label='q=%.2F' % qm, zorder=2)

y_ols_predicted = res_ols.predict(df_p)
plt.plot(x_p, y_ols_predicted, color='red', zorder=1)
#plt.scatter(df.temp, df.dens, alpha=.2)
plt.plot(df.temp, df.dens, 'o', alpha=.2, zorder=0)
plt.xlim((-10, 50))
plt.ylim((0, 0.25))
#plt.legend(loc="upper center")
plt.xlabel('mean air temp')
plt.ylabel('density')
plt.title('')
plt.show()

红线：二阶多项式拟合

黑色虚线：第 5、25、50、75、95 个百分位数

【讨论】：

你能澄清一下模型中“我”的目的是什么吗？ 'dens ~ temp + I(temp ** 2.0)'

这是用于 Patsy 的语法（一个公式“迷你语言”）。见patsy.readthedocs.io/en/v0.1.0/API-reference.html#。 (x ** 2.0) 在公式中不起作用，您需要 I 所以 patsy 不会尝试进行专为分类设计的转换。

我应该担心FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result. 吗？

我不确定代码的哪一部分引发了这个警告...... StatsModels 包中的东西？但是，是的，这看起来可能值得关注。如果它在您自己的代码中，也许一些简单的修改可以解决警告。如果它在 StatsModels 中，可能需要引起 StatsModels 开发人员的注意。