# Is Fig 3.6 in Elements of Statistical Learning correct?

Here is the figure from the textbook:

It shows a decreasing relationship between subset size $$k$$ and mean squared error (MSE) of the true parameters, $$\beta$$ and the estimates $$\hat{\beta}(k)$$. Clearly, this shouldn't be the case - adding more variables to a linear model doesn't imply better estimates of the true parameters. What adding more variables does imply is a lower training error, i.e. lower residual sum of squares.

Is the $$y$$-axis labelled incorrectly? In particular, is it possible that the $$y$$ axis shows e.g. Residual Sum of Squares instead of $$\mathbb{E}|| \hat{\beta}(k) - \beta||^2$$?

EDIT:

Discussions and multiple attempts to reproduce revealed the axis is likely labelled correctly. In particular, it is not RSS since that will be on a completely different scale.

The title question still remains - "Is Figure 3.6 in ESL correct?". My intuition is that MSE should be lowest around the optimal $$k$$ (@SextusEmpiricus's answer suggests that's the case but there correlation is lower). Eyeballing Fig 3.6 we see MSE continues to go down beyond $$k=10$$.

In particular, I'm expecting to see curves similar to those in Figure 3.16:

It does show additional procedures due to that is on a different $$x$$-axis; it also uses different number of samples (300 vs 100). What is relevant here is the shape of e.g. "Forward stepwise" (common in both charts - orange in the first, black in the second) which exhibits quite different behaviour across the two figures.

Final Edit

Here you can find my attempt at replicating Fig3.6; plot shows different levels of correlation and number of non-zero parameters. Source code here.

• Unrelated (to the question), how would the SNR have been calculated? Nov 16 '20 at 11:07
• @SingleMalt the SNR has been calculated as 10*0.4/6.25 (but whether that is a correct figure is another thing). Nov 16 '20 at 11:39
• @SextusEmpiricus got it. The authors could have avoided noting the SNR usage in the figure label, as the “noise” from the 21 extra correlated variables seems more important to one of the points they are trying to show here, which seems to be OVB. Nov 16 '20 at 12:02
• @SingleMalt the omitted variable bias has a very weird effect in this case. Due to the high correlation the overall bias is effectively not so high. When the variable is included it will overestimate (because it needs to compensate the ommision of the others) when the variable is omitted it will be underestimated (because it is zero when it should be not). So in this case the ommision of variables is more like causing higher variance of the error than causing higher bias of the error. Nov 16 '20 at 12:12
• @SingleMalt the problem with this case is that the error due to the omitted variables can go both ways. So there is more error due to the ommision of the variables, but the average error (bias) is not so greatly affected due to the ommision (because the error can go both ways positive/negative). It is the variance of the error that mostly increases. Nov 16 '20 at 13:17

It shows a decreasing relationship between subset size $$k$$ and mean squared error (MSE) of the true parameters, $$\beta$$ and the estimates $$\hat{\beta}(k)$$.

The plot shows the results of alternative subset selection methods. The image caption explains the experimental design: there are 10 elements of $$\beta$$ which are nonzero. The remaining 21 elements are zero. The ideal subset selection method will correctly report which $$\beta$$ are nonzero and which $$\beta$$ are zero; in other words, no features are incorrectly included, and no features are incorrectly excluded.

Omitted variable bias occurs when one or more features in the data generating process is omitted. Biased parameter estimates have expected values which do not equal their true values (this is the definition of bias), so the choice to plot $$\mathbb{E}\|\beta -\hat{\beta}(k) \|^2$$ makes sense. (Note that the definition of bias does not exactly coincide with this experimental setting because $$\beta$$ is also random.) In other words, the plot shows you how incorrect estimates are for various $$k$$ for various subset selection methods. When $$k$$ is too small (in this case, when $$k<10$$) the parameter estimates are biased, which is why the graph shows large values of $$\mathbb{E}\|\beta -\hat{\beta}(k) \|^2$$for small $$k$$.

Clearly, this shouldn't be the case - adding more variables to a linear model doesn't imply better estimates of the true parameters.

Fortunately, that's not what the plot shows. Instead, the plot shows that employing subset selection methods can produce correct or incorrect results depending on the choice of $$k$$.

However, this plot does show a special case when adding additional features does improve the parameter estimates. If one builds a model that exhibits omitted variable bias, then the model which includes those variables will achieve a lower estimation error of the parameters because omitted variable bias is not present.

What adding more variables does imply is a lower training error, i.e. lower residual sum of squares.

You're confusing the demonstration in this passage with an alternative which does not employ subset selection. In general, estimating a regression with a larger basis decreases the residual error as measured using the training data; that's not what's happening here.

Is the $$y$$-axis labelled incorrectly? In particular, is it possible that the $$y$$ axis shows Residual Sum of Squares instead of $$\mathbb{E}\|\beta -\hat{\beta}(k) \|^2$$?

I don't think so; the line of reasoning posited in the original post does not itself establish that the label is incorrect. Sextus' experiments find a similar pattern; it's not identical, but the shape of the curve is similar enough.

As an aside, I think that since this plot displays empirical results from an experiment, it would be clearer to write out the estimator used for the expectation, per Cagdas Ozgenc's suggestion.

Is Figure 3.6 in ESL correct?

The only definitive way to answer this question is to obtain the code used to generate the graph. The code is not publicly available or distributed by the authors.

Without access to the code used in the procedure, it's always possible that there was some mistake in labeling the graph, or in the scale/location of the data or coefficients; the fact that Sextus has had problems recreating the graph using the procedure described in the caption provides some circumstantial evidence that the caption might not be completely accurate. One might argue that these reproducibility problems support a hypothesis that the labels themselves or the graphed points may be incorrect. On the other hand, it's possible that the description is incorrect but the label itself is correct nonetheless.

A different edition of the book publishes a different image. But the existence of a different image does not imply that either one is correct.

• It does shows bias. Note that the y-label is NOT MSE. It's the average difference between the estimated and the true parameters. Usually you don't know the true parameters, but here the data is simulated so the underlying ground truth is known. Nov 15 '20 at 17:35
• I would say this chart is very misleading for a student as it doesn't make it clear regarding with respect to what the expectation is taken. Nov 15 '20 at 21:22
• If $\mathbb{E}\|\beta -\hat{\beta}(k) \|^2$ isn't MSE then what is? Is also says in the caption "shown is the mean-squared error". Nov 15 '20 at 22:22
• I don't see how the y-label $\mathbb{E}\|\beta -\hat{\beta}(k) \|^2$ is bias instead of MSE. Also, eventually, the bias should go to zero when $k=31$, but there will be still error in the estimates. Nov 15 '20 at 23:05
• @Sycorax I would describe bias as $\beta - \mathbb{E}(\hat \beta)$ Nov 15 '20 at 23:13

adding more variables to a linear model doesn't imply better estimates of the true parameters

This is not just estimating variables, but also variable selection. When you only subselect <10 variables, then you are inevitably gonna make an error.

• That is why the error decreases when you are choosing a larger size for the subset. Because more coefficients, which are likely coefficients from the true model, are being estimated (instead of left equal to zero).

• The decrease of the error goes a bit further than $$k=10$$ because of the high correlation between the variables.

The strongest improvement happens before k=10. But with $$k=10$$ you are not there yet, and you are gonna select occasionally the wrong coefficients from the true model.

• Note that after some point, around $$k=16$$, the error goes up when adding more variables.

### Reproduction of the graph

In the R-code at the end I am trying to reproduce the graph for the forward stepwise case. (this is also the question here: Recreating figure from Elements of Statistical Learning)

I can make the figure look similar

But, I needed to make some adjustment to the generation, using $$\beta \sim N(1,0.4)$$ instead of $$\beta \sim N(0,0.4)$$ (and still I do not get the same as the figure which starts at 0.95 and drops down to 0.65, while the MSE computed with the code here is much lower instead). Still, the shape is qualitatively the same.

The error in this graph is not so much due to bias: I wanted to split the mean square error into bias and variance (by computing the coefficient's mean error and variance of the error). However, the bias is very low! This is due to the high correlation between the parameters. When you have a subset with only 1 parameter, then the selected parameter in that subset will compensate for the missing parameters (it can do so because it is highly correlated). The amount that the other parameters are too low will be more or less the amount that the selected parameter will be too high. So on average a parameter will be more or less as much too high as too low.

• The graph above is made with a correlation 0.15 instead of 0.85.
• In addition, I used a fixed $$X$$ and $$\beta$$ (Otherwise the bias would average to zero, more explained about that further).

### Distribution of the error of the parameter estimate

Below you see how the error in the parameter estimate $$\hat\beta_1- \beta_1$$ is distributed as a function of the subset size. This makes it easier to see why the change in the mean square error behaves like it does.

Note the following features

• There is a single peak for small subset sizes. This is because the parameter is often not included in the subset and the estimate $$\hat\beta$$ will be zero making the error $$\hat\beta - \beta$$ equal to $$-\beta$$. This peak decreases in size as the subset size increases and the probability for the parameter to be included increases.
• There is a more or less Gaussian distributed component that increases in size when the single peak decreases in size. This is the error when the parameter is included in the subset. For small subset sizes the error in this component is not centered around zero. The reason is that the parameter needs to compensate for the ommission of the other parameter (to which it is highly correlated). This makes that a computation of the bias is actually very low. It is the variance that is high.

The example above is for fixed $$\beta$$ and $$X$$. If you would change the $$\beta$$ for each simulation then the bias would be every time different. If you then compute the bias as $$\mathbb{E}(\hat \beta - \beta)$$ then you get very close to zero.

library(MASS)

### function to do stepforward regression
stepforward <- function(Y,X, intercept) {
kl <- length(X[1,])  ### number of columns
inset <- c()
outset <- 1:kl

### outer loop increasing subset size
for (k in 1:kl) {
beststep_par <- 0
### inner looping trying all variables that can be added
for (par in outset) {
### create a subset to test
step_set <- c(inset,par)
step_data <- data.frame(Y=Y,X=X[,step_set])
### perform model with subset
if (intercept) {
step_mod <- lm(Y ~ . + 1, data = step_data)
}
else {
step_mod <- lm(Y ~ . + 0, data = step_data)
}
step_RSS <- sum(step_mod$residuals^2) ### compare if it is an improvement if (step_RSS <= beststep_RSS) { beststep_RSS <- step_RSS beststep_par <- par } } bestRSS <- beststep_RSS inset <- c(inset,beststep_par) outset[-which(outset == beststep_par)] } return(inset) } get_error <- function(X = NULL, beta = NULL, intercept = 0) { ### 31 random X variables, standard normal if (is.null(X)) { X <- mvrnorm(300,rep(0,31), M) } ### 10 random beta coefficients 21 zero coefficients if (is.null(beta)) { beta <- c(rnorm(10,1,0.4^0.5),rep(0,21)) } ### Y with added noise Y <- (X %*% beta) + rnorm(300,0,6.25^0.5) ### get step order step_order <- stepforward(Y,X, intercept) ### error computation l <- 10 error <- matrix(rep(0,31*31),31) ### this variable will store error for 31 submodel sizes for (l in 1:31) { ### subdata Z <- X[,step_order[1:l]] sub_data <- data.frame(Y=Y,Z=Z) ### compute model if (intercept) { sub_mod <- lm(Y ~ . + 1, data = sub_data) } else { sub_mod <- lm(Y ~ . + 0, data = sub_data) } ### compute error in coefficients coef <- rep(0,31) if (intercept) { coef[step_order[1:l]] <- sub_mod$$coefficients[-1] } else { coef[step_order[1:l]] <- sub_mod$$coefficients[] } error[l,] <- (coef - beta) } return(error) } ### correlation matrix for X M <- matrix(rep(0.15,31^2),31) for (i in 1:31) { M[i,i] = 1 } ### perform 50 times the model set.seed(1) X <- mvrnorm(300,rep(0,31), M) beta <- c(rnorm(10,1,0.4^0.5),rep(0,21)) nrep <- 500 me <- replicate(nrep,get_error(X,beta, intercept = 1)) ### this line uses fixed X and beta ###me <- replicate(nrep,get_error(X,beta, intercept = 1)) ### this line uses random X and fixed beta ###me <- replicate(nrep,get_error(X,beta, intercept = 1)) ### random X and beta each replicate ### storage for error statistics per coefficient and per k mean_error <- matrix(rep(0,31^2),31) mean_MSE <- matrix(rep(0,31^2),31) mean_var <- matrix(rep(0,31^2),31) ### compute error statistics ### MSE, and bias + variance for each coefficient seperately ### k relates to the subset size ### i refers to the coefficient ### averaging is done over the multiple simulations for (i in 1:31) { mean_error[i,] <- sapply(1:31, FUN = function(k) mean(me[k,i,])) mean_MSE[i,] <- sapply(1:31, FUN = function(k) mean(me[k,i,]^2)) mean_var[i,] <- mean_MSE[i,] - mean_error[i,]^2 } ### plotting curves ### colMeans averages over the multiple coefficients layout(matrix(1)) plot(1:31,colMeans(mean_MSE[1:31,]), ylim = c(0,0.4), xlim = c(1,31), type = "l", lwd = 2, xlab = "Subset size k", ylab = "mean square error of parameters", xaxs = "i", yaxs = "i") points(1:31,colMeans(mean_MSE[1:31,]), pch = 21 , col = 1, bg = 0, cex = 0.7) lines(1:31,colMeans(mean_var[1:31,]), lty = 2) lines(1:31,colMeans(mean_error[1:31,]^2), lty = 3) legend(31,0.4, c("MSE", "variance component", "bias component"), lty = c(1,2,3), lwd = c(2,1,1), pch = c(21,NA,NA), col = 1, pt.bg = 0, xjust = 1, cex = 0.7) ### plotting histogram layout(matrix(1:5,5)) par(mar = c(4,4,2,1)) xpar = 1 for (col in c(1,4,7,10,13)) { hist(me[col,xpar,], breaks = seq(-7,7,0.05), xlim = c(-1,1), ylim = c(0,500), xlab = "", ylab = "", main=paste0("error in parameter ",xpar," for subset size ",col), ) }  • @CagdasOzgenc I will redo my graph and have some ideas to make it more clear (but first I need to rewrite a good function that does the stepwise regression quickly), then I can also see what the effect of setting the X ~ N(1,1) will do. Nov 16 '20 at 7:05 • @CagdasOzgenc I understood that it was about the book. But I guess that it would be nice to create an alternative graph here that allows us to dissect and understand the graph in the book. Nov 16 '20 at 7:50 • I see now where I went wrong with computing bias and variance. I took the errors from all coefficients together but these coefficients will have different distributions of the error and mixing them will reduce the apparent bias (you can even get zero bias when you have two coefficients with exactly opposite bias). Nov 16 '20 at 10:20 • What do you mean? aren't you calculating bias-squared? They cannot cancel each other once they are squared. Also I think X should have 1 mean not$\beta$as it is clear from explanation under the chart that$\beta\$s are 0 mean. Nov 16 '20 at 10:23
• @SextusEmpiricus - I agree, there might be some typing error somewhere in the textbook or they had a very bad seed. As a side note - you do have a small error in your code - you need to sum (not mean) to get the two-norm of beta; averaging only happens across runs. Nov 16 '20 at 17:45

There are good answers here, so I'll try to keep this brief and just add a couple points.

• The point of this figure is to show how close the estimated slopes are to their true values, not how well the model predicts $$y$$ out of sample, or to whether inferences are valid.

adding more variables to a linear model doesn't imply better estimates of the true parameters

• Don't think of this as adding more variables. In all cases, you started with a fixed set of variables determined a-priori. The question is whether you should drop some of those variables for building your final model. Dropping variables based on what you see in your data is generally a bad thing to do. If you retain all variables (assuming you have enough data, which in this case you do) your estimates will be unbiased. Put another way, the variables whose slopes are actually $$0$$ in the data generating process should have slope estimates that are close to $$0$$ in the fitted model. They should be approximately correct. When you drop variables, that's no longer necessarily true.

This case is more complicated, because the variables are all correlated with each other. The correlations mean that the slopes will vary from their true values more widely than they would have if the variables were all mutually orthogonal. As a result, if you pick just the right variables you could reduce the variance somewhat while maintaining the property of unbiasedness. However...

My intuition is that MSE should be lowest around the optimal $$k$$

• That's because your intuition is that stepwise procedures will pick the right variables. Unfortunately, that isn't necessarily what's going to happen. It is very unlikely that you will pick exactly the right variables. And, if you don't pick only the right variables, you will continue to get sampling distributions with higher variance and biased estimates.

Now, let's consider picking the best, say, 15 or 20 variables. What is the probability that we will have included the 10 that we wanted and only thrown away worthless variables that just added noise? It's much better. That's why the curve is lower there.

So a takeaway from this is that if you know how many variables are correct, and you know that they are all included in your dataset, you can focus on retaining some proportion beyond what is needed and you will be likely to only have thrown away garbage. (Of course, I don't find those conditions very realistic, and this discussion only pertains to the slope estimates, not out of sample predictions or statistical inference, so I continue to find stepwise procedures ill-advised.)

What I think is going on is this. This is about variable selection. MSE of estimated betas should be smallest if exactly the correct 10 variables are selected. It should be substantially larger if at least one of these variables is missed. Note that correlation makes this problem worse, because if one of the correct nonzero beta variables is missed, its contribution will be attributed to those that are already in the model because of correlation. This will make their estimators worse, on top of the fact that there is an error from the missing $$\beta$$ itself. It is not true that the effect of correlation is that we can do well with fewer variables than the 10 correct ones, regarding the MSE of the estimators. It may be true for prediction, as the information of a missing variable is compensated for by other correlated variables already in the model. But this is not what the graph is about. The very same effect that may be helpful for prediction will be detrimental for estimation, because the effect of the missing correct nonzero beta variables will be divided among those that are already in the model, affecting their estimation.
This means that the minimum should occur at 10 only if always or almost always exactly the 10 correct variables are selected. But this is very unlikely, because correlation actually makes it very hard to find the correct variables. Chances are that if the procedure selects 11, 12, even 15 variables, still it's not too unlikely that one true nonzero beta variable is missed. True zero beta variables on the other hand will probably have fairly low estimated coefficients anyway, so will not harm the estimator MSE as much as a missed correct nonzero beta variable does. This explains in my view that the estimator MSE goes up only from about $$k=16$$ or even $$k=27$$ or so for stagewise. This seems all fine by me. What it shows is how many variables in this setup need to be selected in order to find all true nonzeroes with large enough probability. 16 seems realistic to me, and it is also clear that stagewise has a hard time in this problem, as it will need many steps to bring initially overestimated parameters down.