Raw or orthogonal polynomial regression?

I want to regress a variable $y$ onto $x,x^2,\ldots,x^5$. Should I do this using raw or orthogonal polynomials? I looked at question on the site that deal with these, but I don't really understand what's the difference between using them.

Why can't I just do a "normal" regression to get the coefficients $\beta_i$ of $y=\sum_{i=0}^5 \beta_i x^i$ (along with p-values and all the other nice stuff) and instead have to worry whether using raw or orthogonal polynomials? This choice seems to me to be outside the scope of what I want to do.

In the stat book I'm currently reading (ISLR by Tibshirani et al) these things weren't mentioned. Actually, they were downplayed in a way.
The reason is, AFAIK,that in the lm() function in R, using y ~ poly(x, 2) amounts to using orthogonal polynomials and using y ~ x + I(x^2) amounts to using raw ones. But on pp. 116 the authors say that we use the first option because the latter is "cumbersome" which leaves no indication that these commands actually to completely different things (and have different outputs as a consequence).
(third question) Why would the authors of ISLR confuse their readers like that?

• @Sycorax I know that poly has something to do with orthogonal polynomials and I(x^2) doesn't (though I don't know the details) - but still, why would the authors of ISLR then recommend a method that does not work ? It seems very misleading if both command seem to do the same, but only one actually is ok. – l7ll7 Jan 26 '17 at 16:15
• @gung I looked at the documentation of poly and spent already a while with this problem, but I can't figure out why poly(x,2) and x+I(x^2) make a difference? Could you please enlighten me here in the comments, if the question is offtopic? – l7ll7 Jan 26 '17 at 16:18
• @gung I did a completey re-edit of my question. This choice raw/orthogonal is confusing me even more - previously I thought this was just a minor R technicality, that I didn't understand, but now it seems to be a fullblown stat problem that hinders me of doing coding a regression that should not be that difficult to code. – l7ll7 Jan 26 '17 at 16:37
• Why the downvote? I just completely re-edited my question. – l7ll7 Jan 26 '17 at 16:38
• @gung That actually confused me more than it helped. Previously I thought that I should just go with orthogonal polynomials, because that seemed to be the right way, but in that answer raw polynomials are used. Suprisingly, everyone on the net is screaming "RTFM", but there is actually not clear answer, when to use what. (Your link also doesn't give an answer to this, just an example, when orth. pol. go wrong) – l7ll7 Jan 27 '17 at 10:31

I'm not 100% certain here, but I believe the answer is less about numeric stability (though that plays a role) and more about reducing correlation.

In essence -- the issue boils down to the fact that when we regress against a bunch of high order polynomials, the covariates we are regressing against become highly correlated. Example code below:

x = rnorm(1000)
raw.poly = poly(x,6,raw=T)
orthogonal.poly = poly(x,6)
cor(raw.poly)
cor(orthogonal.poly)


This is tremendously important. As the covariates become more correlated, our ability to determine which are important (and what the size of their effects are) erodes rapidly. At the limit, if we had two variables that were fully correlated, when we regress them against something, its impossible to distinguish between the two -- you can think of this as an extreme version of the problem, but this problem affects our estimates for lesser degrees of correlation as well. Thus in a real sense -- even if numerical instability wasn't a problem -- the correlation from higher order polynomials does tremendous damage to our estimation routines. This will manifest as larger standard errors that you would otherwise see (see example regression below). For this reason, we might choose to orthogonalize our polynomials before regressing them.

y = x*2 + 5*x**3 - 3*x**2 + rnorm(1000)
raw.mod = lm(y~poly(x,6,raw=T))
orthogonal.mod = lm(y~poly(x,6))
summary(raw.mod)
summary(orthogonal.mod)


If you run this code, interpretation is a touch hard because the coefficients all change and so things are hard to compare. Looking at the T-stats though, we can see that the ability to determine the coefficients was MUCH larger with orthogonal polynomials. For the 3 relevant coefficients, I got t-stats of (560,21,449) for the orthogonal model, and only (28,-38,121) for the raw polynomial model. This is a huge difference for a simple model with only a few relatively low order polynomial terms that mattered.

That is not to say that this comes without costs. There are two primary costs to bear in mind. 1) we lose some interpretability with orthogonal polynomials. We might understand what the coefficient on x**3 means, but interpreting the coefficient on x**3-3x (the third hermite poly -- not necessarily what you will use) can be much harder. Second -- when we say that these are polynomials are orthogonal -- we mean that they are orthogonal with respect to some measure of distance. Picking a measure of distance that is relevant to your situation can be difficult. However, having said that, I believe that the poly function is designed to choose such that it is orthogonal with respect to covariance -- which is useful for linear regressions.

Why can't I just do a "normal" regression to get the coefficients?

Because it is not numerically stable. Remember computer is using fixed number of bits to represent a float number. Check IEEE754 for details, you may surprised that even simple number $0.4$, computer need to store it as $0.4000000059604644775390625$. You can try other numbers here

Using raw polynomial will cause problem because we will have huge number. Here is a small proof: we are comparing matrix condition number with raw and orthogonal polynomial.

> kappa(model.matrix(mpg~poly(wt,10),mtcars))
[1] 5.575962
> kappa(model.matrix(mpg~poly(wt,10, raw = T),mtcars))
[1] 2.119183e+13


You can also check my answer here for an example.

Why are there large coefficents for higher-order polynomial

• You seem to be using single precision floats and quoting them to quadruple precision! How did that happen? Except for GPUs, almost all statistical computation uses at least double precision. E.g., in R the output of print(0.4, digits=20) is 0.40000000000000002. – whuber Mar 9 '18 at 21:43

In a quick simulation with a degree-5 polynomial, N=1000 points, random parameter values (scaled down by $$10^{-d}$$ to keep the response variable within 2 orders of magnitude), and uncorrelated noise $$\sim$$ half the overall variation in response variable, the $$R^2$$'s of the two models were identical. So were the $$adj ~R^2$$'s. Predictive power is the same. But the standard errors of parameter values for the orthogonal model were equal to or orders of magnitude below the raw model.