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.
 A: 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.
A: 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!


