If you can't do it orthogonally, do it raw (polynomial regression)

Question

When performing polynomial regression for $Y$ onto $X$, people sometimes use raw polynomials, sometimes orthogonal polynomials. But when they use what seems completely arbitrary.

Here and here raw polynomials are used. But here and here, orthogonal polynomials seems to give the correct results. What, how, why?!

In contrast to that, when learning about polynomial regression from a textbook (e.g. ISLR), that does not even mention raw or orthogonal polynomials - just the model to be fitted is given.

So when do we have to use what?
And why are the individual p-values for $X$, $X^2$ etc. differing a lot between these two values?

Answer 1

The variables $X$ and $X^2$ are not linearly independent. So even if there is no quadratic effect, adding $X^2$ to the model will modify the estimated effect of $X$.

Let’s have a look with a very simple simulation.

> x <- runif(1e3)
> y <- x + rnorm(length(x))
> summary(lm(y~x))

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.03486    0.06233  -0.559    0.576    
x            1.05843    0.10755   9.841   <2e-16 ***

Now with a quadratic term in the model to fit.

> summary(lm(y~x+I(x^2)))

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.03275    0.09528   0.344    0.731
x            0.65742    0.44068   1.492    0.136
I(x^2)       0.39914    0.42537   0.938    0.348

Of course the omnibus test is still significant, but I think the result we are looking is not this one. The solution is to use orthogonal polynomials.

 > summary(lm(y~poly(x,2)))

 Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.49744    0.03098  16.059   <2e-16 ***
poly(x, 2)1  9.63943    0.97954   9.841   <2e-16 ***
poly(x, 2)2  0.91916    0.97954   0.938    0.348    

Note that the coefficients of x in the first model and of poly(x,2)1 in the second model are not equal, and even the intercepts are different. This is because poly delivers orthonormal vectors, which are also orthogonal to the vector rep(1, length(x)). So poly(x,2)1 is not x but rather (x -mean(x))/sqrt(sum((x-mean(x))**2))...

An important point is that the Wald tests, in this last model, are independent. You can use orthogonal polynomials to decide up to which degree you want to go, just by looking at the Wald test: here you decide to keep $X$ but not $X^2$. Of course you would find the same model by comparing the first two fitted models, but it is simpler this way — if you consider going up to higher degrees, it is really much simpler.

Once you have decided which terms to keep, you may want to go back to raw polynomials $X$ and $X^2$ for interpretability or for prediction.

Answer 2

To give a naive assessment of the situation:

generally: suppose you have two different system of basis functions $\{p_n\}_{n=1}^\infty$, as well as $\{\tilde{p}\}_{n=1}^\infty$ for some function (hilbert-) space, usual $L_2([a,b])$, i.e. the space of all square-integrable functions.

This means that each of the two bases can be used to explain each element of $L_2([a,b])$, i.e. for $y \in L_2([a,b])$ you have for some coefficients $\theta_n$ and $\tilde{\theta}_n \in \mathbb{R}$, $n=1,2,\dots$ (in the $L_2$-sense): $$ \sum_{n=1}^\infty \tilde{\theta}_n \tilde{p}_n = y= \sum_{n=1}^\infty \theta_n p_n.$$

However, on the other hand, if you truncate both sets of basis functions at some number $k<\infty$, i.e. you take $$\{p_n\}_{n=1}^k$$ as well as $$\{\tilde{p}\}_{n=1}^k,$$ these truncated sets of basis functions are very likely two describe "different parts" of $L_2([a,b])$.

However, here in the special case where one basis, $\{\tilde{p}\}_{n=1}^\infty$, is just an orthogonalization of the other basis, $\{p_n\}_{n=1}^\infty$, the overall prediction of $y$ will be the same for each truncated model ($\{p\}_{n=1}^k$ and their orthogonalized counterpart will describe the same $k$-dimensional subspace of $L_2([a,b])$).

But each individual basis function from the two "different" bases will yield a different contribution to this predcition (obviously as the functions/predictors are different!) resulting in different $p$-values and coefficients.

Hence, in terms of prediction there is (in this case) no difference.

From a computational point of view a model matrix consisting of orthogonal basis functions have nice numerical/computational properties for the least squares estimator. While at the same time from the statistical point of view, the orthogonalization results in uncorrelated estimates, since $var(\hat{\tilde{\theta}}) = I \sigma²$ under the standard assumptions.

---

The natural question arises if there is a best truncated basis system. However the answer to the question is neither simple nor unique and depends for example on the definition of the word "best", i.e. what you are trying to archive.