[R] Polynomial fitting

[Barry Rowlingson]

Dimitris Rizopoulos wrote:
> try this:
> 
> y <- c(15.51, 12.44, 31.5, 21.5, 17.89, 27.09, 15.02, 13.43, 18.18, 
> 11.32)
> x <- seq(3.75, 6, 0.25)
> coef(lm(y ~ x + I(x^2) + I(x^3)))

Or use the 'poly' function:

  > coef(lm(x~poly(y,3)))
  (Intercept) poly(y, 3)1 poly(y, 3)2 poly(y, 3)3
   4.8750000  -0.6233293   0.0312415  -0.6464533

  But hmmm those aren't the same coefficients?

  Because you need to use 'raw=TRUE' if you really want to fit the raw 
powers rather than orthogonal polynomials (which are better behaved than 
fitting the raw powers):

  > coef(lm(y~poly(x,3,raw=TRUE)))
             (Intercept) poly(x, 3, raw = TRUE)1 poly(x, 3, raw = TRUE)2
             -774.258364              472.172611              -91.633939
poly(x, 3, raw = TRUE)3
                5.800653

And there's your coefficients. See help(poly) for more.

Barry