Why is the use of high order polynomials for regression discouraged?

Question

I've read many times on this site that high order polynomials (generally more than third) shouldn't be used in linear regression, unless there is a substantial justification to do so.

I understand the issues about extrapolation (and prediction at the boundaries).

Since extrapolation isn't important to me...

1. Are high order polynomials also a bad way of approximating the underlying function within the range of the data points? (i.e. interpolation)
2. If so, what problems are arising?

I don't mind being redirected to a good book or paper about this. Thanks.

Answer 1

I cover this in some detail in Chapter 2 of RMS. Briefly, besides extrapolation problems, ordinary polynomials have these problems:

1. The shape of the fit in one region of the data is influenced by far away points
2. Polynomials cannot fit threshold effects, e.g., a nearly flat curve that suddenly accelerates
3. Polynomials cannot fit logarithmic-looking relationships, e.g., ones that get progressively flatter over a long interval
4. Polynomials can't have a very rapid turn

These are reasons that regression splines are so popular, i.e., segmented polynomials tend to work better than unsegmented polynomials. You can also relax a continuity assumption for a spline if you want to have a discontinuous change point in the fit.

Answer 2

Yes, polynomials are also problematic in interpolation, because of overfitting and high variability.

Here is an example. Assume your dependent variable $y$ is uniformly distributed on the interval $[0,1]$. You also have a "predictor" variable $x$, also uniformly distributed on the interval $[0,1]$. However, there is no relationship between the two. Thus, any regression of $y$ on any power of $x$ will be overfitting.

Now, assume we draw 20 data points from this data generating process and fit $y\sim x^n$ for $n=1, 2, 3, 5$. Here are the fits:

As you see, the fit gets "wigglier" for higher $n$.

However, one key problem is that (of course) the fit will depend on the data we have randomly sampled from our data generating process. After all, we could have drawn 20 quite different pairs $(x,y)$. Let's repeat the exercise another three times, with new random samples each time. Below, the top row is the same as the previous plot, and the three other rows are just the fits based on new samples:

The main problem is visible when you compare the left column (linear fits) and the right column ($x^5$ fits): the fit for a lower order polynomial is much less variable and dependent on the randomness in our data sampling than the fit for the high order polynomial. If we want to interpolate $y$ for some $x$ even somewhere in the middle of the interval, using a higher order polynomial will yield a fit that is much more wobbly than a lower order polynomial.

R code:

nn <- 20
xx <- seq(0,1,by=.01)

png("sims_1.png",width=800,height=200)
    opar <- par(mfrow=c(1,4),mai=c(.3,.3,.3,.1),las=1)
    set.seed(1)
    obs <- data.frame(x=runif(nn),y=runif(nn))
    for ( exponent in c(1,2,3,5) ) {
        model <- lm(y~poly(x,exponent),obs)
        plot(obs$x,obs$y,pch=19,xlab="",ylab="",main=paste("Exponent:",exponent),xlim=c(0,1),ylim=c(0,1))
        lines(xx,predict(model,newdata=data.frame(x=xx)),col="red",lwd=2)
    }
dev.off()

png("sims_2.png",width=800,height=800)
    opar <- par(mfrow=c(4,4),mai=c(.3,.3,.3,.1),las=1)
    for ( jj in 1:4 ) {
        set.seed(jj)
        obs <- data.frame(x=runif(nn),y=runif(nn))
        for ( exponent in c(1,2,3,5) ) {
            model <- lm(y~poly(x,exponent),obs)
            plot(obs$x,obs$y,pch=19,xlab="",ylab="",main=paste("Exponent:",exponent),xlim=c(0,1),ylim=c(0,1))
            lines(xx,predict(model,newdata=data.frame(x=xx)),col="red",lwd=2)
        }
    }
    par(opar)
dev.off()

Answer 3

If your goal is interpolation, you typically want the simplest function that describes your observations and avoid overfitting.

Given that it is unusual to see physical laws and relationships which contain powers higher than 3, using higher order polynomials when there is no intuitive justification is taken to be a sign of overfitting.

See also this similar question.

Answer 4

Runge's phenomenon can lead to high-degree polynomials being much wigglier than the variation actually suggested by the data. An appeal of splines as a substitute for high-degree polynomials, particularly natural splines, is to allow nonmonotonicity and varying slopes without varying too wildly. I would be hard-pressed to find an example of real data for which a fourth- or higher-degree polynomial seems a more reasonable choice than a spline.

Answer 5

As an inveterate contrarian, I feel the need to amend the premise that high order polynomials shouldn't be used for interpolation. I would argue that the correct statement is "high order polynomials make poor interpolants unless properly regularized".

Indeed, it is quite popular (at least in academic circles) to conduct interpolation with polynomials of infinite degree, which is referred to as Kriging when a particular regularizer is used (the so called Reproducing Kernel Hilbert Space regularization, which intuitively speaking penalizes both a function's $L_2$ norm and its tendency to oscillate). And, of course, by infinite degree, we mean it is possible to construct a sequence of arbitrarily high degree polynomials which converge to this reasonable interpolant.

Additionally, recent work has shown that by increasing the degree to 2 or 3 times the size of the data (but keeping it finite), together with appropriate regularization, can still perform well. I can't find it just this moment, but this work was based on a similar discovery in the context of neural networks.

Generically, regularization prevents many of the issues brought up by other answerers.