Comparing RSS from linear and higher power models in a training and test data

I'm currently trying to complete the book 'Introduction to Statistical Learning' by James et al. and I'm stuck in one of the exercises, trying to get some logic out of the results.

The question is this one (see in line):

1. I collect a set of data (n = $100$ observations) containing a single predictor and a quantitative response. I then fit a linear regression model to the data, as well as a separate cubic regression, i.e. $Y = β_0 +β_1X +β_2X^2 +β_3X^3 +\varepsilon$.

(a) Suppose that the true relationship between X and Y is linear, i.e. $Y = β_0 + β_1X + \varepsilon$. Consider the training residual sum of squares (RSS) for the linear regression, and also the training RSS for the cubic regression. Would we expect one to be lower than the other, would we expect them to be the same, or is there not enough information to tell? Justify your answer.

Given that the true relationship is linear, I thought adding a cubic term wouldn't improve the fit, it would take up degrees of freedom, and simple would increase the RSS. I tried to test this empirically:

training <- function(p, n) {
set.seed(1)
x <- rnorm(n)
y <- 5 + 2*x + rnorm(n, sd = 0.5)

sigma_linear <- summary(lm(y[p] ~ x[p]))$sigma sigma_polynomial <- summary(lm(y[p] ~ x[p] + I(x[p]^2) + I(x[p]^3)))$sigma

c(linear = sigma_linear, polynomial = sigma_polynomial)
}

training(1:20, 100)
#    linear polynomial
# 0.4177607  0.4203941

If you run this, the linear RSS is smaller than the cubic one, which is what I would expect. But remove the set.seed option from the function and replicate this a couple of times and in many cases the cubic one will be higher. This might be simple random noise and the cubic term doesn't add or remove anything from the model. I don't understand why the cubic model would give a better 'fit'. In fact, in that exact link, both higher power terms are insignificant and the RSE is higher for the cubic model.

Now, using the remaining part of the data:

training(20:100, 100)
#    linear polynomial
# 0.4957765  0.4984665

The cubic sigma is still higher. But again, replicate without the set.seed option and it will sometimes be smaller or higher.

Why do then people expect for the first cubic model to have a lower RSS in the linear model?

• Are you sure that you are not just fitting noise? Oct 21 '16 at 17:51
• Hi Antoni, I think I understood correctly. Thanks for you answer. If I am right, then, adding higher order terms to the 'test' data will have the same effect, i.e. lowering the RSS, although it will be higher than in the training data, right? Oct 23 '16 at 1:43
• I think this is an issue that applies to fitting models to data far beyond machine learning, and it is one of the strong selling points for settling for the deceivingly rudimentary straight line to understand relationships in scientific research. In the particular instance of ML, the ability to fit the training set is practically unlimited, and the RSS can be decreased virtually to zero, but it will not carry out to the test set. Oct 23 '16 at 1:50
• There is a wonderful series of youtube series by Yaser Abu-Mostafa that takes you to the most interesting topics at the root of ML - the actual ability to learn, in the computer sense, from the data. Much more fascinating IMHO than any particular applications. Oct 23 '16 at 1:54

Here is the result of $\small 1,000$ simulations of your basic linear data-generating process, each time with $100$ data points, and each time fitted with polynomial regression models of degree $1$ to $5$, and getting the corresponding $5$ different SSR for each iteration:

mat = matrix(rep(0,5000), nrow=1000)
for(i in 1:1000){
n = 100
x = rnorm(n)
y = 5 + 2 * x + rnorm(n, 0.5)
for(j in 1:5){
mat[i,j] = sum(residuals(lm(y ~ poly(x,j,raw=T)))^2)
}
} Clearly, even if we use the same exact linear dataset for each iteration, the tendency is to a progressive decrease in $\text{RSS}$ with increasing polynomial degree in the model.

In fact these differences are statistically significant when running an ANOVA between the different polynomial models (p-value: 6.444e-10).

BACKGROUND:

In general, the higher the order of the polynomial model, the better the fit is going to be to the data points in the training set. You can see a funny example of it in this post. Logically, the problem is overfitting, and the consequent out-of-sample error.

\begin{align} \text{SST}_{\text{otal}} &= \color{red}{\text{SSE}_{\text{xplained}}}+\color{blue}{\text{SSR}_{\text{esidual}}}\\ \end{align} or \begin{align}\small \sum(y_i-\bar y)^2 &=\small \color{red}{\sum(\hat y_i-\bar y)^2}+\color{blue}{\sum(y_i-\hat y_i)^2} \end{align}

In the OP, the R code calls for the Residual Standard Error instead of the sum of squares residual (SSR or RSS), which amounts to:

$$\small \text{RSE} = \sqrt{\frac{\sum(y_i - \hat y_i)^2}{\text{df}}}$$

However, it doesn't change what follows.

Here is an example with the dataset mtcars plotting the predicted polynomial lines on the data points of miles per gallon (mpg) versus weight (wt) of different car models, resulting from increasing polynomial order models from $1$ to $10$: Although the relationship is fairly linear, as can be inferred from general knowledge about cars, and from inspecting the data points, as the the degree of the polynomial model increases the predicted points fall closer and closer to the actual data points because there is increased ability to accommodate individual variations resulting from noise in the data.

Concomitantly, there is a drop in the RSS: all of it at the expense of overfitting and increased out-of-sample error.