python - 如何计算拟合 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太好了，因此从 40增加到41的预期变化X等于Mg。

如何计算X0.25 个单位的边际变化的 95% 置信区间？

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

标签： pythonstatisticsstatsmodels

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仍然可以用于此。​</p>

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
==============================================================================