Fitting Polynomial Regression in R

A linear relationship between two variables x and y is one of the most common, effective and easy assumptions to make when trying to figure out their relationship. Sometimes however, the true underlying relationship is more complex than that, and this is when polynomial regression comes in to help.

Let see an example from economics: Suppose you would like to buy a certain quantity q of a certain product. If the unit price is p, then you would pay a total amount y. This is a typical example of a linear relationship. Total price and quantity are directly proportional. To plot it we would write something like this:

p <- 0.5
q <- seq(0,100,1)
y <- p*q
plot(q,y,type='l',col='red',main='Linear relationship')

The plot will look like this:

Now, this is a good approximation of the true relationship between y and q, however when buying and selling we might want to consider some other relevant information, like: Buying significant quantities it is likely that we can ask and get a discount, or buying more and more of a certain good we might be pushing the price up.
This can lead to a scenario like this one where the total cost is no longer a linear function of the quantity:

y <- 450 + p*(q-10)^3
plot(q,y,type='l',col='navy',main='Nonlinear relationship',lwd=3)

With polynomial regression we can fit models of order n > 1 to the data and try to model nonlinear relationships.

How to fit a polynomial regression

First, always remember use to set.seed(n) when generating pseudo random numbers. By doing this, the random number generator generates always the same numbers.

set.seed(20)

Predictor (q). Use seq for generating equally spaced sequences fast

q <- seq(from=0, to=20, by=0.1)

Value to predict (y):

y <- 500 + 0.4 * (q-10)^3

Some noise is generated and added to the real signal (y):

noise <- rnorm(length(q), mean=10, sd=80)
noisy.y <- y + noise

Plot of the noisy signal:

plot(q,noisy.y,col='deepskyblue4',xlab='q',main='Observed data')
lines(q,y,col='firebrick1',lwd=3)

This is the plot of our simulated observed data. The simulated datapoints are the blue dots while the red line is the signal (signal is a technical term that is often used to indicate the general trend we are interested in detecting).

Our model should be something like this: y = a*q + b*q² + c*q³ + cost

Let’s fit it using R. When fitting polynomials you can either use

model <- lm(noisy.y ~ poly(q,3))

Or

model <- lm(noisy.y ~ x + I(X^2) + I(X^3))

However, note that q, I(q^2) and I(q^3) will be correlated and correlated variables can cause problems. The use of poly() lets you avoid this by producing orthogonal polynomials, therefore I’m going to use the first option.

summary(model)

Call:
lm(formula = noisy.y ~ poly(q, 3))

Residuals:
     Min       1Q   Median       3Q      Max 
-212.326  -51.186    4.276   61.485  165.960 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  513.615      5.602   91.69   <2e-16 ***
poly(q, 3)1 2075.899     79.422   26.14   <2e-16 ***
poly(q, 3)2 -108.004     79.422   -1.36    0.175    
poly(q, 3)3  864.025     79.422   10.88   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 79.42 on 197 degrees of freedom
Multiple R-squared:  0.8031,	Adjusted R-squared:  0.8001 
F-statistic: 267.8 on 3 and 197 DF,  p-value:  0 

By using the confint() function we can obtain the confidence intervals of the parameters of our model.
Confidence intervals for model parameters:

confint(model, level=0.95)

                2.5 %     97.5 %
(Intercept)  502.5676  524.66261
poly(q, 3)1 1919.2739 2232.52494
poly(q, 3)2 -264.6292   48.62188
poly(q, 3)3  707.3999 1020.65097

Plot of fitted vs residuals. No clear pattern should show in the residual plot if the model is a good fit

plot(fitted(model),residuals(model))

Overall the model seems a good fit as the R squared of 0.8 indicates. The coefficients of the first and third order terms are statistically significant as we expected. Now we can use the predict() function to get the fitted values and the confidence intervals in order to plot everything against our data.

Predicted values and confidence intervals:

predicted.intervals <- predict(model,data.frame(x=q),interval='confidence',
                               level=0.99)

Add lines to the existing plot:

lines(q,predicted.intervals[,1],col='green',lwd=3)
lines(q,predicted.intervals[,2],col='black',lwd=1)
lines(q,predicted.intervals[,3],col='black',lwd=1)

Add a legend:

legend("bottomright",c("Observ.","Signal","Predicted"), 
       col=c("deepskyblue4","red","green"), lwd=3)

Here is the plot:

We can see that our model did a decent job at fitting the data and therefore we can be satisfied with it.

A word of caution: Polynomials are powerful tools but might backfire: in this case we knew that the original signal was generated using a third degree polynomial, however when analyzing real data, we usually know little about it and therefore we need to be cautious because the use of high order polynomials (n > 4) may lead to over-fitting. Over-fitting happens when your model is picking up the noise instead of the signal: even though your model is getting better and better at fitting the existing data, this can be bad when you are trying to predict new data and lead to misleading results.

A gist with the full code for this example can be found here.

Thank you for reading this post, leave a comment below if you have any question.