How to interpret coefficients from a polynomial model fit?

Question

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?

Answer 1

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.

Answer 2

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.

Answer 3

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.

Answer 4

'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!