Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

I'm trying to create a second order polynomial fit to some data I have. Let's say I plot this fit with ggplot():

ggplot(data, aes(foo, bar)) + geom_point() + 
       geom_smooth(method = "lm", formula = y ~ poly(x, 2))

I get:

So, a second order fit works quite well. I calculate it with R:

summary(lm(data$bar ~ poly(data$foo, 2)))

And I get:

lm(formula = data$bar ~ poly(data$foo, 2))

...

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)         3.268162   0.008282 394.623   <2e-16 ***
poly(data$foo, 2)1 -0.122391   0.096225  -1.272    0.206
poly(data$foo, 2)2  1.575391   0.096225  16.372   <2e-16 ***

....

Now, I would assume the formula for my fit is:

$$ \text{bar} = 3.268 - 0.122 \cdot \text{foo} + 1.575 \cdot \text{foo}^2 $$

But that just gives me the wrong values. For example, with $\text{foo}$ being 3 I would expect $\text{bar}$ to become something around 3.15. However, inserting into above formula I get:

$$ \text{bar} = 3.268 - 0.122 \cdot 3 + 1.575 \cdot 3^2 = 17.077 $$

What gives? Am I incorrectly interpreting the coefficients of the model?

share|improve this question
1  
This question is answered in several threads that can be found by searching our site for orthogonal polynomial –  whuber May 1 at 16:05
1  
@whuber If I had known that the problem was with "orthogonal polynomials", I probably would have found an answer. But if you don't know what to search for, it's a little hard. –  user13907 May 1 at 16:17
1  
You could also find answers by searching on poly, which appears prominently in your code. I put such information in comments for two reasons: (1) the links may help future readers as well as yourself and (2) they may help show you how to exploit our (somewhat idiosyncratic) search system. –  whuber May 1 at 16:21
    
I see, thank you. –  user13907 May 1 at 16:53
3  
You posted a question relating to your use of poly without typing ?poly in R first? That says 'Compute Orthogonal Polynomials' at the top in large friendly letters. –  Glen_b May 1 at 17:21

1 Answer 1

up vote 18 down vote accepted
+50

My detailed answer is below, but the general (i.e. real) answer to this kind of question is: 1) experiment, screw around, look at the data, you can't break the computer no matter what you do, so . . . experiment; or 2) RTFM.

Here is some R code which replicates the problem identified in this question, more or less:

# This program written in response to a Cross Validated question
# http://stats.stackexchange.com/questions/95939/
# 
# It is an exploration of why the result from lm(y_x+I(x^2))
# looks so different from the result from lm(y~poly(x,2))

library(ggplot2)


epsilon <- 0.25*rnorm(100)
x       <- seq(from=1, to=5, length.out=100)
y       <- 4 - 0.6*x + 0.1*x^2 + epsilon

# Minimum is at x=3, the expected y value there is
4 - 0.6*3 + 0.1*3^2

ggplot(data=NULL,aes(x, y)) + geom_point() + 
       geom_smooth(method = "lm", formula = y ~ poly(x, 2))

summary(lm(y~x+I(x^2)))       # Looks right
summary(lm(y ~ poly(x, 2)))   # Looks like garbage

# What happened?
# What do x and x^2 look like:
head(cbind(x,x^2))

#What does poly(x,2) look like:
head(poly(x,2))

The first lm returns the expected answer:

Call:
lm(formula = y ~ x + I(x^2))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.53815 -0.13465 -0.01262  0.15369  0.61645 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.92734    0.15376  25.542  < 2e-16 ***
x           -0.53929    0.11221  -4.806 5.62e-06 ***
I(x^2)       0.09029    0.01843   4.900 3.84e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2241 on 97 degrees of freedom
Multiple R-squared:  0.1985,    Adjusted R-squared:  0.182 
F-statistic: 12.01 on 2 and 97 DF,  p-value: 2.181e-05

The second lm returns something odd:

Call:
lm(formula = y ~ poly(x, 2))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.53815 -0.13465 -0.01262  0.15369  0.61645 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.24489    0.02241 144.765  < 2e-16 ***
poly(x, 2)1  0.02853    0.22415   0.127    0.899    
poly(x, 2)2  1.09835    0.22415   4.900 3.84e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2241 on 97 degrees of freedom
Multiple R-squared:  0.1985,    Adjusted R-squared:  0.182 
F-statistic: 12.01 on 2 and 97 DF,  p-value: 2.181e-05

Since lm is the same in the two calls, it has to be the arguments of lm which are different. So, let's look at the arguments. Obviously, y is the same. It's the other parts. Let's look at the first few observations on the right-hand-side variables in the first call of lm. The return of head(cbind(x,x^2)) looks like:

            x         
[1,] 1.000000 1.000000
[2,] 1.040404 1.082441
[3,] 1.080808 1.168146
[4,] 1.121212 1.257117
[5,] 1.161616 1.349352
[6,] 1.202020 1.444853

This is as expected. First column is x and second column is x^2. How about the second call of lm, the one with poly? The return of head(poly(x,2)) looks like:

              1         2
[1,] -0.1714816 0.2169976
[2,] -0.1680173 0.2038462
[3,] -0.1645531 0.1909632
[4,] -0.1610888 0.1783486
[5,] -0.1576245 0.1660025
[6,] -0.1541602 0.1539247

OK, that's really different. First column is not x, and second column is not x^2. So, whatever poly(x,2) does, it does not return x and x^2. If we want to know what poly does, we might start by reading its help file. So we say help(poly). The description says:

Returns or evaluates orthogonal polynomials of degree 1 to degree over the specified set of points x. These are all orthogonal to the constant polynomial of degree 0. Alternatively, evaluate raw polynomials.

Now, either you know what "orthogonal polynomials" are or you don't. If you don't, then use Wikipedia or Bing (not Google, of course, because Google is evil---not as bad as Apple, naturally, but still bad). Or, you might decide you don't care what orthogonal polynomials are. You might notice the phrase "raw polynomials" and you might notice a little further down in the help file that poly has an option raw which is, by default, equal to FALSE. Those two considerations might inspire you to try out head(poly(x, 2, raw=TRUE)) which returns:

            1        2
[1,] 1.000000 1.000000
[2,] 1.040404 1.082441
[3,] 1.080808 1.168146
[4,] 1.121212 1.257117
[5,] 1.161616 1.349352
[6,] 1.202020 1.444853

Excited by this discovery (it looks right, now, yes?), you might go on to try summary(lm(y ~ poly(x, 2, raw=TRUE))) This returns:

Call:
lm(formula = y ~ poly(x, 2, raw = TRUE))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.53815 -0.13465 -0.01262  0.15369  0.61645 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)              3.92734    0.15376  25.542  < 2e-16 ***
poly(x, 2, raw = TRUE)1 -0.53929    0.11221  -4.806 5.62e-06 ***
poly(x, 2, raw = TRUE)2  0.09029    0.01843   4.900 3.84e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2241 on 97 degrees of freedom
Multiple R-squared:  0.1985,    Adjusted R-squared:  0.182 
F-statistic: 12.01 on 2 and 97 DF,  p-value: 2.181e-05

There are at least two levels to the above answer. First, I answered your question. Second, and much more importantly, I illustrated how you are supposed to go about answering questions like this yourself. Every single person who "knows how to program" has gone through a sequence like the one above sixty million times. Even people as depressingly bad at programming as I am go through this sequence all the time. It's normal for code not to work. It's normal to misunderstand what functions do. The way to deal with it is to screw around, experiment, look at the data, and RTFM. Get yourself out of "mindlessly following a recipe" mode and into "detective" mode.

share|improve this answer
1  
I think this deserves a +6. I'll try to remember in a couple days when that becomes possible. FTR, I think it needn't be quite so sarcastic, but it does a good job of showing what orthogonal polynomials are / how they work, & showing the process you use to figure such things out. –  gung May 1 at 15:25
2  
Great answer, thank you. Although I am a little offended by a "RTFM" (but maybe that's just me): The problem is that in all I've read, at least with regard to doing linear regression in R, people sometimes do this, others do that. Frankly, I do not understand the Wikipedia entry on orthogonal polynomials. It doesn't occur to me why one would use this for regression if the coefficients you get are "wrong". I am not a mathematician — I try to follow the recipes because I'm not a learned cook, but I need to eat something nonetheless. –  user13907 May 1 at 15:59
5  
@user13907, that's not just you. This is indeed a good answer that deserves to be up-voted, but it would benefit from having a nicer tone. –  Waldir Leoncio May 1 at 17:58
4  
You don't really need to understand what orthogonal polynomials are here---you just need to understand that they are not what you want. Why might someone want orthogonal polynomials? Submit cov(poly(x,2)) to find that the covariance between the two terms in the polynomial is zero (up to roundoff error). This is the key property of orthogonal polynomials---their terms have zero covariance with each other. Sometimes it is convenient for your RHS variables to have zero correlation with each other. Their coefficients are not wrong, really, they just have to be interpreted differently. –  Bill May 1 at 21:18
    
Oh, okay, that explanation in plain English now makes sense. Thank you. –  user13907 May 2 at 6:01

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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