On Mon, Apr 1, 2013 at 12:42 PM, Gabor Grothendieck <ggrothendieck@gmail.com
> wrote:

> On Mon, Apr 1, 2013 at 1:20 PM, Paul Johnson <pauljohn32@gmail.com> wrote:
> > Here's my little discussion example for a quadratic regression:
> >
> > http://pj.freefaculty.org/R/WorkingExamples/regression-quadratic-1.R
> >
> > Students press me to know the benefits of poly() over the more obvious
> > regression formulas.
> >
> > I think I understand the theory on why poly() should be more numerically
> > stable, but I'm having trouble writing down an example that proves the
> > benefit of this.
> >
> > I am beginning to suspect that because R's lm uses a QR decomposition
> under
> > the hood, the instability that we used to see when we inverted (X'X) is
> no
> > longer so serious as it used to be. When the columns in X are
> >
> > x1 x2 x3 x3squared
> >
> > then the QR decomposition is doing about the same thing as we would do if
> > we replaced x3 and x3squared by the poly(x3, 2) results.
> >
> > Or not. I'm not arguing that, I'm thinking out loud, hoping you'll
> correct
> > me.
> >
>
> # 1. One benefit is if you fit a higher degree polynomial using
> poly(x, n) the lower degree coefficients will agree with the fit using
> a lower n.
> # For example, compare these two:
> set.seed(123)
> y <- rnorm(10); x <- 1:10
> lm(y ~ poly(x, 2))
> lm(y ~ poly(x, 3))
>
> # however, that is not true if you don't use orthogonal polynomials.
> Compare these two:
> lm(y ~ poly(x, 2, raw = TRUE))
> lm(y ~ poly(x, 3, raw = TRUE))
>
>
Thanks, Gabor:

As usual, you are very helpful. I can't thank you enough.

I'm pasting your comments to the end of that R example file I mentioned
before, changing the examples to fit with the variable names I use there.

I think your reason #2 is the prize winner, I can't see anything wrong with
that one.

Reason #1 is not quite as strong, only holds strictly if the poly() is the
only predictor.  Users usually have other variables in the model.

If you--or anybody else--has other ideas about this, I am delighted to hear
your ideas.  I still hold out hope that the "numerical stability" argument
will find some traction.

pj



> # 2. A second advantage is that the t statistics are invariant under
> shift if you use orthogonal polynomials
> # compare these two:
> summary(lm(y ~ poly(x, 2)))
> summary(lm(y ~ poly(x+1, 2)))
>
> # however, that is not true if you don;t use orthogonal polynomials,
> Compare these two:
> summary(lm(y ~ poly(x, 2, raw = TRUE)))
> summary(lm(y ~ poly(x+1, 2, raw = TRUE)))
>
> --
> Statistics & Software Consulting
> GKX Group, GKX Associates Inc.
> tel: 1-877-GKX-GROUP
> email: ggrothendieck at gmail.com
>



-- 
Paul E. Johnson
Professor, Political Science      Assoc. Director
1541 Lilac Lane, Room 504      Center for Research Methods
University of Kansas                 University of Kansas
http://pj.freefaculty.org               http://quant.ku.edu

	[[alternative HTML version deleted]]

