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?

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.
