1

I have fitted a simple 2nd order polynomial to time series data in the following form:

polyfit <- lm(y ~ poly(x,2))

I wish to extract the respective coefficients (A, B and C) from the fitted polynomial of the form y = Ax^2 + Bx + C. I naturally thought the answer would be found in polyfit$coefficients within the polyfit object but these coefficients are not correct. I have tried some very simple data sets and compared with excel and whilst the poly curve fits are identical in R and excel, the A,B and C coefficints obtained from excel are correct but those obtained from the polyfit object arent? Have I extracted the incorrect information from the polyfit object? It would be more convenient for me to extract the coefficients directly from R for my purposes? Can anyone help?

RoyalTS
  • 9,545
  • 12
  • 60
  • 101
Phil W
  • 115
  • 1
  • 7
  • Provide a reproducible example with data (fake or real) and you will be much more likely to get the answer you're looking for. – ndoogan May 10 '13 at 02:46
  • 1
    If you want a regular polynomial you have to do `poly(x,2,raw=TRUE)` for otherwise R will fit an orthogonal polynomial – RoyalTS May 10 '13 at 02:46
  • Thank-you RoyalTS all works perfectly now, cheers Phil – Phil W May 10 '13 at 04:15
  • @RoyalTS `raw = TRUE` also changes the fitting procedure and thus might lead to a different polynomial. People argue that raw = FALSE gives the better fit. So is there an easy way to obtain formulas for the orthogonal polynomials? – jarauh Dec 08 '17 at 09:53

1 Answers1

5

By default poly fits orthogonal polynomials. Essentially it uses the following ideas...

x <- 1:10

# in this case same as x-mean(x)
y1 <- residuals(lm(x ~ 1)) 
# normalize to have unit norm
y1 <- y1/sqrt(sum(y1^2))

y2 <- residuals(lm(y1^2 ~ y1))
y2 <- y2/sqrt(sum(y2^2))

y3 <- residuals(lm(x^3 ~ y2 + y1))
y3 <- y3/sqrt(sum(y3^2))

cbind(y1, y2, y3)    
poly(x, 3)

to construct a set of orthonormal vectors that still produce the same predictions. If you just want to get a polynomial in the regular way then you want to specify raw=TRUE as a parameter.

y <- rnorm(20)
x <- 1:20
o <- lm(y ~ poly(x, 2, raw = TRUE))
# alternatively do it 'by hand'
o.byhand <- lm(y ~ x + I(x^2))
Dason
  • 60,663
  • 9
  • 131
  • 148
  • Thank-you Dason I was missing the "raw=TRUE" parameter, all works perfectly now, much obliged Phil – Phil W May 10 '13 at 03:57
  • @PhilW While `raw = TRUE` gives the desired output, it also changes the fitting procedure and thus might lead to a different polynomial. People argue that `raw = FALSE` gives the better fit. So is there an easy way to obtain formulas for the orthogonal polynomials? – jarauh Dec 08 '17 at 09:49