Why is r squared lowered when adding polynomial features?

Question

I am trying to find a best fit line f(x) = ? for a random set of x,y coordinates.

Linear Regression with polynomial features works well for around 10 different polynomials but beyond 10 the r squared actually starts to drop!

If the new features are not useful to the Linear Regression I would assume that they would be given a coefficient of 0 and therefore adding features should not hurt the overall r squared.

I reproduced this problem when housing price predictions when creating a large amount of interaction features.

I have my python code below:

Create Random Data

import numpy as np
import matplotlib.pyplot as plt


def pol(x):
    return x * np.cos(x)

x = np.linspace(0, 12, 100)
rng = np.random.RandomState(1234)
rng.shuffle(x)
x = np.sort(x[:25])
y = pol(x) + np.random.randn(25)*2


plt.scatter(x, y, color='green', s=50, marker='.')

plt.show()

Regress and Check Each R Squared

from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures

for p in range(1,30):
    plot_range = [i/10 for i in range(0,120)]
    poly = PolynomialFeatures(p)
    X_fin = poly.fit_transform([[samp] for samp in x])
    X_fin_plot = poly.fit_transform([[samp] for samp in plot_range])
    reg = LinearRegression().fit(X_fin, y)

    from sklearn.metrics import mean_squared_error, r2_score
    print(p,r2_score(y, reg.predict(X_fin)))

Display Last Regression Line

plt.scatter(x, y, color='green', s=50, marker='.')
plt.plot(plot_range,reg.predict(X_fin_plot))
plt.show()

I also have two plots to compare. The first is with 10 polynomial features and the second is with 40. Notice how the second misses the majority of the first points.

Answer 1

You've got 25 points, so there is a perfect fitting polynomial of degree 24. That doesn't happen, so something is breaking in the OLS solver, but I'm not sure of what exactly or how to detect that. It's not too surprising though that you may have numerical issues when p gets large: you've got an x-value near 0.1 and others past 10; raising them to the 24th power pushes them very far apart, and probably generates many more significant digits than python is keeping around.

I've put together a demonstration:
https://github.com/bmreiniger/datascience.stackexchange/blob/master/53818.ipynb
Scaling the x-values helps, though we still don't find something visually matching the perfect polynomial fit.

See also https://stats.stackexchange.com/questions/350130/why-is-gradient-descent-so-bad-at-optimizing-polynomial-regression

Answer 2

My original answer was not correct, so here is a corrected answer:

When you use PolynomialFeatures(), you don't get the intended polynomials. Instead you get polynomials plus an interaction term:

from sklearn.preprocessing import PolynomialFeatures import numpy as 
np    z = np.array([[0, 1],
                    [2, 3],
                    [4, 5]]) 
poly = PolynomialFeatures(2)
print(poly.fit_transform(z))

Output is:

[[ 1.  0.  1.  0.  0.  1.]
 [ 1.  2.  3.  4.  6.  9.]
 [ 1.  4.  5. 16. 20. 25.]]

A raw polynomial should look like:

new_z = np.hstack((z**(i+1) for i in range(2)))
print(new_z)

Output is:

[[ 0  1  0  1]
 [ 2  3  4  9]
 [ 4  5 16 25]]

Here is a quick R implementation of your problem with raw polynomials:

x = c(0.12121212, 1.09090909, 3.27272727, 3.51515152, 4, 4.24242424,
  4.72727273, 4.84848485, 5.09090909, 6.18181818, 6.78787879, 7.15151515,
  7.39393939, 7.63636364, 8.24242424, 8.60606061, 9.09090909, 9.81818182,
  9.93939394, 10.3030303, 10.54545455, 10.66666667, 11.39393939, 11.63636364,
  11.87878788)

y = c(-2.87011136,1.77132943,-1.23698978,-3.09768628,-2.11919042,-4.11234626,
  -1.1684339, 1.34601699, -2.37623758,4.20290438, 6.16349341, 3.60661197,
  2.58898819, 3.80785471, -2.96359566, -5.672873, -9.71694313, -7.62778351,
  -8.95730409, -8.04664475, -5.18464423, -6.54562138, 3.45527603, 6.11936457,
  9.30106747)

regdata = data.frame(x,y)
colnames(regdata) <- c("x","y")

r2list = list()
r2adjlist = list()
plist = list()

for (p in seq(1:29)){
  reg = lm(y~poly(x,p, raw=T), data=regdata)
  print(paste0("Poly: ", p))
  print(paste0("  R2      ", summary(reg)$r.squared))
  print(paste0("  R2_adj. ", summary(reg)$adj.r.squared))
  r2list[[p]] <-  summary(reg)$r.squared
  r2adjlist[[p]] <- summary(reg)$adj.r.squared
  plist[[p]] <- p
}

plot(plist, r2list,xlab="Polynomial", ylab="R2")
lines(plist, r2list)

The R2 contingent on the degree of the polynomial is shown below:

So your initial intuition was (of course) correct, but your treatment of data was not correct.