如何计算拟合 statsmodels OLS 模型的预测区间？

Question

我得到了以下模型：

因此，X增加一个单位的预期输出变化为：

假设我假设 X 的值为 40。

如何计算以 0.25 个单位 增加 X 的效果的 95% 置信区间？

下面是一个可复制的例子。

# Generate data
import pandas as pd
from scipy import stats as st
df = pd.DataFrame({'const':1,'X':st.norm(loc=40, scale=5).rvs(1000)})
df['X_sq'] = df['X'].pow(2)
df['y'] = 1200 + df['X'] + df['X_sq'] + st.norm().rvs(1000)
df = df[['y','const','X','X_sq']]

# Declare and fit model
y = df['y']
X = df[['const','X','X_sq']]
m1 = OLS(endog=y, exog=X).fit()

# Assume a value for `Xi`
x = 40

# Predicted marginal effect of increasing `Xi` in ONE UNIT
Mg = m1.params['X'] + (2 * m1.params['X_sq'] * x)

很好，所以 y 的预期变化是 X 从 40 增加到 41 等于 Mg。

如何计算 0.25 个单位的 X 边际变化的 95% 置信区间？

作为提示，我认为这可以通过 m1.t_test()

来完成

Answer 1

t_test 起作用是因为统计数据 Mg 在参数方面是线性的。注意，Mg 是导数，即 x 点的边际变化。要获得离散变化，我们可以将 MG 乘以 dx = x1 - x 以获得线性近似值或使用 y 的离散变化，这在参数中也是线性的。

我们可以使用由字符串或显式约束矩阵定义的限制条件。

我使用老式的字符串插值，并在模拟之前添加了一个种子以获得可复制的结果。

导数 Mg 给出的边际变化

np.random.seed(987125348)

"X + %f * X_sq" % (2 * x)
'X + 80.000000 * X_sq'

m1.t_test("X + %f * X_sq" % (2 * x))
<class 'statsmodels.stats.contrast.ContrastResults'>
                             Test for Constraints                             
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
c0            81.0089      0.007   1.24e+04      0.000      80.996      81.022
==============================================================================

Mg
81.00891173785777

具有显式限制矩阵：

m1.t_test([0, 1, 2 * x])
<class 'statsmodels.stats.contrast.ContrastResults'>
                             Test for Constraints                             
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
c0            81.0089      0.007   1.24e+04      0.000      80.996      81.022
==============================================================================

t 检验值为 80

m1.t_test("X + %f * X_sq = 80" % (2 * x))
<class 'statsmodels.stats.contrast.ContrastResults'>
                             Test for Constraints                             
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
c0            81.0089      0.007    154.100      0.000      80.996      81.022
==============================================================================

0.25 变化的置信区间

将 x 从 40 增加到 40.25 有什么影响？

预测值的变化可以写成参数为线性的函数，所以t_test仍然可以用于此。

x0 = 40
x1 = 40.25
m1.predict([0, x1, x1**2]) - m1.predict([0, x0, x0**2])
array([20.314731])

离散变化

m1.t_test([0, (x1 - x0), (x1**2 - x0**2)])
<class 'statsmodels.stats.contrast.ContrastResults'>
                             Test for Constraints                             
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
c0            20.3147      0.002   1.24e+04      0.000      20.312      20.318
==============================================================================

在 x0 = 40 处使用导数的线性近似

dx = x1 - x0
m1.t_test([0, 1 * dx, 2 * x0 * dx])
<class 'statsmodels.stats.contrast.ContrastResults'>
                             Test for Constraints                             
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
c0            20.2522      0.002   1.24e+04      0.000      20.249      20.255
==============================================================================

如何计算拟合 statsmodels OLS 模型的预测区间？

How to calculate a prediction interval for a fitted statsmodels OLS model?

python

statistics

statsmodels

导数 Mg 给出的边际变化

0.25 变化的置信区间