# Selecting best model based on linear, quadratic and cubic fit of data

I have a Java code that performs a linear regression on a set of data using the Gauss-Jordan elimination. It calculates a linear, quadratic and cubic functions using the least squares method.

My problem is choosing a function out of the three that best models my data. This is ofcourse, without plotting the curve. So lets say I have a set of data {x:1,2,3,4} {y:3,8,910}, and I get a linear, quadratic and cubic function for it using the least squares method.

How do I choose which function is the right model for my data?

• I have tried to uncover what I think is your real problem, but the last sentence is still cryptic for me... Please correct me if I'm wrong. – user88 Jan 28 '12 at 14:26
• What do you intend to use the model for? It's hard to know how to pick the 'best' model unless you specify an objective. – onestop Jan 28 '12 at 17:51
• To add to the comment from @onestop, if your data comes from a mechanistic or biological process, normally the function is based on that process. – Michelle Jan 28 '12 at 19:07
• The model is to be used to determine the o-notation complexity of a sorting algorithm. The x data will be the size, while the the y data will be the time it takes for it to be sorted. Error with the y data posted above, it should be {y:3,8,9,10}. – user1173951 Jan 28 '12 at 20:04
• Can you fit the cubic and then ask whether the coefficients of the cubic and quadratic terms are significantly different from zero? – Dzamo Norton Jan 29 '12 at 6:32

The general term for what you are asking about is model selection. You have a set of possible models, in this case something like \begin{aligned} y&=\beta_1x + \beta_0\\ y&=\beta_2x^2 + \beta_1x + \beta_0 \\ y&=\beta_3x^3 + \beta_2x^2 + \beta_1x + \beta_0 \\ \end{aligned} and you want to determine which of these models is most parsimonious with your data. We generally worry about parsimony rather than best-fitting (i.e, highest $R^2$) since a complex model could "over-fit" the data. For example imagine your timing data is generated by a quadratic algorithm, but there's a little bit of noise in the timing (random paging by the OS, clock inaccuracy, cosmic rays, whatever). The quadratic model might still fit reasonably well, but it won't be perfect. However, we can find a (very high order) polynomial that goes through each and every data point. This model fits perfectly but will be terrible at making future predictions and, obviously, doesn't match the underlying phenomenon either. We want to balance model complexity with the model's explanatory power. How does one do this?

There are many options. I recently stumbled upon this review by Zucchini, which might be a good overview. One approach is to calculate something like the AIC (Akaike Information Criterion), which adjusts each model's likelihood to take the number of parameters into account. These are often relatively easy to compute. For example, AIC is: $$AIC = 2k -2ln(L)$$ where L is the likelihood of the data given the model and k is the number of parameters (e.g., 2 for linear, 3 for quadratic, etc). You compute this criterion for each model, then choose the model with the smallest AIC.

Another approach is to use cross-validation (or something like that) to show that none of your models are over-fit. You could then select the best-fitting model.

That's sort of the general case. However, as @Michelle noted above, you probably don't want to be doing model selection at all if you know something about the underlying phenomemon. In this case, if you have the code or know the underlying algorithm, you should just trace through it to determine the algorithm's order.

Also, keep in mind that the Big-O order of the algorithm isn't technically defined in terms of the best-fit to the observed run time; it's more of a limiting property. You could feasibly have an algorithm with a massive linear component and a small quadratic component to its runtime, something like $$t(N) = 0.0000001n^2 + 999999999n$$ I would bet that a runtime-vs-input size plot for that would be pretty linear-looking over the ranges you're likely to test, but I believe the algorithm would technically be considered $O(n^2)$

• +1 for good general explanation, while still returning to the specific question. – Assad Ebrahim Oct 20 '12 at 14:00
• Yes - $t(N)$ will be considered quadratic complexity. The reason is because for large enough $n$ (say $n$ > 1 billion), the quadratic term will far outstrip the linear term, and this holds for any such equation, no matter how small the quadratic coefficient or how large the linear coefficient. – Assad Ebrahim Oct 20 '12 at 14:02
• @Matt Krause "However, we can find a (very high order) polynomial that goes through each and every data point. " any chance you have a reference for the proof for this? Ideally for the general case of p predictors. Thanks! – jmb Sep 11 '18 at 13:55
• @jmb, the magic phrase here is "interpolating polynomial". Here's a proof mathforcollege.com/nm/mws/gen/05inp/… and here's a math.se thread on the same topic: math.stackexchange.com/questions/837902/… – Matt Krause Sep 11 '18 at 14:55

Model selection will result in an estimate of the residual variance that is biased low. This will bias all other aspects of inference and inflate $R^2$. The unbiased estimate of $\sigma^2$ has in the denominator $n - p - 1$ where $p$ is the number of pre-specified parameters excluding the intercept. It is not clear that model selection will help. One can just use a regression spline with default knot locations (usually based on quantiles of the predictor distribution), choosing the number of knots based on what complexity the effective sample size will support. The R rms package makes this easy, and course notes at http://biostat.mc.vanderbilt.edu/rms will help. Removal of "non-significant" terms will hurt inference. The use of AIC can sometimes backfire if the modeling is not highly structured. I use AIC to tell me the optimum number $k$ of knots for a large number of predictors if I restrict each predictor to have $k$ knots. This structure minimizes the damage caused by model uncertainty.

I'm not 100% sure but if the model is time to sort data you should probably include a $x\log(x)$ or something like this (perhaps its $x^2\log(x)$) as I think this term appears in theoretical time complexity calculations for sorting data.