Selecting the best regression model using R2 and AIC - what is the best approach?

Question

I have a dataset in which I have one dependent and 3 independent variables (y ~ x1 + x2 + x3). For exploratory analysis, I have fitted the following models (using R):

mod1 <- lm(y ~ x1 + x2 + x3, data = data)
mod2 <- lm(y ~ x1, data = data)
mod3 <- lm(y ~ x2, data = data)
mod4 <- lm(y ~ x3, data = data)
mod5 <- lm(log(y) ~ x1, data = data)
mod6 <- lm(y ~ log(x1), data = data)

#Polynomial regression
mod7 <- lm(y ~ poly(x1, 1, raw = TRUE), data = data) #Fits a first order regression mod
mod8 <- lm(y ~ poly(x1, 2, raw = TRUE), data = data) #Fits a second order regression mod
mod9 <- lm(y ~ poly(x1, 3, raw = TRUE), data = data) #Fits a third order regression mod
mod10 <- lm(y ~ poly(x1, 4, raw = TRUE), data = data) #Fits a 4th order regression mod
mod11 <- lm(y ~ poly(x1, 5, raw = TRUE), data = data) #Fits a 5th order regression mod
mod12 <- lm(y ~ poly(x1, 6, raw = TRUE), data = data) #Fits a 6th order regression mod

mod13 <- lm(y ~ poly(x2, 1, raw = TRUE), data = data) #Fits a first order regression mod
mod14 <- lm(y ~ poly(x2, 2, raw = TRUE), data = data) #Fits a second order regression mod
mod15 <- lm(y ~ poly(x2, 3, raw = TRUE), data = data) #Fits a third order regression mod
mod16 <- lm(y ~ poly(x2, 4, raw = TRUE), data = data) #Fits a 4th order regression mod
mod17 <- lm(y ~ poly(x2, 5, raw = TRUE), data = data) #Fits a 5th order regression mod
mod18 <- lm(y ~ poly(x2, 6, raw = TRUE), data = data) #Fits a 6th order regression mod

mod19 <- lm(y ~ poly(x3, 1, raw = TRUE), data = data) #Fits a first order regression mod
mod20 <- lm(y ~ poly(x3, 2, raw = TRUE), data = data) #Fits a second order regression mod
mod21 <- lm(y ~ poly(x3, 3, raw = TRUE), data = data) #Fits a third order regression mod
mod22 <- lm(y ~ poly(x3, 4, raw = TRUE), data = data) #Fits a 4th order regression mod
mod23 <- lm(y ~ poly(x3, 5, raw = TRUE), data = data) #Fits a 5th order regression mod
mod24 <- lm(y ~ poly(x3, 6, raw = TRUE), data = data) #Fits a 6th order regression mod

#Quadratic and cubic on the first independent variable
mod25 <- lm(y ~ x1 + x1^2, data = data)
mod26 <- lm(y ~ x1 + x1^2 + x1^3, data = data)

The resulting R2 and AIC values are displayed below:

       df        AIC         R2
mod1   5     -25.591466     0.67
mod2   3     -13.471396     0.02
mod3   3     -22.552733     0.46 
mod4   3     -17.774707     0.26
mod5   3     -3.854236      0.01
mod6   3     -13.180926     0.00
mod7   3     -13.471396     0.01
mod8   4     -12.613856     0.09
mod9   5     -10.876091     0.10  
mod10  6     -9.063492      0.11
mod11  7     -7.159273      0.12
mod12  7     -7.159273      0.12
mod13  3     -22.552733     0.46
mod14  4     -36.891799     0.82
mod15  5     -41.786514     0.88
mod16  6     -44.197137     0.92
mod17  7     -46.293560     0.94
mod18  8     -47.957872     0.95
mod19  3     -17.774707     0.26
mod20  4     -16.409778     0.29 
mod21  5     -15.417082     0.34
mod22  6     -18.144996     0.52
mod23  7     -23.829747     0.77  
mod24  8     -28.124779     0.81
mod25  4     -12.613856     0.09
mod26  5     -10.876091     0.11

Judging by the lowest AIC and highest R2 scores, the higher-order regression models seem to be performing the best. However, I'm not sure that just selecting the highest-order polynomial regression model is the best choice due to over-fitting a model.

Mod 1 seems to be the best fit without over-fitting, being that it is a multiple regression model, has a R2 value of 0.67, and an overall p-value of 0.0058. Does this reasoning seem sound? I would appreciate any insights/thoughts on this matter, as well as on regression model selection as a whole. Thanks!

Answer 1

Your concern about over-fitting is valid. But you can't check for overfitting by checking the performance of the model on the training set. Overfitting happens when the fit to the training set is better than the fit to a held out test set.

If you have enough data, set some aside as an independent test set. If you don't, then the best you can do is cross-validation.

Do you have any a priori grounds to expect one of these models to be better? It is pretty unusual for a real mechanism to involve a high order polynomial.