Goodness of fit. How to evaluate if polynomial of order n+1 gives statistically better fit than polynomial of order n?

I fit polynomials with increasing order to some data. What is the best way to evaluate if the additional parameter of polynomial of order n+1 provides a statistically significant better fit than the previous polynomial (order n) with a given confidence level (e.g. 99%)?

So far, I used the F-test in the following way: $$F = \frac{[rss_n - rss_{n+1}] / (p_{n+1} - p_n)}{rss_{n+1} / (N-p_{n+1})}$$

where n is the order of the polynomial, rss$_n$ is its related Residual Sum of Squares, N is the amount of measurements and p$_n$ is the number of parameters of polynomial of order n. This formula comes from: https://en.wikipedia.org/wiki/F-test#Regression_problems. To give context, in my case p$_{n+1}$-p$_n$ is always equal to 1 and N is typically in the order of a few thousands.

Then I checked if F>F_critical(df$_1$,df$_2$,99%), where df$_1$=p$_{n+1}$ - p$_n$=1 and df$_2$=N-p$_{n+1}$ (df=degrees of freedom). If not, the polynomial of order n+1 does not provide a better fit with 99% confidence and I select the previous polynomial as the best fit.

Is this the way to deal with this problem? Are there other ways to do this (maybe using $\chi^2$)?

I have doubts because in some cases with my data it happens that each polynomial until n=9 keeps providing a better fit with 99% confidence and this seems a bit suspicious (particularly because the fitted polynomial seem fairly similar visually). Here is a figure showing polynomials of order 9 (dark green) and 6 (light green). According to the formula above, each polynomial provides a statistically better fit than the previous with 99% confidence interval. However, I think there's a great deal of overfitting instead. Data for this test (black dots in the figure) can be downloaded here.

• Your F statistic is a bit strange: (a) the numerator is negative, (b) a nth order polynom has $n+1$ parameters and (c) the denominator usually refers to the bigger model, not to the smaller. – Michael M Jul 4 '15 at 10:11
• Hi Michael, thanks for your comments. b) That's right, I rewrote the formula so that this is clear. Note that in my case the difference in models' parameters is always 1 as I compare only 'consecutive' polynomials. c) That was a typo from my side, now corrected. a) If the higher order polynomial provides a worse fit, its rss will be higher, which in turns implies the F-test can be negative? – CuriousMind Jul 4 '15 at 12:19
• Now it is a beauty! RSS will never increase if you add a parameter, so you are safe. – Michael M Jul 4 '15 at 13:57

You are right to be suspicious. That's because there is no way for the F-test to account for overfitting, so as long as the additional polynomial terms improve the in-sample fit even slightly, the null hypothesis of zero improvement will be false. In large data sets, you will very rarely fail to reject this hypothesis, even when the model is grossly overfitted.

By the way, an F test between models that differ by only one term is equivalent to a T test of the hypothesis that the one additional coefficient is zero. That is, you're just doing forward selection on an infinite set of increasingly irrelevant features. Even aside from the issue that you will probably overfit with exploratory polynomial regression, there are good reasons to avoid adding and keeping predictors based on one-at-a-time T testing.

If you genuinely believe that a polynomial is the appropriate model, Long and Trivedi (1991) provide a detailed and lucid review of some misspecification tests.

• Hi @ssdecontrol, thanks for your answer. I'm rather sure that a polynomial will do in my case, I'll check the document you suggest but if you have a pointer to a specific test it would be great. I'm kind of surprised there is no straightforward test to apply for this actually, but maybe I'm just underestimating the task. For example, another thing I tried was to simply look for convergence of the adjusted R$^2$, although 'convergence' was found empirically via some kind of threshold. This also occasionally led to high order polynomials with no visual difference from lower order ones. – CuriousMind Jul 4 '15 at 12:33
• R-squared, even if adjusted, can only ever give you information about in-sample fit. Out-of-fold R-squared in cross-validation will be much more useful – shadowtalker Jul 4 '15 at 13:00
• Although it's not the question, I'd amplify the warning here about overfitting. Polynomials up to degree 9, even when they fit the observed data well, have typically no interesting or useful substantive interpretation and will go haywire outside the observed range. If goodness of fit to data is the main concern, you will do as well or better with some more localised method, say spline or local polynomial. – Nick Cox Jul 4 '15 at 14:04
• Overfitting is exactly my concern, justified by the fact that I don't visually see a substantial improvement of polynomial of order 9 w.r.t polynomial of order 6 for example. I'm preparing a small test dataset to share and a figure to show this. Using the F-test formula above, I get that n=9 gives a statistically better fit than n=8 (and so on) with 99% confidence interval, but I don't trust this. – CuriousMind Jul 4 '15 at 14:23
• @ssdecontrol, I'm checking the LOOCV and see if I get different results than with the F-test above. Will update the discussion as soon as I have more info. – CuriousMind Jul 4 '15 at 14:57

I played around a bit with the data you posted, doing a 2-fold cross-validation with polynomial fits up to order 9. Someone else should check me on if I'm getting the CV procedure right, but I split the data into two groups; ran polynomial fits on each separately; and then took

$\sum_i \left[ (y^0_i - \hat{y}^1(x^0_i))^2 + (y^1_i - \hat{y}^0(x^1_i))^2 \right]$

as the residual sum of squares. That is, I took the polynomial fit from one group, plugged in the x-values of the other group, and compared it to the y-values of the other group; and vice versa. This quantity itself was minimized for an order 5 polynomial; and the RSS actually increased for the 3,6,7, and 9-degree polynomials. I haven't done a proper F-test comparing, say, order 5 to order 3, but it seems like order 5 should be good enough. Here are the results from a final fit of the whole set to order 5. I think it speaks to the cross-validation that the 5th-order fits to the whole dataset are similar to the fits from only half the data at a time. Just about every fit I looked at has some weird edge effects, though, so you do probably want to go with a spline fit or something to manage that behavior.

Edit: I revised how I was doing the validation, which should be more in line with the Wikipedia entry. I also started doing F-tests against the last best fit, rather than just the previous order polynomial. The results...n=9 shows up as the best order every time. That might actually be your best fit to this data?

I'm posting my code here (in Python) and attaching two new plots. The first one shows the RSS values from each fold and the mean RSS across folds, as a function of polynomial order. The second one just shows the data with the cross-validated polynomial fits for each order. The 'best' one, based on that new F-Test, is highlighted. Let me know if this helps!  • Hi @Jareth, thanks for your help. I tried this as well but order 9 is still chosen most of the times. I ran the same algorithm a number of times (the folds are generated randomly each time) and still get that as 'preferred' choice while order 5 is chosen very seldomly. The fact that a different polynomial is chosen (although 9 is chosen most often) let me think it depends a lot on how the folds are generated. Have you tried your algorithm multiple times? Does it provide a more stable solution? – CuriousMind Jul 5 '15 at 23:44
• Additionally, I also tried a k-fold with k=5. In this case I divide data into k subsets and I use the first 4 to produce 4 fits using polynomials from n=1 to 9. I then use these fits to generate data at the validation x values (the 5th subset) and calculate the RSS against the validation y values. I repeat this whole procedure k times, which means that for each polynomial I produce 20 RSS estimates. I minimize their sums to select the best polynomial. Still get number 9 most of the times. – CuriousMind Jul 6 '15 at 0:07
• I think this comes down to selection criteria. I did a few runs, and orders 5 and 6 were usually the minimum RSS values with orders 8 and 9 once each. In every case, the order 9 polynomial improves things over the order 8 polynomial based on an F-Test. I think what you need to do is generate an F-test comparing the current polynomial fit to the last best polynomial fit, rather than to the previous order. I also think that in the 5-fold validation, you produce 1 fit that you test on 4 different independent datasets. Can someone more experienced verify that? – Jareth Holt Jul 6 '15 at 14:43
• Hi, I based my algorithm on the description on the wikipedia (en.wikipedia.org/wiki/…), which of course might be wrong. It'd be great if anyone could confirm if the method is sound. I changed slightly the selection criterion for the k-fold algorithm and now orders other than 9 are also chosen (including 5 and 6 relatively often) although there is little consistency between different runs, which is worrying. – CuriousMind Jul 6 '15 at 21:28
• My 2-fold algorithm instead should be like yours (including the selection criterion), but it almost always selects order 9 as the solution (so there must be something different there which I haven't found yet). Are the subsets in your program completely disjointed? I noticed in my case there is about 30 non-unique elements over about 700 observations used as training sets. – CuriousMind Jul 6 '15 at 21:28