Training r-squared decreases after adding higher degree polynomial terms to regression model

Question

I was playing around with some examples to get some experience using the PolyFeatures tool from Scikit-Learn, and I ran into something strange. I iteratively added higher and higher degree polynomial features to my regression model, and this would occasionally cause the r-squared value for the model to decrease, which should not be possible.

I originally noticed this while working with the Boston Housing dataset, but here is a simple example demonstrating the issue:

import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures

np.random.seed(1)
n = 500
x1 = np.random.uniform(0, 3, n)
x2 = np.random.uniform(0, 3, n)
x3 = np.random.uniform(0, 3, n)
y = 3 + 0.01 * x1**3 + 0.02 * x2**2 + 0.03 * x2*x3 + np.random.normal(0, 0.2, n)

X = np.vstack((x1, x2, x3)).transpose()

for d in range(1, 9):
    poly = PolynomialFeatures(d)
    Xp = poly.fit_transform(X)

    mod = LinearRegression()
    mod.fit(Xp, y)

    print('Degree', d, '- Training r-Squared:', mod.score(Xp, y))

The output of this code is:

Degree 1 - Training r-Squared: 0.2773006611069333
Degree 2 - Training r-Squared: 0.3168358821057937
Degree 3 - Training r-Squared: 0.33258321401873814
Degree 4 - Training r-Squared: 0.3160261669178669
Degree 5 - Training r-Squared: 0.3729512734983266
Degree 6 - Training r-Squared: 0.3234788901084178
Degree 7 - Training r-Squared: 0.24399386671590273
Degree 8 - Training r-Squared: 0.42981336522995917

Notice that r-squared drops on three occasions as the degree of the model increases (from 3-4, 5-6, and 6-7).

Any ideas why this is happening? Thanks in advance!

---

I came up with a simpler example involving simple linear regression. I ran it in Python and got the same unexpected results. Then I ran it in R, and everything looked normal. I was able to get the expected behavior in Python after simple scaling the features.

Answer 1

In fact the reason has to do with colinearity among the features, but not in the way suggested by the other answer. Rather, it has to do with the default parameter settings of the PolynomialFeatures and LinearRegression classes. If you read the documentation of the respective classes, you will see that PolynomialFeatures class has an "include_bias" attribute, which by default is set to True. This means that the first column of the returned feature matrix will be all ones. Similarly, the LinearRegression class has a "fit_intercept" attribute which is also set to True by default. This means that before the linear regression is performed, an extra column of 1s will be appended to the feature matrix. However the feature matrix already contains a column of 1s, so we will have colinearity in the form of two identical columns.

We can verify this by examining the values of the coefficients of the regression:

import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures

np.random.seed(1)
n = 500
x1 = np.random.uniform(0, 3, n)
x2 = np.random.uniform(0, 3, n)
x3 = np.random.uniform(0, 3, n)
y = 3 + 0.01 * x1**3 + 0.02 * x2**2 + 0.03 * x2*x3 + np.random.normal(0, 0.2, n)

X = np.vstack((x1, x2, x3)).transpose()

d=9
poly = PolynomialFeatures(d)
Xp = poly.fit_transform(X)

mod = LinearRegression()
mod.fit(Xp, y)
print(mod.intercept_,np.abs(mod.coef_).max()

the output is 285693288.6867925 285693277.157223,thus suggesting an issue of numerical overflow, as would be expected when there is perfect colinearity among the predictors. Now consider instead the following code, where the only difference is that we have not included an intercept term in the LinearRegression class.

import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures

np.random.seed(1)
n = 500
x1 = np.random.uniform(0, 3, n)
x2 = np.random.uniform(0, 3, n)
x3 = np.random.uniform(0, 3, n)
y = 3 + 0.01 * x1**3 + 0.02 * x2**2 + 0.03 * x2*x3 + np.random.normal(0, 0.2, n)

X = np.vstack((x1, x2, x3)).transpose()

d=9
poly = PolynomialFeatures(d)
Xp = poly.fit_transform(X)

mod = LinearRegression(fit_intercept=False)
mod.fit(Xp, y)
print(mod.intercept_,np.abs(mod.coef_).max()

The output is 0.0 247.7975362818573, so the coefficients take much more reasonable values and numerical overflow is likely not an issue (albeit, the problem is still fairly ill-conditioned).

Finally, when I run your exact code, except setting fit_intercept=False, i get the following output:

Degree 1 - Training r-Squared: 0.2773006611069332
Degree 2 - Training r-Squared: 0.31683588210579394
Degree 3 - Training r-Squared: 0.33258321401873825
Degree 4 - Training r-Squared: 0.3496077641849682
Degree 5 - Training r-Squared: 0.37578112218309345
Degree 6 - Training r-Squared: 0.40396955460657147
Degree 7 - Training r-Squared: 0.4387038406549014
Degree 8 - Training r-Squared: 0.4998373846906773

The r squared values now increase monotonically, as they should.

Answer 2

Although, adding new features increase the R-squared value in general. But adding higher polynomial features can induce multi collinearity into the model. This makes the parameter coefficients highly unstable tilting from positive to negative. Since R-squared = ((Total Sum of squares before regression) - (Residual sum of squares))/Residual sum of squares. In this case, the variation in RSS might reduce the R-squared value (from 3-4, 5-6 etc). Though my answer is an assumption,we can calculate the TSS and RSS from the data and verify wherever we observe declining in the R-squared.