How to calculate a prediction interval for a fitted statsmodels OLS model?
I was given the following model:
Thus, the expected change in output of increasing X by one unit is given by:
Let's say I assume a value of 40 for X.
How can I calculate a 95% confidence interval for the effect of increasing X in 0.25 units?
What follows is a replicable example.
# 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)
Great, so the expected change in y followed by increasing X from 40 to 41 is equal to Mg.
How can I calculate a 95% confidence interval for a marginal change in X of 0.25 units?
As a hint, I think this can be done with m1.t_test()
t_test works because the statistic Mg is linear in parameters.
Note, Mg is the derivative, i.e. a marginal change at point x.
To get a discrete change, we can multiply MG by dx = x1 - x to get a linear approximation or use the discrete change in y which is also linear in the parameters.
We can use a restriction defined by either a string or an explicit constraint matrix.
I use old fashioned string interpolation, and I added a seed before the simulation to get replicable results.
Marginal change given by derivative 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
with explicit restriction matrix:
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
==============================================================================
with t-test that value is 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
==============================================================================
Confidence interval for 0.25 change
What is the effect of increasing x from 40 to 40.25?
The change in the predicted value can be written as a function that is linear in parameters, so t_test can still be used for this.
x0 = 40
x1 = 40.25
m1.predict([0, x1, x1**2]) - m1.predict([0, x0, x0**2])
array([20.314731])
discrete change
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
==============================================================================
linear approximation using derivative at 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
==============================================================================