My detailed answer is below, but the general (i.e. real) answer to this kind of question is: 1) experiment, mess around, look at the data, you can't break the computer no matter what you do, so ... experiment; or 2) read the documentation.
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.
poly
without typing?poly
in R first? That says 'Compute Orthogonal Polynomials' at the top in large friendly letters. $\endgroup$?poly
to understand the syntax. Admittedly, I have only little knowledge of the concepts behind it. I didn't know that there was something else (or such a big difference between "normal" polynomials and orthogonal polynomials), and the examples I saw online all usedpoly()
for fitting, especially withggplot
– so why wouldn't I just use that and be confused if the outcome was "wrong"? Mind you, I am not skilled in mathematics—I'm merely applying what I've seen others do, and trying to understand it. $\endgroup$