# Can you add polynomial terms to multiple linear regression?

I am a little confused about when you should or shouldn't add polynomial terms to a multiple linear regression model. I know polynomials are used to capture the curvature in the data, but it always seems to be in the form of:

$$y = x_1 + x_2 + x_1^2 + x_2^2 + x_1x_2 + c$$

What if you know that there is a linear relationship between $$y$$ and $$x_1$$, but a non-linear relationship between $$y$$ and $$x_2$$? Can you use a model in the form of:

$$y = x_1 + x_2 + x_2^2 + c$$

I guess my question is, is it valid to drop the $$x_1^2$$ term and the $$x_1x_2$$ term, or do you have to follow the generic form of a polynomial regression model?

• Just for completeness note that if you have $x^2$ in the model you must have $x$ too. Search this site for principle of marginality for more info. I know you did not suggest doing it but the info might be helpful. Commented Sep 12, 2019 at 15:51
• It’s really weird to see something like $y = ax + bx^2$ and call it a linear regression model, but it definitely is. Remember that it has to do with linearity in the coefficients. I suggest Mathematical Monk’s video on this: m.youtube.com/watch?v=rVviNyIR-fI. Fifteen minutes might seem like a while to spend on this, but you’ll never forget that linear regression can involve nonlinear behavior. Maybe as a question to the OP: is $y= \sqrt{ax + bx^2}$ a linear regression model?
– Dave
Commented Sep 13, 2019 at 9:05
• @mdewey, your statement is too general to be correct and/or make sense in all settings. I can easily imagine a well-defined model with $x^2$ but without $x$. E.g. cases where the data generating process is $y=\beta_0+\beta_1 x^2+\varepsilon$ or where it is $y=\beta_0+\beta_1 z+\varepsilon$ where $x=\sqrt{z}$ (and $z>0$), or where the model is simply the better approximation (as compared to the one including $x$) to whatever data generating process we are facing. Commented Sep 13, 2019 at 12:29

In addition to @mkt's excellent answer, I thought I would provide a specific example for you to see so that you can develop some intuition.

Generate Data for Example

For this example, I generated some data using R as follows:

set.seed(124)

n <- 200
x1 <- rnorm(n, mean=0, sd=0.2)
x2 <- rnorm(n, mean=0, sd=0.5)

eps <- rnorm(n, mean=0, sd=1)

y = 1 + 10*x1 + 0.4*x2 + 0.8*x2^2 + eps


As you can see from the above, the data come from the model $$y = \beta_0 + \beta_1*x_1 + \beta_2*x_2 + \beta_3*x_2^2 + \epsilon$$, where $$\epsilon$$ is a normally distributed random error term with mean $$0$$ and unknown variance $$\sigma^2$$. Furthermore, $$\beta_0 = 1$$, $$\beta_1 = 10$$, $$\beta_2 = 0.4$$ and $$\beta_3 = 0.8$$, while $$\sigma = 1$$.

Visualize the Generated Data via Coplots

Given the simulated data on the outcome variable y and the predictor variables x1 and x2, we can visualize these data using coplots:

library(lattice)

coplot(y ~ x1 | x2,
number = 4, rows = 1,
panel = panel.smooth)

coplot(y ~ x2 | x1,
number = 4, rows = 1,
panel = panel.smooth)


The resulting coplots are shown below.

The first coplot shows scatterplots of y versus x1 when x2 belongs to four different ranges of observed values (which are overlapping) and enhances each of these scatterplots with a smooth, possibly non-linear fit whose shape is estimated from the data.

The second coplot shows scatterplots of y versus x2 when x1 belongs to four different ranges of observed values (which are overlapping) and enhances each of these scatterplots with a smooth fit.

The first coplot suggests that it is reasonable to assume that x1 has a linear effect on y when controlling for x2 and that this effect does not depend on x2.

The second coplot suggests that it is reasonable to assume that x2 has a quadratic effect on y when controlling for x1 and that this effect does not depend on x1.

Fit a Correctly Specified Model

The coplots suggest fitting the following model to the data, which allows for a linear effect of x1 and a quadratic effect of x2:

m <- lm(y ~ x1 + x2 + I(x2^2))


Construct Component Plus Residual Plots for the Correctly Specified Model

Once the correctly specified model is fitted to the data, we can examine component plus residual plots for each predictor included in the model:

library(car)

crPlots(m)


These component plus residual plots are shown below and suggest that the model was correctly specified since they display no evidence of nonlinearity, etc. Indeed, in each of these plots, there is no obvious discrepancy between the dotted blue line suggestive of a linear effect of the corresponding predictor, and the solid magenta line suggestive of a non-linear effect of that predictor in the model.

Fit an Incorrectly Specified Model

Let's play the devil's advocate and say that our lm() model was in fact incorrectly specified (i.e., misspecified), in the sense that it omitted the quadratic term I(x2^2):

m.mis <-  lm(y ~ x1 + x2)


Construct Component Plus Residual Plots for the Incorrectly Specified Model

If we were to construct component plus residual plots for the misspecified model, we would immediately see a suggestion of non-linearity of the effect of x2 in the misspecified model:

crPlots(m.mis)


In other words, as seen below, the misspecified model failed to capture the quadratic effect of x2 and this effect shows up in the component plus residual plot corresponding to the predictor x2 in the misspecified model.

The misspecification of the effect of x2 in the model m.mis would also be apparent when examining plots of the residuals associated with this model against each of the predictors x1 and x2:

par(mfrow=c(1,2))
plot(residuals(m.mis) ~ x1, pch=20, col="darkred")
abline(h=0, lty=2, col="blue", lwd=2)
plot(residuals(m.mis) ~ x2, pch=20, col="darkred")
abline(h=0, lty=2, col="blue", lwd=2)


As seen below, the plot of residuals associated with m.mis versus x2 exhibits a clear quadratic pattern, suggesting that the model m.mis failed to capture this systematic pattern.

Augment the Incorrectly Specified Model

To correctly specify the model m.mis, we would need to augment it so that it also includes the term I(x2^2):

m <- lm(y ~ x1 + x2 + I(x2^2))


Here are the plots of the residuals versus x1 and x2 for this correctly specified model:

par(mfrow=c(1,2))
plot(residuals(m) ~ x1, pch=20, col="darkred")
abline(h=0, lty=2, col="blue", lwd=2)
plot(residuals(m) ~ x2, pch=20, col="darkred")
abline(h=0, lty=2, col="blue", lwd=2)


Notice that the quadratic pattern previously seen in the plot of residuals versus x2 for the misspecified model m.mis has now disappeared from the plot of residuals versus x2 for the correctly specified model m.

Note that the vertical axis of all the plots of residuals versus x1 and x2 shown here should be labelled as "Residual". For some reason, R Studio cuts that label off.

Yes, what you're suggesting is fine. It's perfectly valid in a model to treat the response to one predictor as linear and a different one as being polynomial. It's also completely fine to assume no interactions between the predictors.

• Hi. Just an associated doubt. If both $x$ and $x^2$ are regressors, and say $x$ is strictly positive, would there be multicollinearity problem? Is it possible that the coefficients would have larger standard errors? Commented Sep 13, 2019 at 10:09
• @Dayne Good question! This is an issue that's discussed well here and here
– mkt
Commented Sep 13, 2019 at 10:56
• @mkt Those links give totally separate approaches than Ingolifs's idea about orthogonal polynomials. Any thoughts about the orthogonal polynomial approach?
– Dave
Commented Sep 13, 2019 at 11:03
• @Dave Don't know too much about it, I'm afraid. May be good to ask a new question about comparing the approaches.
– mkt
Commented Sep 13, 2019 at 12:24

You should take care to use Orthogonal polynomials if you're going to add polynomial terms.

Why? Without them you have a problem resembling colinearity. In certain regions, $$x^2$$ will look quite similar to $$x$$, and a parabola will do a decent job of fitting a straight line.

Observe:

These are polynomials of $$x,x^2,x^3$$.

Between 0 and 1.5 all three curves increase monotonically and while they curve differently from each other, they will give similar quality fits when x is positively correlated with y. By using all three in your code

y ~ x + x^2 + x^3


you are essentially using redundant shapes to fit your data with.

Orthogonal polynomials essentially give you added wiggle room when fitting, and each polynomial is essentially independent of the others.

Three polynomials of degree 1,2 and 3 generated by the poly() function in R.

Perhaps instead of explicitly thinking of them as polynomials, you instead think of them as 'trend components' or something:

$$x$$ represents 'more is always better' (or worse if the coefficient is negative). If you were doing a regression on music quality vs cowbell, you'd need this component.

$$x^2$$ represents a kind of goldilocks zone. If you were doing a regression on food tastiness vs amount of salt, this component would be salient.

$$x^3$$ is probably unlikely to be a dominant component on its own (the only example I could think of is How Much People Know vs How Much They Think They Know), but its presence will influence the shape and symmetry of the $$x$$ and $$x^2$$ terms.

There is a lot of hardout maths involved in orthogonal polynomials, but thankfully you only really need to know two things:

1. Orthogonal polynomials are only orthogonal over a certain region. The example I gave involves polynomials that are only orthogonal between 0 and 1.5.
2. If you're using R, use the poly() function to make your polynomials. poly(x,n) where n is the degree of the highest polynomial. It will make them orthogonal for you over the domain of your data $$x$$.
• This is extremely interesting and not something I’d heard before. Do you have a reference saying that this is necessary or useful, however? And useful for prediction or parameter inference? And do you know a Python command for “poly”?
– Dave
Commented Sep 13, 2019 at 10:35
• Various facets of orthogonal polynomial regression has been addressed here and here and lots more. Commented Sep 13, 2019 at 12:08
• Great answer, and thanks for bringing it up better than the post I had in mind. :) Commented Sep 13, 2019 at 12:18
• @Ingolifs: This is the best explanation of orthogonal polynomials I've read so far on this site! I agree with Dave that it would be helpful to comment in your answer on the usefulness of orthogonal polynomials for prediction or parameter inference. Commented Sep 13, 2019 at 15:16
• +1 but be careful with poly because if you try to predict with it it will refit the polynomials on the prediction sample; i.e. we get junk. Commented Sep 13, 2019 at 19:34

There's no rule that says you have to use all your variables. If you're trying to predict income, and your feature variables are SSN, years of schooling, and age, and you want to drop the SSN because you expect any correlation between it and income to be spurious, that's your judgment call to make. A model isn't invalid simply because there are other variables that you theoretically could have included, but didn't. Deciding what polynomial terms to include is just one of many decisions regarding feature selection.

While polynomial models often start with all terms being included, that's just so that all of them can be evaluated as to how much they are adding to the model. If it looks like a particular term is mostly just overfitting, it can be dropped in later iterations of the model. Regularization, such as lasso regression, can drop less useful variables automatically. Generally, it's better to start will a model that has too many variables, and whittle it down to the ones that are most useful, than to start with only the variables you think the model should rely on, and possibly miss out on a relationship you weren't expecting.