# Why is r squared lowered when adding polynomial features?

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.  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.

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()
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.

• "But at some point, adding more polys simply does not increase fit (even with a lot of data) because the function does not properly fetch the „data generating process“ (too much zig-zag, very „local“ fit, see right part of your second plot)." I do see that it is minimizing the squared error at the end of the plot, but I thought that the OLS would find the minimum of the residual function. Even with the higher level polynomials, the minimum of the cost function should not increase, as you can just set the new polynomial features' coefficients to 0 (Even without the help of lasso). – NS0 Jun 17 '19 at 0:16
• @NS0: You are right... I updated my answer. – Peter Jun 17 '19 at 12:27
• I'm not sure this is right still. OP's input is one-dimensional, so there are no interaction terms to speak of (and they should be included in a multivariate polynomial anyway). With just 25 points, there should be a perfectly fitting polynomial of degree 24. That not happening suggests maybe the OLS solver (in OP's python, but not your R) is breaking? – Ben Reiniger Jun 19 '19 at 19:43
• I think it must be a numerical computation issue rather than an actual error, but I'm not sure. The interaction terms in your example are the products $x_1 x_2$, but OP has only the one variable. Indeed, X_fin involves just the first p powers of the sample x-values (and has shape (25,p)). – Ben Reiniger Jun 19 '19 at 21:10
• I was thinking the same thing as @BenReiniger – NS0 Jun 24 '19 at 15:42