11
$\begingroup$

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?

|cite|improve this question
$\endgroup$
  • 1
    $\begingroup$ You should give some thought as to which p-values are different when you fit the same model to the same data using raw & orthogonal polynomials, & their interpretation. What about the model predictions? $\endgroup$ – Scortchi - Reinstate Monica♦ Jan 27 '17 at 11:07
  • $\begingroup$ @Scortchi I added the relevant information to my question. $\endgroup$ – l7ll7 Jan 27 '17 at 11:21
  • 4
    $\begingroup$ Another good reason to use orthogonal polynomials is numerical stability; the associated design matrix for fitting in the monomial basis can be quite ill-conditioned for high-degree fitting since the higher-order monomials are "very nearly linearly dependent" (a concept that could be made more mathematically precise), while the design matrix for orthogonal polynomials are a bit better behaved. I discussed the equispaced abscissas (Gram) case here, but the deal is similar in the non-equispaced case. $\endgroup$ – J. M. is not a statistician Jan 27 '17 at 14:28
  • $\begingroup$ (Nevertheless, one should not fit to high-degree polynomials without a good reason for doing so.) $\endgroup$ – J. M. is not a statistician Jan 27 '17 at 14:29
7
$\begingroup$

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.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ +1 Finally a clear answer! Thank you! Before I accepted, can you please tell me, are there any other statistics, like R^2 or the F-statistic that I should better read of orthogonal plot summary than the raw one? Besides plotting the variables, is the fit using raw polynomials good for anything else in this scenario? $\endgroup$ – l7ll7 Jan 27 '17 at 13:20
  • $\begingroup$ And when I have multiple predictors, does the same hold true? $\endgroup$ – l7ll7 Jan 27 '17 at 13:20
  • $\begingroup$ How would you "use orthogonal polynomials to decide whether you want to include a quadratic term or not"? $\endgroup$ – Scortchi - Reinstate Monica♦ Jan 27 '17 at 13:26
  • 1
    $\begingroup$ Point is, the test of the highest-order effect, the quadratic in this case, is the same whether you use raw or orthogonal polynomials. So why bother with orthogonal polynomials? $\endgroup$ – Scortchi - Reinstate Monica♦ Jan 27 '17 at 13:34
  • 4
    $\begingroup$ Well, of course you simply shouldn't do those marginal tests in that model; you should re-fit after discarding the highest-order effect. Orthogonal polynomials spare you the bother, allowing an easy step-down procedure - perhaps you could illustrate with a cubic term. $\endgroup$ – Scortchi - Reinstate Monica♦ Jan 27 '17 at 13:44
3
$\begingroup$

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.

|cite|improve this answer
$\endgroup$
  • 1
    $\begingroup$ (+1) No difference in terms of prediction; & it might be said no difference in terms of any meaningful inference. $\endgroup$ – Scortchi - Reinstate Monica♦ Jan 27 '17 at 13:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.