# Choosing regression model

How can I choose between these two regression models?

R outputs:

Regression 1

Call:
lm(formula = log.h ~ year)

Residuals:
Min       1Q   Median       3Q      Max
-1.24004 -0.45221 -0.05301  0.44744  1.42060

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.588e+02  1.565e+01  -10.15 9.47e-13 ***
year         8.809e-02  7.859e-03   11.21 4.65e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6396 on 41 degrees of freedom
Multiple R-squared: 0.7539, Adjusted R-squared: 0.7479
F-statistic: 125.6 on 1 and 41 DF,  p-value: 4.654e-14  Regression 2:

Call:
lm(formula = log.h ~ year + I(year^2) + I(year^3))

Residuals:
Min      1Q  Median      3Q     Max
-1.4044 -0.3966  0.1141  0.3447  1.2321

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.075e+06  4.663e+05  -2.304   0.0266 *
year         1.622e+03  7.027e+02   2.308   0.0264 *
I(year^2)   -8.162e-01  3.529e-01  -2.313   0.0261 *
I(year^3)    1.369e-04  5.909e-05   2.317   0.0259 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.58 on 39 degrees of freedom
Multiple R-squared: 0.8075, Adjusted R-squared: 0.7927
F-statistic: 54.54 on 3 and 39 DF,  p-value: 5.093e-14  Thoughts: The point is the first model is good, $R^{2}$ is good, p-values are great, but there are hints of heteroskedasticity of residuals vs predicted values.

The second one has a better adjusted $R^{2}$, a great $p$-value of regression, and the scatterplots have got better (not drastically, but better). On the other hand, the $p$-values of the coefficients are a lot higher than before.

Considering all this, how can I choose between them?

• Skimming your plots, @PeterFlom is right: there is clearly some structure in the residuals that isn't being captured by the 1st model fit. The 2nd is better, but isn't quite perfect. Using a spline seems ideal here. Jul 10, 2013 at 22:36
• whatever your variable h is,it doesn't make sense adding square and cubic of time variable without any theoretical justification. Jul 10, 2013 at 22:45
• I would not interpret the pattern of residuals in the first model as heteroscedasticity, but instead as a possible non-linear trend. Applying a robust smoother (such as LOESS) to the residuals will help you see that. If the smoother is close to a straight horizontal line you can suspect heteroscedasticity (and test it by smoothing the roots of the absolute values of the residuals), but if it wiggles around appreciably, you have detected a systematic variation in the average residual.
– whuber
Jul 11, 2013 at 1:59
• the lowess is't close to a horizontal line untile 1990, then bent of about 40 degrees!
– Ant
Jul 11, 2013 at 9:55
• You often need to tune the degree of smoothing with lowess (loess). (Some people apply cross-validation; some programs do that automatically, or so I guess.) Jul 11, 2013 at 10:11

I agree with the general tenor of other comments (including an aversion to stepwise). Like others I have reservations about polynomials here, and would tend to prefer splines, but you can do one simple thing to make your polynomial results less unfriendly.

The one plot most missing is a basic plot of the data, as @Peter Flom emphasized, but some reverse engineering makes it clearer what is happening.

Your first model is evidently close to

$-158.8 + 0.08809\ {\rm year}$,

and from your predicted values it can be guessed that year ranges from about 1975 to 2010.

Your second model is evidently close to

$-1075000 + 1622\ {\rm year} -0.8162\ {\rm year}^2 + 0.0001369\ {\rm year}^3$,

although digits not shown are essential to get predicted values close to those shown.

These polynomials are close numerically but not algebraically. A simple trick that would make it easier to work with polynomials is to choose an origin within the range of the data. There is no reason to privilege the origin of the Western calendar. Working with (say) (year $- 2000$) would give the intercept the simple interpretation of value in $2000$ and should make it easier to compare linear, quadratic and cubic fits.

• I would really like to thank everyone, every answer has improved my understanding of regression.. I was mislead by the high p-value values, which, as i now discover, are not that significant (since translating -2000 years solves the problem of p-values, but of course the regression it's exactly the same)
– Ant
Jul 11, 2013 at 9:59
• I appreciate the acceptance, but want to underline that other answers contain excellent strategic and tactical advice. Jul 11, 2013 at 10:10
• that is why I upvoted all of them :-)
– Ant
Jul 11, 2013 at 10:33

First, I wouldn't judge solely by $R^2$ and not at all by p values

Second, like @gung I would not use stepwise or its variants.

Third, your plots are nice, but the first plot I'd look at is a simple scatterplot of log.h and year, perhaps with a spline or other smooth added. Something like

plot(log.h~year)
lines(lowess(log.h~year))


to see what I could see.

Fourth, I would think about what the relationship could mean. Since I don't know what log.h is, I have no idea. But does accelerating change over time make sense? Or decelerating?

Fifth, I'd be careful of overfitting the sample, and especially of extrapolating beyond it with a polynomial model, which is often the aim when year is a predictor.

Sixth, I'd consider a spline fit.

• Indeed, it makes quite sense that it's accelerating over time with a greater speed than linear one.. I guess I'm going for the cubic fit
– Ant
Jul 11, 2013 at 9:57

Three things:

1. The adjusted $R^2$ value of your 2nd model is larger than your first. This means it explains more of the variation in the data than the first one. Unlike the standard $R^2$ value it is not impacted by the number of variables added to the model.

2. To more rigorously compare them you can run an ANOVA test like so:

anova(regression1.model, regression2.model)


Look for the p-value. The closer to zero the better, though common practice is it should be at least less than .05.

3. Make sure it makes sense. If you created regression 2 because you thought it might be a better fit and the above other points are met, you've found yourself an improved model.

This is a hard question to answer based on what you've shown. And that's why it's hard for you. If you had more information included it would be easier.

One might question whether the polynomial is contributing a meaningful amount of extra predictive power given the cost of the additional parameters and assess using AIC (AIC function in R). It doesn't seem like they add much (they'll always add some) but you can check.

Then there's the question of whether the curved wavy line you've made is overfit. Will it predict new values? You could look up cross validation techniques but for something like that.

Probably the more important thing to do is assess what you're expecting your predictor to do at a theoretical level. Should it have a linear effect? Are you expecting it to be wavy (you polynomial)? Why is it wavy? If you can't address these and other theoretical questions then neither of these statistical models is particularly good or bad.

To help you on the use of AIC you might want to look at a couple of papers.

$R^{2}$ will always increase as you add more variables to the model. It seems what you are asking is what variables to include in the model. This is called variable selection. You can accomplish this in R by the step function. This function bases the best model on a statistic called AIC. Make sure you use the full model, the one with the most variables.

step(yourfullmodel)


the final output will give you your the best subset model under Call:

• Welcome to the site, @willy. Note that stepwise variable selection routines are not considered valid anymore. If this doesn't make sense to you / you want to know why, it may help to read my answer here: algorithms-for-automatic-model-selection. Jul 10, 2013 at 21:59
• He does have the adjusted R^2 however which does not suffer from the same flaws as the standard R^2. Jul 11, 2013 at 1:51