# Why do we even bother running regression models?

I'm working through regression with Intro to Statistical Learning by Hastie, Witten, James and Tibshirani. They break down regression into stages: data cleaning and processing, model building and selection and model testing.

There are many pieces to this where the decision is algorithmic. I look at two models and pick the one with the better adjusted r squared or, I look at collinearity and remove values with really high VIF scores.

Given the sheer amount of knowledge we have on best practices for each of these steps and the if-this-then-that nature of running regression analysis, why are we doing them at all?

Are there programs that will make every possible model and give you a list of the best models to choose from? It could try splines and polynomials and so on.

• There won't be a single model that optimizes across all possible scoring criteria, so there is always going to be a need for researcher intervention. Selecting a single criteria (e.g. R-squared) guarantees meaningless, a-theoretical models that have as many co-variates as observations. Or, suppose you select on p-values for some coefficient; then you have just automated p-hacking. The whole point is that theory should underpin model selection – Parseltongue Jan 18 at 20:02
• Unless you have a truly vast number of data relative to your variables & the modeling options, your suggestion will be too opportunistic & lead to models that don't fare well in the long run. It may help to read my answer here: Algorithms for automatic model selection. – gung Jan 18 at 20:04
• As a counterpoint to that example, please read some of our posts on "overfitting" and perhaps some related to "regularization" or "Lasso." – whuber Jan 18 at 20:45
• Cross validating every possible model or a selection out of some standard repertoire comes at a cost. You are better of when some experimenter with (good) prior knowledge selects some more sensible model. It will reduce the variance/error in the estimates. If you 'know' that a linear model is working then you will be better of with fitting a linear model rather than having the computer spit out some higher order polynomial that happens to be a better prediction for your particular experiment. – Martijn Weterings Jan 18 at 23:31
• If you run many models, regularized or other, and choose the one that performs best on your test sample, you will overfit to the test sample. – Stephan Kolassa Jan 19 at 8:29

Are there programs that will make every possible model and give you a list of the best models to choose from? It could try splines and polynomials and so on.

...

My understanding is you can avoid overfitting by estimating test error with techniques like cross validation.

Cross validating every possible model or a selection out of some standard repertoire comes at a cost. You are better of when some experimenter with (good) prior knowledge selects some more sensible model. It will reduce the variance/error in the estimates. If you 'know' that a linear model is working then you will be better of with fitting a linear model rather than having the computer spit out some higher order polynomial that happens to be a better prediction for your particular experiment.

The main use of cross-validation is to test the 'performance' of a model. The performance is influenced by two factors (1) the bias in the deterministic part of the model (2) the error in the estimate of unknown model coefficients. (related to the balance of bias and variance is this question about the use of the loss function).

Cross-validation is providing a scheme where you make a balance between the bias of the model and the variance of the model (simply by some sort of trial and error). Improvements are made when

1. you have an unbiased model (possibly with the lowest possible variance out of all possible unbiased models) but CV finds a biased model that makes up with lower variance,
2. or when you have a biased model (or no clue) and CV find a less biased model and/or less variable model.

But this is not necessarily the best balance.

• CV can increase overfitting By introducing all those alternative models, you may on the one hand find a model with less bias (especially usefull when you do not know much about the bias of your model or models), but on the other hand you have a higher tendency to overfit because of the higher degrees of freedom.
• CV can decrease overfitting Cross-validation may reduce the tendency to overfit because it includes biased models that may have lower variation in their model parameters. Strictly speaking you do not need to use cross-validation to get this lower variation; one could use some lower variance parameter estimate straight away (e.g. all sorts of shrinkage estimators), but when you do not know by how much the parameters need to shrink then you can use cross-validation.

Based on prior knowledge, an experimenter may decide to not use cross-validation. For instance when he knows that there is not gonna be much bias. Also, cross-validation is not necessary when the aim of the experiment is to estimate parameters or to test a theory.

So cross-validation and regression are used in different settings. You will find the cross-validation more often in settings where the bias of models is being tested. Either one wishes to remove bias when there is no good theoretical idea about the right model and when one would expect bias in the deterministic part of the model, or alternatively one actually has a low bias model but one wishes to introduce bias in order to decrease the variation by means of some sort of parameter shrinking.

Cross validation is more like a model selection tool than a model fitting tool. Although in machine learning these two may become blurred, for instance 'training a neural network' might be seen as both model fitting and model selection.

### Example

Example (inspired by an answer here) illustrating the above two points:

Let the true model be some 4-order polynomial with Gaussian noise $$\epsilon \sim N(0,\sigma_\epsilon^2)$$:

$$y = 0 + 1 x + 0.5 x^2 + 0.1 x^3 + 0.05 x^4 + \epsilon$$

We fit three type of models

• a biased 2nd order polynomial with gaussian noise
• a unbiased 4th order polynomial with gaussian noise
• a cross validation of 10 models (polynomials of order 1 to 10)

we do this 2000 times (while varying the noise level $$\sigma_\epsilon$$) by leave one out cross validation (which is, in fact, just equivalent to regression) with training data of sample size 20.

The validation of the 2000 times 3 models is done with new data of sample size 1000. The image below displays the mean of the squared residuals in the validation data scaled by the known variance $$\sigma_\epsilon^2$$.

The path of individual curves

• The unbiased model (4th order polynomial) has a relative constant error. This error consist of the error in the validation data and the error in the fitted model parameters. These scale both similarly with $$\sigma_\epsilon$$.
• The biased model (2th order polynomial) does not not have a relative constant error and initially the error is larger. This is because the 2th order polynomial is biased and this bias weighs more strongly when the error terms $$\epsilon$$ are smaller.
• The cross validation model is a bit in between the two biased and unbiased models.

Pairwise comparisons of the curves

• biased vs unbiased: When $$\sigma_\epsilon$$ is small then the unbiased curve is below the biased curve (the unbiased model performs better). This is because the bias adds a relatively larger error. When $$\sigma_\epsilon$$ is large then the biased curve is below the unbiased curve (the biased model performs better). This is because the biased model has a smaller variance of the model coefficients. It will be relatively less over-fitting (modeling the noise rather than the true model).

• unbiased vs cross validated: When $$\sigma_\epsilon$$ is small then the unbiased curve performs is below the cross validated curve (the unbiased model performs better). This is because the cross validated model is more likely to introduce additional variance in the fitted model by having all those extra biased models (which you might consider to be wrong and superfluous in this setting). When $$\sigma_\epsilon$$ is large then the cross validated curve is below the unbiased curve (the cross validated model performs better). This is because the cross validated model will be having a smaller variance by introducing the biased models with smaller variance (which you might consider still wrong yet useful in this setting).

• biased vs cross validated: When $$\sigma_\epsilon$$ is small then the biased curve is above the cross validated curve (the biased model performs worse). This is because the cross validated model is more likely to introduce some model with less bias. When $$\sigma_\epsilon$$ is large then the biased curve is below the cross validated curve (the biased model performs better). This is because the biased model does well in reducing the variance of coefficients and the cross validated model (having more degrees of freedom) will be having a larger variance.

In this last case you might devise some examples where the cross validated model performs better than the biased model at low $$\sigma_\epsilon$$. Then upon increasing the errors $$\sigma_\epsilon$$ the biased model does better, but after increasing further the cross validated model does better again.

So the comparison 'unbiased vs cross validated' relates to using cross validation to introduce bias that reduces the variance (setting: finding a model with better prediction performance), and the comparison 'biased vs cross validated' relates to using cross validation to reduce the bias (setting: finding a model that is more correct to the true model). In both cases there is a good side and a bad side.

                    |                           MODEL
|
|        biased               |   unbiased
|        or unknown bias      |
|                             |
----------------------------------------------------------------------
Cross-validation    |     - reduce bias           |   - introduce bias
reason          |     - to select model       |   - to reduce variance
|       with low bias         |     and reduce prediction error
----------------------------------------------------------------------
NOISE LEVEL       |
|       CV better                  CV worse
low            |       due to selecting           due to more variance
|       lower bias                 (more degrees of freedom)
|
----------------------------------------------------------------------
|
high           |       CV mixed                   CV better
|       more variance              due to selecting
|       (more degrees of freedom)  model with lower variance
|       less variance
|       (better model)


### Baysian viewpoint

I am sure that the example above could be formulated as a Bayesian problem where you have some prior distribution on the model parameters and after gathering data you have a posterior distribution on the model parameters.

Then the use of cross-validation may be seen as selecting that particular prior and you may wonder whether that selection (e.g. assigning probability to higher order polynomials) is 'correct'.

### Code example

Code example to produce the above images in R.

set.seed(1)

nn <- 2000
nt <- 20
nv <- 1000
np <- 10

perfb <- rep(NA,nn)
perf2 <- rep(NA,nn)
perf3 <- rep(NA,nn)
perf4 <- rep(NA,nn)
sig <- 2^seq(0,12,length.out=nn)
bestcvs <- rep(0,np)

for (i in 1:nn) {
# training and test data (using leave one out)
x <- rnorm(nt,0,4)
y <- 0 + 1*x + 0.5*x^2 + 0.1*x^3 + 0.05*x^4 + rnorm(nt,0,sig[i])

# training and testing polynomials from 1 to 20
meanerrors <- rep(NA,np)
for (j in 1:np) {
fitt <- lm(y~1+stats::poly(x,j,raw=TRUE))
meanerrors[j] <- loocv(fitt)
}

# polynomial with lowest testing error
best <- which.min(meanerrors)
bestcvs[best] <- bestcvs[best]+1

fittb <- lm(y~1+stats::poly(x,best,raw=TRUE))
fitt2 <- lm(y~1+stats::poly(x,2,raw=TRUE))
fitt3 <- lm(y~1+stats::poly(x,3,raw=TRUE))
fitt4 <- lm(y~1+stats::poly(x,4,raw=TRUE))

# validation data
x <- rnorm(nv,0,4)
y <- 0 + 1*x + 0.5*x^2 + 0.1*x^3 + 0.05*x^4 + rnorm(nv,0,sig[i])

# predictions from models
Xf <- cbind(rep(1,40),as.matrix(stats::poly(x,20,raw=TRUE)))
ypb <- Xf[,1:(best+1)] %*% as.numeric(coef(fittb))
yp2 <- Xf[,1:(3)] %*% as.numeric(coef(fitt2))
yp3 <- Xf[,1:(4)] %*% as.numeric(coef(fitt3))
yp4 <- Xf[,1:(5)] %*% as.numeric(coef(fitt4))

# validation performance
perfb[i] <- mean((y-ypb)^2)
perf2[i] <- mean((y-yp2)^2)
perf3[i] <- mean((y-yp3)^2)
perf4[i] <- mean((y-yp4)^2)
}

# plot
plot(sig,(perfb)/sig^2,log="xy",  pch=21,col=rgb(0,0,0,0.05), bg=rgb(0,0,0,0.05),
xlab = expression(paste(sigma[epsilon])),
ylab = expression(paste(SSE/n/sigma[epsilon]^2)),
ylim=c(0.8,10^6))
# plotting smoothened lines
lines(sig,exp(sgolayfilt(log(perfb/sig^2),p=1,n=201)),col=rgb(0,0,0))
lines(sig,exp(sgolayfilt(log(perf2/sig^2),p=1,n=201)),col=rgb(1,0,0))
lines(sig,exp(sgolayfilt(log(perf4/sig^2),p=1,n=201)),col=rgb(0,0,1))
title("performance cross validated model")
plot(sig,(perf2)/sig^2,log="xy",  pch=21,col=rgb(1,0,0,0.05), bg=rgb(1,0,0,0.05),
xlab = expression(paste(sigma[epsilon])),
ylab = expression(paste(SSE/n/sigma[epsilon]^2)),
ylim=c(0.8,10^6))

• Very interesting and insightful post! I just wonder if a noise level of 5000 is making a lot of sense since the signal level is only about 4? – Vickyyy Jan 22 at 18:23
• That is right. With high noise levels you are basically making good predictions by just using the mean of the sample since the higher order models are not gonna make a difference. This extreme example is just for the sake of making the illustration of the cases more clear. The principle remains in less extreme cases. – Martijn Weterings Jan 22 at 18:43