# 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?

• This question is answered in several threads that can be found by searching our site for orthogonal polynomial – whuber May 1 '14 at 16:05
• @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 '14 at 16:17
• 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 '14 at 16:21
• 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 '14 at 17:21
• @Glen_b Yeah, well, I did type in ?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 used poly() for fitting, especially with ggplot – 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. – user13907 May 1 '14 at 20:14

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:

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


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.

• 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 '14 at 15:25
• 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 '14 at 15:59
• @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 '14 at 17:58
• 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 '14 at 21:18
• Oh, okay, that explanation in plain English now makes sense. Thank you. – user13907 May 2 '14 at 6:01

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.

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.seed(3)
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))

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

# predict y using modI
coef(modI)[1] + coef(modI)[2] * 3^1 + coef(modI)[3] * 3^2

# predict y using modp
# calculate the new x values using predict.poly()
x_poly <- stats:::predict.poly(object = poly(x,2), newdata = 3)
coef(modp)[1] + coef(modp)[2] * x_poly[1] + coef(modp)[3] * x_poly[2]


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.

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