How to interpret coefficients from a polynomial model fit? 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?
 A: There's an interesting approach to interpretation of polynomial regression by Stimson et al. (1978). It involves rewriting
$Y = \beta_{0} + \beta_{1} X + \beta_{2} X^{2} + u$
as
$Y = m + \beta_{2} \left( f - X \right)^{2} + u$
where $m = \beta_{0} - \left. \beta_{1}^{2} \right/ 4 \beta_{2}$ is the minimum or maximum (depending on the sign of $\beta_{2}$) and $f = \left. -\beta_{1} \right/ 2 \beta_{2}$ is the focal value. It basically transforms the three-dimensional combination of slopes into a parabola in two dimensions. Their paper gives an example from political science.
A: 'poly' performs Graham-Schmidt ortho-normalization on the polynomials 1, x, x^2, ..., x^deg
For example this function does the same thing as 'poly' without returning 'coef' attributes of course.
MyPoly <- 
function(x, deg)
{
    n <- length(x)
    ans <- NULL
    for(k in 1:deg)
    {
        v <- x^k
        cmps <- rep(0, n)
        if(k>0) for(j in 0:(k-1)) cmps <- cmps + c(v%*%ans[,j+1])*ans[,j+1]
        p <- v - cmps
        p <- p/sum(p^2)^0.5
        ans <- cbind(ans, p)
    }
    ans[,-1]
}

I landed on this thread because I was interested in the functional form. So how do we express the result of 'poly' as an expression?  Just invert the Graham-Schmidt procedure. You'll end up with a mess!
A: If you just want a nudge in the right direction without quite so much judgement: poly() creates orthogonal (not correlated) polynomials, as opposed to I(), which completely ignores correlation between the resultant polynomials. Correlation between predictor variables can be a problem in linear models (see here for more information on why correlation can be problematic), so it's probably better (in general) to use poly() instead of I(). 
Now, why do the results look so different? Well, both poly() and I() take x and convert it into a new x. In the case of I(), the new x is just x^1 or x^2. In the case of poly(), the new x's are much more complicated. If you want to know where they come from (and you probably don't), you can get started here or the aforementioned Wikipedia page or a textbook. 
The point is, when you're calculating (predicting) y based on a particular set of x values, you need to use the converted x values produced by either poly() or I() (depending which one was in your linear model). So:
library(ggplot2)    

# set the seed to make the results reproducible.
set.seed(3)

#### simulate some data ####
# epsilon = random error term
epsilon <- 0.25*rnorm(100)
# x values are just a sequence from 1 to 5
x       <- seq(from=1, to=5, length.out=100)
# y is a polynomial function of x (plus some error)
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

# visualize the data (with a polynomial best-fit line)
ggplot(data=NULL,aes(x, y)) + geom_point() + 
   geom_smooth(method = "lm", formula = y ~ poly(x, 2))


#### Model the data ####
# first we'll try to model the data with just I()
modI <- lm(y~x+I(x^2)) 
# the model summary looks right
summary(modI)

# next we'll try it with poly()
modp <- lm(y ~ poly(x, 2))
# the model summary looks weird
summary(modp)

#### make predictions at x = 3 based on each model ####
# predict y using modI
# all we need to do is take the model coefficients and plug them into the formula: intercept + beta1 * x^1 + beta2 * x^2
coef(modI)[1] + coef(modI)[2] * 3^1 + coef(modI)[3] * 3^2

(Intercept)
   3.122988
# predict y using modp
# this takes an extra step.
# first, calculate the new x values using predict.poly()
x_poly <- stats:::predict.poly(object = poly(x,2), newdata = 3)
# then use the same formula as above, but this time instead of the raw x value (3), use the polynomial x-value (x_poly)
coef(modp)[1] + coef(modp)[2] * x_poly[1] + coef(modp)[3] * x_poly[2]

(Intercept)
   3.122988 
In this case, both models return the same answer, which suggests that correlation among predictor variables is not influencing your results. If correlation were a problem, the two methods would predict different values. 
