Is MSE decreasing with increasing number of explanatory variables?

Question

I am wondering, if there is a negative correlation between Mean Squared Error

\begin{equation} MSE = \frac{1}{n} \sum (\hat{Y}_i - Y_i)^2 \end{equation}

and the number of explanatory variables. Intuitively, I would guess that more independent variables should be able to explain more variation in my dependent variable. However, I was not able to find any literature about the topic.

Question 1: Is MSE decreasing with increasing number of explanatory variables? If yes, why?

Question 2: Is there any literature about the topic available?

Answer 1

I am assuming that you are talking about an ordinary least squares regression scenario and are referring to in-sample MSE, and that $Y$ is an n-by-1 vector and $X$ is an n-by-p matrix of orthogonal predictors (or variables, by your terminology). Remember that the columns of any matrix $X$ can be orthogonalized; this will become important for making an intuitive leap later on. Let's also assume that the columns of $X$ have variance $1/n$ and are centered, such that their means are zero.

Granted the foregoing, the answer to (1) is yes. Here's why.

MSE = $(1/n)\|Y-\hat{Y}\|^2$

$=(1/n)\|Y-Xβ\|^2$

$=(1/n)\|Y-X(X^TX)^{-1}X^TY\|^2$

Now, $(X^TX)^{-1}$ is simply a p-by-p identity matrix (this follows from the orthogonality we imposed earlier). It then follows that $(X^TX)^{-1}X^T=X^T$, and we have

MSE = $(1/n)\|Y-XX^TY\|^2$

So, what can we say about $XX^T$? We know it is an n-by-n matrix, and in the special case of p=n, it is an n-by-n identity matrix. That is, for p=n,

MSE = $(1/n)\|Y-Y\|^2 = 0$

Which we know to be intuitively correct. Furthermore, we know that MSE is at its maximum when we lack any predictors and $X$ is simply a column of ones; this is how we would fit an intercept-only model. In such a case, $XX^T$ is an n-by-n matrix of ones. As p gets larger, the off-diagonal elements of $XX^T$ shrink, eventually reaching zero when p=n.

This is not a rigorous proof and it does not, in fact, demonstrate that MSE is monotonically decreasing with p, but I think it provides a good intuitive foundation for understanding the behavior of least squares fitting.

Edit: If you want to extend this analysis to estimating MSE out of sample, then you would consider the following:

$\hat{MSE}=\hat{bias}^2+\hat{var}$

$\hat{bias}^2$ is monotonically decreasing with p, and $\hat{var}$ is monotonically increasing with p. There are some relationships between p, n, and $\hat{MSE}$, for that I recommend Wessel van Wieringen's lecture notes on ridge regression as well as Elements of Statistical Learning, as mentioned in another answer to your original question. Hopefully that answers (2).

Edit: I thought about this some more and is are two additional points I'd like to make. The first is the specific conditions under which an additional predictor will reduce in-sample MSE. Those conditions are:

1) The additional predictor does not lie entirely within the column space of $X$; that is, it cannot be obtained via any linear combination of the existing predictors, and

2) The component of the new predictor lying outside the column space of $X$ is not orthogonal to $Y$.

The second point is that we can do a simple thought experiment showing that the addition of new predictors does, in general, tend to decrease in-sample MSE. Imagine we have solved our linear regression and obtained $β$, that is, a p-by-1 vector of model coefficients. Now imagine that we add an additional predictor. Unless BOTH of the two aforementioned conditions are satisfied, the (p+1) value of $β$ will be zero, and the model is exactly the same as it was prior to the addition of the new predictor (same MSE). In general, though, both of those conditions will be satisfied, and therefore the (p+1) value of $β$ will be something other than zero. Since both the zero-appended $β$ and nonzero-appended $β$ lie within the solution space of the least squares regression with p+1 predictors, we conclude that the p+1 model must have lower MSE than the p model if the new coefficient is anything other than zero.

Answer 2

It is unfortunate that Empirical Risk Minimization is a topic where the Internet is full of incorrect information including https://en.wikipedia.org/wiki/Mean_squared_error, where Wikipedia changes its mind whether it is an estimate or a population parameter within the same article.

First of all MSE is a population metric, therefore it is incorrect to talk about in-sample or out-of-sample MSE.

Let's assume we have $N$ $(y,X)$ pairs for training where $y$ is a sacalar and $X$ is a vector containing predictors. Let's also assume that $\hat{y}_N(X)$ is a predictor found using a fitting procedure using those $N$ samples, then $MSE$ of prediction of fitted model is:

$MSE(\hat{Y}) = E[(\hat{y}(X)-y)^2] = \int (\hat{y}(X)-y)^2 p(y,X) dydX$

On the other hand we can have an estimate of MSE. There are two flavors of this: the in sample and out-of-sample estimate. Please take note of the hat I am putting over $MSE$ below.

$\hat{MSE_{in}} = \frac{1}{N}\sum_i^N ((\hat{y}(X_i)-y_i)^2$

Let's assume we have another $M$ test samples we haven't used during the fitting procedure:

$\hat{MSE_{out}} = \frac{1}{M}\sum_i^M ((\hat{y}(X_i)-y_i)^2$

Now what can we say about the expected values of these estimators?

It is easy to see that

$E[\hat{MSE_{out}}] = MSE$

simply push the expectation inside the summation and use that fact that $\hat{y}$ is independent of the test samples. Also we use the $IID$ assumption over the samples.

Same cannot be concluded for in sample estimate because $\hat{y}$ is not independent of the training samples (obviously as we used the training samples to find $\hat{y}$), hence:

$E[\hat{MSE_{in}}] \ne MSE$

In fact it is downward biased (https://web.stanford.edu/~hastie/Papers/ESLII.pdf, section 7.4):

$E[\hat{MSE_{out}}] = E[\hat{MSE_{in}}] + \frac{2}{N}\sum_1^N Cov(y_i,\hat{y_i})$

This is the source of the nasty overfitting problem.

Having all that cleared now let's look at how these estimates change according to number of predictors and number of samples used in the fitting process

The difference between the red and blue lines illustrates the bias between the in-sample and out-of-sample estimates of MSE. We also added a new concept: "MSE of proposed model". This is not the MSE of the fit but the asymptotic MSE as if we had infinite training samples. As the number of training samples increase both in-sample estimate and out-of-sample estimate converge to the flat line given that the training procedure is consistent!!! (http://bengio.abracadoudou.com/lectures/theory.pdf, slide 10)

Now let's examine the effect of adding predictors.

For nested models adding predictors always decreases (at least doesn't increase) $E[\hat{MSE_{in}}]$ given that the fitting procedure reached the global minimum. For nonlinear models optimization usually gets stuck in local minimum, and this may appear to have increased the $\hat{MSE_{in}}$, but this is not an effect of the added predictor, but a side-effect of the fitting procedure.

For non-nested cases it can go either way. I believe the original question was towards a nested linear model.

What happens to $\hat{MSE_{out}}$ is more interesting and deserves a demonstration. It is also more important as it provides an unbiased estimate of $MSE$ and is also a key objective in Cross Validation procedures.

Suppose that we have data generated by a polynomial of degree 4 with additive Gaussian noise. Below is a bunch of samples to illustrate how it looks like.

Now I will plot $\hat{MSE_{out}}$ for fitted polynomials from degree 2 to degree 6. vertical axis is logarithm of error and horizontal axis is again number of samples used in fitting procedure. I always used 1000 test samples to estimate $\hat{MSE_{out}}$.

The correct model is degree 4 and plotted in red color. magenta is degree 2, green degree 3, yellow degree 5, and blue degree 6. As you can see when the number of training samples is less than 10, the best $MSE$ is achived with a polynomial degree 2, after that polynomial degree 3 takes over, and when we reach 15, degree 4 dominates and from that point on it becomes unbeatable. This is a known phenomenon where true model degree may not give the best prediction ability on small samples.

When will a less true model predict better than a truer model?

Another interesting observation is that even though degree 2 and degree 3 do a better job at small samples, as number of training samples increase not only they are dominated by degree 4 (the true model), but also degree 5 and 6. This is natural as degree 5 and degree 6 subsume degree 4.

Answer 3

The short answers are:

1. Yes. A more precise answer should be "non-increasing", as mentioned in comment. For example, if we include a complete random noise as independent variable, it will not make MSE decrease, but make MSE the same.
2. See the links mentioned below.

Assuming we only have one data set, and want to build a model to have as low MSE as possible. Then, we will have a over-fitting problem. Basically we the model we build are too specific on this given data set and the model is not able to generalize.

People also call this as a "false high R square problem", see this post for reasons why people use adjusted R square. Why is adjusted R-squared less than R-squared if adjusted R-squared predicts the model better?

The intuition behind adjusted R square or other regularization algorithms are trying to consider both performance (such as MSE) and number of parameters used int he model. As a result, although increasing number of independent variables will not hurt MSE, but if the improvements are "marginal", the recommendations would still be not adding it.

---

Answer 4

Ironically, the answer begins far away in a nether world of philosophy.

Error is rarely intrinsically an "error" at all, except in quantum physics. It is the result of a virtually infinite number of factors, unobserved, contributing to an almost chaotic effect. If we had a correctly specified, infinite dimensional model, there would be no error in deterministic systems, like weather patterns, stock prices, planetary motion, and so on. We would have effects for the proverbial butterfly that flaps its wings in China, causing the Tsunami in Guam. And if we reliably estimate a high dimensional, though finite, model, the "errors" are still smaller than would be seen in a lower dimensional model.

So the answer to 1 is yes, in general a well-estimated, high dimensional models confer better predictiveness, and thus a lower MSE, than a model with fewer predictors.

The major caveat is reliability. Pitifully often, people use internal validation to check their model. Mean squared error is a quantity which we estimate using split-sample validation (good), cross-validation (good), bootstrapping (good), or internal validation (bad). In internal validation we use residuals which are, as a consequence of the model fitting procedure, perfectly orthogonal to the design matrix, to estimate error. This is problematic because any external, independent dataset will not have that added benefit, in that the linear combinations which comprise predicted values optimally reduce the error in that dataset. We would have to refit the model in every new instance of data, making the whole purpose of prediction a moot point.

Coincidentally, if you overfit a model then use internal validation: it is true that including more features will lead to a lower estimate of MSE, but in fact, the MSE of the model (as determined by validation in new, independent data) will be higher. This is because overfitting is out-of-sample variance, and the MSE is the sum of variance and squared bias.

The answer to 2 is myriad. But regarding prediction and especially the point above about overfitting, and better understanding bias/variance tradeoff, reliability, and prediction: an accessible (free) text is Elements of Statistical Learning.