Getting a second-order polynomial trend line from a set of data

Question

Alright, so I have about a thousand datapoints that I'm plotting on a chart (scatter plot).

Here's a few of the records:

x              y
2.426032708   10
9.449923509   24
4.409997771   24
-3.125392294  11
13.04820475   26

When looking at the full data on a chart, I can visually see a 2nd order polynomial trendline is the way to go. With my full set of data, how can I calculate the formula of the line?

Everything I've found in google has told me how to do it in Excel, and that's not what I want. I want to get the formula so I can predict results at other points.

(Please be gentle. I'm a programmer, not a stats guy. If you need anything else, don't hesitate to ask)

Answer 1

As a programmer, you will find an algorithm to be even better than a formula, especially if the algorithm is based on simple steps. There is one, "sequential matching," that requires no matrix algebra and no specialized mathematical procedures (like Cholesky decomposition, QR decomposition, or matrix pseudo-inversion). It is easy to code, easy to test, extensible, and reasonably efficient.

---

A good algorithm for this situation is derived from the idea of "matching" developed by Tukey and Mosteller, as described in my answer at "How to Normalize Regression Coefficients". For quadratic regression it involves these simple ideas:

1. There are three independent variables: a constant (usually set to $1$ for simplicity and easy interpretation), $x$ itself, and $x^2$. For this purpose, you compute $x^2$ once and for all and then treat it henceforth as if it had no mathematical relationship to $x$ whatsoever.

2. After you have fit (or "matched") a single dependent variable $y$ to a single independent variable $x$, thereby producing a formula in the form $y = \beta_{y\cdot x} x + \text{ error },$ you take the fit away from the dependent variable by subtracting the fit, writing $y - \beta_{y\cdot x} x$ for what is left over. This (the "error" term above) is the residual. As a matter of notation, let the residual after matching $y$ to $x$ be called $y_{\cdot x}.$

3. To fit a dependent variable $y$ to $n\ge 2$ independent variables $x_1, x_2, \ldots, x_n$, you separately match $y$ and the last $n-2$ variables to the first one, producing dependent residual $y_{\cdot x_1}$ and independent residuals $x_{2\cdot x_1}, x_{3\cdot x_1}, \ldots, x_{n\cdot x_1}.$ Now proceed recursively to fit the dependent residual to the $n-1$ independent residuals.

Provided, then, that you have a routine to match a dependent variable to a single independent variable, you can figure out the residual for any multivariate linear fitting problem with practically no more coding. If you know the residuals, you know the prediction, so at this point it's just a matter of finding the coefficients. The most direct way to proceed is to do the algebra to work out the proper combination of all the appropriate $\beta$'s. This is worked out for the case $n=2$ in the answer previously referenced. The R code below shows it for quadratic regression. Rather than coding it in a loop over the $x_i$, I have unrolled that loop to exhibit every one of the steps that is needed, showing the basic simplicity of the algorithm. I have also avoided using any R-specific idioms, apart from its readiness to combine two vectors (such as x and y) component-by-component when they are added, subtracted, multiplied, or divided.

The output of this sample program consists of the coefficients of $1$, $x$, and $x^2$, named beta.0, beta.1, and beta.2, respectively, followed by the same coefficients as computed with R's built-in regression function:

> c(beta.0, beta.1, beta.2)
[1] 18.5575094  2.1555036 -0.1092891

> coef(lm(y ~ x + I(x^2)))
(Intercept)           x      I(x^2) 
 18.5575094   2.1555036  -0.1092891 

The perfect agreement attests to the correctness and precision of this sequential matching process.

#
# Data.
#
n <- 32
set.seed(17)
x <- runif(n, -5, 20)
y <- floor(30 - ((x-10)/3)^2 + rnorm(n, sd=3))
#
# Linear regression ("matching")
#
fit.ls <- function(y, x) sum(x*y) / sum(x*x) # Returns the coefficient
#
# The additional variables are a constant `one` and the squared x's.
#
one <- rep(1, length(x))
x2 <- x*x
#
# Step 1: Match everything to `one`.
#
beta.x.1 <- fit.ls(x, one)
beta.x2.1 <- fit.ls(x2, one)
beta.y.1 <- fit.ls(y, one)
#
#         Compute the residuals.
x.1 <- x - one * beta.x.1
x2.1 <- x2 - one * beta.x2.1
y.1 <- y - one * beta.y.1
#
# Step 2: Match the residuals to `x`.
#
beta.x2.1.x <- fit.ls(x2.1, x.1)
beta.y.1.x <- fit.ls(y.1, x.1)
#
#         Compute the residuals.
#
x2.1.x <- x2.1 - x.1 * beta.x2.1.x
y.1.x <- y.1 - x.1 * beta.y.1.x
#
# Step 3: Match the residuals to x^2.
#
beta.y.1.x.x2 <- fit.ls(y.1.x, x2.1.x)
y.1.x.x2 <- y.1.x - x2.1.x * beta.y.1.x.x2
#
# Combine the coefficients into the full formula for `y` in terms of
# `one`, `x`, and x^2.
#
beta.2 <- beta.y.1.x.x2
beta.1 <- beta.y.1.x - beta.2*beta.x2.1.x
beta.0 <- beta.y.1 - beta.y.1.x*beta.x.1 - beta.2*(beta.x2.1 - beta.x2.1.x*beta.x.1)
c(beta.0, beta.1, beta.2)
#
# Compare this to the output of `R`'s built-in multiple regression.
#
coef(lm(y ~ x + I(x^2)))
#
# Plot the results.
#
par(mfrow=c(1,1))
plot(x, y)
curve(beta.0 + beta.1*x + beta.2*x^2, lwd=2, col="Red", add=TRUE)

Answer 2

Your model will be:

$$y_i = \beta_0 + \beta_1 x_i + \beta_2 x_i^2$$

Where $\beta_0$, $\beta_1$ and $\beta_2$ are parameters to be estimated from the data. Standard practice is to find values of these parameters such that the sum of squares:

$$ \sum_{i=1}^{n}\left[ y_i - (\beta_0 + \beta_1 x_i + \beta_2x_i^2) \right]^2$$

is minimized. In words, we are looking for coefficients of the polynomial such that the fitted values of the polynomial are as close to the observations as possible. In matrix/vector notation what we want is the vector $\vec{\beta}$ which satisfies:

$$\vec{y} = \matrix{X}\vec{\beta} $$

where $\vec{\beta} = [\beta_0,\beta_1,\beta_2]^T$, $\vec{y}=[y_1,\cdots,y_n]^T$ and

$$ \matrix{X} = \left[\array{1,x_1,x_1^2\\1,x_2,x_2^2\\\cdots\\1,x_n,x_n^2}\right] $$

As we cannot invert the matrix $\matrix{X}$ (it's not square for one thing), we solve the equation as follows:

$$ \matrix{X}^{T}\vec{y} = \matrix{X}^{T}\matrix{X}\vec{\beta} \\ (\matrix{X}^{T}\matrix{X})^{-1}\matrix{X}^{T}\vec{y} = \vec{\beta} $$

Actually, computationally speaking, inverting the matrix $\matrix{X}^{T}\matrix{X}$ is not the most efficient way of solving for $\vec{\beta}$ (see the Wikipedia entry on linear least squares) but it will give the right results.

Once you have found the vector of polynomial coefficients then you can evaluate the polynomial for any new value of $x$ you wish to predict.

An example of how to find values of $\vec{\beta}$ in R, manually:

X <- cbind(1,x,x^2)
beta <- solve(t(X) %*% X) %*% t(X) %*% y

Or just using the lm() function to fit a linear model which will exploit those efficient computational methods (and provide many more details on the model fit):

fit <- lm(y~x+I(x^2))
beta <- coef(fit)

A note for the uninitiated: when we specify the model in the call to lm() using R model syntax an intercept term (i.e. the column of $1$s in $\matrix{X}$) is automatically implicitly included. Hence we can just write y~x+I(x^2) rather than y~1+x+I(x^2).