In the previous sections, we assumed a linear relationship between explanatory and response variables. One way to account for the violation of linearity assumption is to use a polynomial regression model by adding polynomial terms:
Here, d denotes the degree of the polynomial. Although we can use polynomial regression to model a nonlinear relationship, it is still considered a multiple linear regression model because of the linear regression coefficients w. In the following subsections, we will see how we can add such polynomial terms to an existing dataset conveniently and fit a polynomial regression model.
We will now learn how to use the PolynomialFeatures transformer class from scikit-learn to add a quadratic term (d = 2) to a simple regression problem with one explanatory variable. Then, we compare the polynomial to the linear fit following these steps:
- Add a second degree polynomial term:
from sklearn.preprocessing import PolynomialFeatures >>> X = np.array([ 258.0, 270.0, 294.0, 320.0, 342.0, ... 368.0, 396.0, 446.0, 480.0, 586.0]) ... [:, np.newaxis] >>> y = np.array([ 236.4, 234.4, 252.8, 298.6, 314.2, ... 342.2, 360.8, 368.0, 391.2, 390.8]) >>> lr = LinearRegression() >>> pr = LinearRegression() >>> quadratic = PolynomialFeatures(degree=2) >>> X_quad = quadratic.fit_transform(X)
- Fit a simple linear regression model for comparison:
>>> lr.fit(X, y) >>> X_fit = np.arange(250,600,10)[:, np.newaxis] >>> y_lin_fit = lr.predict(X_fit)
- Fit a multiple regression model on the transformed features for polynomial regression:
>>> pr.fit(X_quad, y) >>> y_quad_fit = pr.predict(quadratic.fit_transform(X_fit))
- Plot the results:
>>> plt.scatter(X, y, label='training points') >>> plt.plot(X_fit, y_lin_fit, ... label='linear fit', linestyle='--') >>> plt.plot(X_fit, y_quad_fit, ... label='quadratic fit') >>> plt.legend(loc='upper left') >>> plt.show()
In the resulting plot, we can see that the polynomial fit captures the relationship between the response and explanatory variable much better than the linear fit:
>>> y_lin_pred = lr.predict(X)
>>> y_quad_pred = pr.predict(X_quad)
>>> print('Training MSE linear: %.3f, quadratic: %.3f' % (
... mean_squared_error(y, y_lin_pred),
... mean_squared_error(y, y_quad_pred)))
Training MSE linear: 569.780, quadratic: 61.330
>>> print('Training R^2 linear: %.3f, quadratic: %.3f' % (
... r2_score(y, y_lin_pred),
... r2_score(y, y_quad_pred)))
Training R^2 linear: 0.832, quadratic: 0.982As we can see after executing the code, the MSE decreased from 570 (linear fit) to 61 (quadratic fit); also, the coefficient of determination reflects a closer fit of the quadratic model (
) as opposed to the linear fit (
) in this particular toy problem.
After we learned how to construct polynomial features to fit nonlinear relationships in a toy problem, let's now take a look at a more concrete example and apply those concepts to the data in the Housing dataset. By executing the following code, we will model the relationship between house prices and LSTAT (percent lower status of the population) as using second degree (quadratic) and third degree (cubic) polynomials and compare it to a linear fit:
>>> X = df[['LSTAT']].values
>>> y = df['MEDV'].values
>>> regr = LinearRegression()
# create quadratic features
>>> quadratic = PolynomialFeatures(degree=2)
>>> cubic = PolynomialFeatures(degree=3)
>>> X_quad = quadratic.fit_transform(X)
>>> X_cubic = cubic.fit_transform(X)
# fit features
>>> X_fit = np.arange(X.min(), X.max(), 1)[:, np.newaxis]
>>> regr = regr.fit(X, y)
>>> y_lin_fit = regr.predict(X_fit)
>>> linear_r2 = r2_score(y, regr.predict(X))
>>> regr = regr.fit(X_quad, y)
>>> y_quad_fit = regr.predict(quadratic.fit_transform(X_fit))
>>> quadratic_r2 = r2_score(y, regr.predict(X_quad))
>>> regr = regr.fit(X_cubic, y)
>>> y_cubic_fit = regr.predict(cubic.fit_transform(X_fit))
>>> cubic_r2 = r2_score(y, regr.predict(X_cubic))
# plot results
>>> plt.scatter(X, y, label='training points', color='lightgray')
>>> plt.plot(X_fit, y_lin_fit,
... label='linear (d=1), $R^2=%.2f$' % linear_r2,
... color='blue',
... lw=2,
... linestyle=':')
>>> plt.plot(X_fit, y_quad_fit,
... label='quadratic (d=2), $R^2=%.2f$' % quadratic_r2,
... color='red',
... lw=2,
... linestyle='-')
>>> plt.plot(X_fit, y_cubic_fit,
... label='cubic (d=3), $R^2=%.2f$' % cubic_r2,
... color='green',
... lw=2,
... linestyle='--')
>>> plt.xlabel('% lower status of the population [LSTAT]')
>>> plt.ylabel('Price in $1000s [MEDV]')
>>> plt.legend(loc='upper right')
>>> plt.show()The resulting plot is as follows:
As we can see, the cubic fit captures the relationship between house prices and LSTAT better than the linear and quadratic fit. However, we should be aware that adding more and more polynomial features increases the complexity of a model and therefore increases the chance of overfitting. Thus, in practice it is always recommended to evaluate the performance of the model on a separate test dataset to estimate the generalization performance.
In addition, polynomial features are not always the best choice for modeling nonlinear relationships. For example, with some experience or intuition, just looking at the MEDV-LSTAT scatterplot may lead to the hypothesis that a log-transformation of the LSTAT feature variable and the square root of MEDV may project the data onto a linear feature space suitable for a linear regression fit. For instance, my perception is that this relationship between the two variables looks quite similar to an exponential function:
Since the natural logarithm of an exponential function is a straight line, I assume that such a log-transformation can be usefully applied here:
Let's test this hypothesis by executing the following code:
# transform features
>>> X_log = np.log(X)
>>> y_sqrt = np.sqrt(y)
# fit features
>>> X_fit = np.arange(X_log.min()-1,
... X_log.max()+1, 1)[:, np.newaxis]
>>> regr = regr.fit(X_log, y_sqrt)
>>> y_lin_fit = regr.predict(X_fit)
>>> linear_r2 = r2_score(y_sqrt, regr.predict(X_log))
# plot results
>>> plt.scatter(X_log, y_sqrt,
... label='training points',
... color='lightgray')
>>> plt.plot(X_fit, y_lin_fit,
... label='linear (d=1), $R^2=%.2f$' % linear_r2,
... color='blue',
... lw=2)
>>> plt.xlabel('log(% lower status of the population [LSTAT])')
>>> plt.ylabel('$sqrt{Price ; in ; $1000s ; [MEDV]}$')
>>> plt.legend(loc='lower left')
>>> plt.show()After transforming the explanatory onto the log space and taking the square root of the target variables, we were able to capture the relationship between the two variables with a linear regression line that seems to fit the data better (
) than any of the polynomial feature transformations previously:
