# The confusion occurs due to the two different forms of statsmodels predict() method.
# This is just a consequence of the way the statsmodels folks designed the api.
# Both forms of the predict() method demonstrated and explained below.
import statsmodels.api as sm
import numpy as np
import matplotlib.pyplot as plt
# random normal x, y with y = x + 10
x = np.random.randn(100)
y = x + np.random.randn(100) + 10
fig, ax = plt.subplots(figsize=(8, 4))
ax.scatter(x, y, alpha=0.6, color='blue')
# will add the regression line to the plot and display below.
# Plot a linear regression line through the points in the scatter plot, above.
# Using statsmodels.api.OLS(Y, X).fit().
# To include a regression constant, one must use sm.add_constant() to add a column of '1s'
# to the X matrix. Basically, this tells statsmodels to calculate a constant for the regression line.
#
# FYI, the sklearn.linear_model.LinearRegression model includes a fit_intercept parameter
# and does not require the X matrix to have a column of ones.
x_matrix = sm.add_constant(x)
model = sm.OLS(y, x_matrix)
# regression_results is an object: statsmodels.regression.linear_model.RegressionResults.
regression_results = model.fit()
# nice summary
# print(regression_results.summary())
#
# There are two forms of the predict() method:
# There is sm.OLS(y, x).predict(). This predict() method has no knowledge of the fitted coefficients
# produced by model.fit(), above. This is just the way the scipy developers decided to implement
# the linear model.
# The fitted coefficients for the linear model are in RegressionResults and RegressionResults
# has its own predict().
# If you use model.predict(), you need to pass the coefficients. This form of the predict() method
# basically calculates y = ax + b where you pass the coefficients, a and b. 988
This is why the error
# occurs.
#
# Use RegressionResults.predict() which acts as you expect, except that you still must add
# a column of ones to the x-values used to predict y. This is because the original model
# was fit with a regression constant (intercept). Generally linear models are fit with an intercept
# unless there is some problem-specific reason not to.
#
x_pred = np.linspace(x.min(), x.max(), 50)
# put the X matrix in 'standard' form, i.e. with a column of ones.
x_matrix = sm.add_constant(x_pred)
y_pred = regression_results.predict(x_matrix)
# Line from RegressionResults.predict() in 'black'.
ax.plot(x_pred, y_pred, color='black')
# one more graphic to add before display below.
# So why is model.predict() provided? The calculation for a linear model is a trivial
# linear numpy calculation. For more complex models, this will not be the case
# and model.predict() can be useful.
# Here is how to use it.
coef = regression_results.params[1] # get the fitted model coefficients
const = regression_results.params[0]
# These two lines of code produce the same results, as expected.
# y_pred = coef * x_pred + const
y_pred = model.predict(params=[const, coef], exog=x_matrix)
# Line from model.predict() in 'purple' overlays the 'black' line from above.
ax.plot(x_pred, y_pred, color='purple')
plt.show()
plt.close() 988