---
jupytext:
formats: ipynb,md:myst
text_representation:
extension: .md
format_name: myst
format_version: 0.13
jupytext_version: 1.11.5
kernelspec:
display_name: Python 3 (ipykernel)
language: python
name: python3
---
```{code-cell} ipython3
:tags: [remove-cell]
from ndslib.config import jupyter_startup
jupyter_startup()
```
(ml-selection)=
# Model selection
Machine learning practitioners often emphasize the importance of having very
large training datasets. It's all well and good to say that if we want to make
good predictions, we should collect an enormous amount of data; the trouble is
that sometimes that's not feasible. For example, if we're neuroimaging
researchers acquiring structural brain scans from people with an autism
diagnosis, we can't magically make tens of thousands of scans appear out of thin
air. We may only be able to pool together one or two thousand data points, even
with a monumental multi-site effort like the ABIDE initiative.
If we can't get a lot more data, are there other steps we can take to improve
our predictions? We've already seen how cross-validation methods can indirectly
help us make better predictions by minimizing overfitting and obtaining better
estimates of our models' true generalization abilities. In this section, we'll
explore another general strategy for improving model performance: selecting
better models.
## Bias and variance
We've talked about the tradeoff between overfitting and underfitting, and the
related tradeoff between flexibility and stability in estimation. Now we'll
introduce a third intimately-related pair of concepts: bias and variance.
Understanding the difference between bias and variance, and why we can often
trade one for the other, is central to developing good intuitions about when and
how to select good estimators.
It's probably easiest to understand the concepts of bias and variance through
visual illustration. Here's a classic representation:
![](./figures/bias_variance_dartboard.png)
```{eval-rst}
.. index::
single: Bias variance tradeoff
```
Comparing the top and bottom rows gives us a sense of what we mean by
*variance*: it's the degree of scatter of our observations around their central
tendency. In the context of statistical estimators, we say that an estimator has
high variance if the parameter values it produces tend to vary widely over the
space of datasets to which it could be applied. Conversely, a low-variance estimator
tends to produce similar parameter estimates when applied to different datasets.
Bias, by contrast, refers to a statistical estimator's tendency to produce
parameter estimates that cluster in a specific part of the parameter space. A
high-bias estimator produces parameter estimates that are systematically shifted
away from the correct estimate—even as the size of the training dataset grows
asymptotically. A low-bias estimator, conversely, will converge on the true
parameter estimate as the sample size increases.
### The bias-variance tradeoff
One interesting property of bias and variance is that they tend to trade off.
That is, one can often decrease the variance of an estimator by deliberately
increasing its bias or vice versa. To see this intuitively, imagine the extreme
case of an estimator that always generates an identical prediction for all
observations. Such an estimator would have low variance (actually, *no*
variance) and extremely high bias. Conversely, an estimator that produces
extremely variable estimates must be one capable of freely exploring various
parts of parameter space—and hence, shows a weaker preference for, or bias
towards, certain parts of the space.
This idea has a deep connection to the overfitting/underfitting and
flexibility/stability distinctions we drew earlier. All three distinctions
ultimately reflect the same core idea. A *high-variance* (low-bias) estimator is
a *flexible* one that can identify a wide variety of patterns in the data; this
also makes it more likely to *overfit*. Conversely, a high-bias (low-variance)
estimator is a *stable* one that tends to produce similar parameter estimates no
matter what data it's presented with; such an estimator will frequently
*underfit* because it lacks the flexibility needed to capture complex patterns
in the data.
Neuroimaging researchers coming from a traditional scientific background often
think of bias as a dirty word -— i.e., as something to be minimized at all
costs. From a machine learning standpoint, however, there's nothing
intrinsically bad about bias. The bias of a particular model is nothing more
than its systematic tendency to produce parameter estimates that cluster within
a particular part of the space. Not only is this not a bad thing, but it's often
essential to producing good predictions. The bias of a given estimator -- its
tendency to give certain kinds of solutions -- is often referred to as its
*inductive bias*. The key question is then: *to what degree does the inductive
bias of a given estimator match the processes that give rise to our data?*
```{eval-rst}
.. index::
single: Inductive bias
```
In cases where there's a match between an estimator and the data it's applied
to, a biased estimator will tend to perform better than a less biased one.
Conversely, if an estimator has the *wrong* bias, it will tend to perform more
poorly than an estimator with less bias.
This is a fairly abstract idea, so let's make it concrete. Suppose you park your
car at the airport parking lot when you head off for a weekend vacation.
Consider two possible strategies you could use to find your car on your return.
First, you might rely on your memory of where you left the car; perhaps you
think you left it somewhere near that light pole over there in the `B` area of
the lot. Second, you could just randomly roam around the parking lot looking for
a car that looks like yours.
Which of these is the better car-finding strategy? Well, it depends! If you have
a good memory, the high-bias strategy -— only searching within a particular part of
the lot —- will lead you to your car more quickly than a random walk around the
entire lot. But if your memory is poor, the biased strategy is likely to be
counterproductive, because your car won't be where you remember leaving it, no
matter how much time you spend walking around that specific area. In the latter
case, you would have been better off randomly walking around until you found
your car.
## Regularization
```{eval-rst}
.. index::
single: Regularization
```
Our car parking example is a bit contrived, but it illustrates a core principle at the heart of good model selection: *if we have useful prior knowledge about the domain, we should pick an estimator that builds that knowledge into its assumptions*. For example, if we have a set of 1,000 potential features we could use in our model, but have good reason to believe that only a small fraction of those features make a meaningful contribution to the outcome variable, it seems like a good idea to use an estimator that can somehow *select* a subset of features. Using an estimator that acts as if all of our features are potentially important (as linear regression does) is likely to lead to overfitting.
The principle of introducing background knowledge into our model to improve
predictions and prevent overfitting is known in machine learning and statistics
as *regularization*. When we say we're using a regularized estimator, what we're
typically saying is that we've taken an existing estimator (e.g., an OLS
estimator) and constrained its behavior with additional information, so that
it's *biased* towards certain parts of the parameter space.
(penalized_regression)=
### Penalized regression
To researchers who learned statistics in biomedical or social science training
programs, it may be tempting to think of least-squares estimation as a kind of
automatic default, with all other procedures requiring some special
justification to use. But there's nothing special about OLS. Well, it's
*kind of* special in that it always minimizes the squared error in our training
sample. But as we've already seen, we usually care about minimizing test error,
not training error (plus, sometimes we don't want to quantify error in terms of
a least-squares criterion). So if we can get our hands on an estimator that's
biased relative to OLS, but tends to reduce error in the *test* dataset, we
probably still want to use it -— even if it doesn't reduce error in the *training*
data quite as well as OLS.
It turns out that there are many such estimators. One of the most widely used
class of regularized estimators for regression problems is *penalized regression*. These methods generally start with garden variety OLS and then add
a little twist in the form of a *penalty* parameter.
Consider the cost function we seek to minimize in OLS:
$Cost = RSS = \sum^N_i(y_i - \sum^P_j\beta_jx_{ij})^2$
Here, $RSS$ is the residual sum of squares, $N$ is the number of samples, and
$P$ is the number of features. Our goal is to obtain the set of $\beta$
coefficients that minimize the $RSS$.
Now consider a slight variation:
$Cost = \sum^N_i(y_i - \sum^P_j\beta_jx_{ij})^2 + \alpha\sum^P_j|\beta_j|$
The only difference here is that, in addition to the $RSS$, our cost function
includes a second term, $\alpha\sum^P_j|\beta_j|$ -— the sum of absolute values
of the $\beta$ coefficients, weighted by the *penalty* parameter, $\alpha$. The
addition of the penalty parameter is the reason we refer to this variation on
linear regression as penalized regression. The addition of the penalty may seem
like a small change, but it has major consequences, which we'll explore shortly.
The above form of penalized regression is known as *lasso regression*, and it's
very widely used. There's also another very common form of penalized regression
known as *ridge regression*. It looks a lot like lasso regression, except
instead of computing the penalty as the sum of absolute coefficients, we take
the sum of their squares:
$Cost = \sum^N_i(y_i - \sum^P_j\beta_jx_{ij})^2 + \alpha\sum^P_j\beta_j^2$
The difference between lasso and ridge regression may seem even smaller, but it
again has important implications. Let's explore the behavior of each of these
methods.
```{eval-rst}
.. index::
single: Regularization
```
```{eval-rst}
.. index::
single: Penalized regression
```
```{eval-rst}
.. index::
single: Ridge regression
```
```{eval-rst}
.. index::
single: Lasso
```
#### Lasso regression
```{eval-rst}
.. index::
single: Lasso regression
```
Let's take the lasso first. Let's think about what the addition of the penalty
term to the standard RSS criterion does, conceptually. By making the overall
cost depend on a (weighted) sum of absolute coefficients, we're saying that
there's a certain cost to having large coefficients in our model. If the net
negative contribution of a given feature to the RSS is smaller than its net
positive contribution to the penalty, then it will be beneficial to "shrink"
that feature's coefficient to the point where the two terms perfectly cancel
each other out. Thus, penalized regression *biases* the estimated regression
coefficients by shrinking at least some of them.
For reasons outside the scope of our discussion here, the lasso doesn't just
shrink coefficients; it shrinks them to zero. That is, as we increase the
penalty, more features drop out of our model, leaving us with a simpler
prediction equation involving fewer features. For this reason, the lasso is
often described as implicitly including a *feature selection* step.
```{eval-rst}
.. index::
single: Feature selection
```
To see how this works in practice, let's return to our ABIDE-II data and
visualize what happens to our regression coefficients estimated with lasso
regression as we increase the penalty parameter. Here, we'll plot the
*coefficient paths*: how each coefficient changes as $\alpha$ changes
(judiciously picking a reasonable range of $\alpha$ values). We've implemented
a `plot_coef_path` function that does this. Note also that we standardize our
features to mean-zero and unit-variance because lasso and ridge are sensitive
to scale.
```{code-cell} ipython3
from ndslib.viz import plot_coef_path
from ndslib.data import load_data
abide_data = load_data("abide2")
features = abide_data.filter(like='fs')
phenotypes = abide_data.iloc[:, :6]
from sklearn.linear_model import Lasso
from sklearn.preprocessing import scale
number_of_features = 200
X = scale(features.sample(number_of_features, axis=1, random_state=99))
y = phenotypes['age']
import numpy as np
alpha = np.logspace(-3, 1, 100)
ax = plot_coef_path(Lasso, X, y, alpha, max_iter=2000)
ax.set_title("Coefficient paths for lasso regression")
ax.set_xlim(1e-3, 10);
```
In the figure produced by the `plot_coef_path` helper function, each line
represents a different feature. We're using 200 randomly selected features (from
the full set of 1,440) as predictors, so there are 200 lines. The x-axis
displays the penalty parameter used for lasso estimation; the y-axis displays
the resulting value of the coefficients. Notice how, at a certain point (around
`alpha = 0.01`), coefficients start to disappear. They're not just small;
they're *gone* (i.e., they've been "shrunk" to 0). This built-in feature
selection is a very useful property of the lasso. It allows us to achieve an
arbitrarily sparse solution simply by increasing the penalty.
(ml-selection_ex1)=
#### Exercises
1. Change the `number_of_features` variable to be much larger (e.g., 400) or
much smaller (e.g., 100). How are the coefficient paths affected by this
change? Why do you think these changes happen?
2. What happens if you use the site-residualized ages instead of the age variable?
Of course, as you know by now, *there's no free lunch in machine learning*.
There must be some price we pay for producing more interpretable solutions with
fewer coefficients, right? There is. We'll look at it in a moment; but first,
let's talk about ridge regression.
#### Ridge regression
```{eval-rst}
.. index::
single: Ridge regression
```
Ridge regression, as noted above, is mathematically very similar to lasso regression. But it behaves very differently. Whereas the lasso produces sparse solutions, ridge regression, like OLS, always produces dense ones. Ridge *does* still shrink the coefficients (i.e., their absolute values get smaller as the penalty increases). It also pushes them towards a normal distribution. The latter property is the reason a well-tuned ridge regression model usually outperforms OLS: you can think of ridge regularization as basically a way of saying *it may **look** like a few of the OLS coefficients are way bigger than others, but that's probably a sign we're fitting noise. In the real world, outcomes usually reflect the contributions of a lot of different factors. So let's squash all of the extreme values towards zero a bit so that all of our coefficients are relatively small and bunched together.*
```{code-cell} ipython3
from sklearn.linear_model import Ridge
# Coefficient paths for ridge regression, predicting age from 30 features
alpha = np.logspace(-5, 5, 100)
ax = plot_coef_path(Ridge, X, y, alpha)
ax.set_title("Coefficient paths for ridge regression")
ax.set_xlim(1e-5, 1);
```
It may look to the eye like the coefficients are shrinking all the way to zero
as the penalty increases, but we promise you they aren't. They may get very
small, but unlike lasso regression, ridge coefficients will never shrink to
nothing.
#### Effects of regularization on predictive performance
We've looked at two different flavors of penalized regression, and established
that they behave quite differently. So, which one should we use? Which is the
*better* regularization approach? Well, as with everything else *it depends*. If
the process(es) that generated our data are sparse (or at least, have a good
sparse approximation) -— meaning, a few of our features make very large
contributions, and the rest don't matter much —- then lasso regression will tend to
perform better. If the data-generating process is dense—i.e., lots of factors
make small contributions to the outcomes—then ridge regression will tend to
perform better.
In practice, we rarely know the ground truth when working with empirical data.
So we're forced to rely on a mix of experience, intuition, and validation
analyses.
Let's do some of the latter here. To try to figure out what the optimal lasso
and ridge penalties are for our particular problem, we'll make use of Scikit
Learn's `validation_curve` utility, which allows us to quickly generate training
and testing performance estimates for an estimator as a function of one of its
parameters. As alluded to in {numref}`ml-validation`, the `validation_curve` is
very similar to the `learning_curve` we've seen before, except that instead of
systematically varying the dataset size, we systematically vary one of the
estimator's parameters (in this case, `alpha`). Let's generate a validation
curve for lasso regression, still using 200 random features to predict age.
We'll add OLS regression $R^2$ for comparison and use our `plot_train_test`
helper function to visualize it.
```{code-cell} ipython3
from sklearn.model_selection import validation_curve, cross_val_score
from sklearn.linear_model import LinearRegression
from ndslib.viz import plot_train_test
x_range = np.logspace(-3, 1, 30)
train_scores, test_scores = validation_curve(
Lasso(max_iter=5000), X, y, param_name='alpha',
param_range=x_range, cv=5, scoring='r2')
ols_r2 = cross_val_score(
LinearRegression(), X, y, scoring='r2', cv=5).mean()
ax = plot_train_test(
x_range, train_scores, test_scores, 'Lasso',
hlines={'OLS (test)': ols_r2})
ax.set_title("Age prediction performance")
ax.set_xlim(1e-3, 10)
```
Notice that, in the training dataset, $R^2$ decreases monotonically as alpha
increases. This is necessarily true (Why do you think that is? *Hint*: think of
how $R^2$ is defined, and its relationship to the least-squares component of the
lasso cost function.)
In contrast, the error in the test dataset is lowest within a specific range of
`alpha` values. This optimum will vary across datasets, so you'll almost always
need to tune the alpha parameter (e.g., using a validation curve like the one
above) to get the best performance out of your model. This extra degree of
complexity in the model selection process is one of the reasons you might *not*
want to always opt for a penalized version of OLS, even if you think it can
improve your predictions. If you muck things up and pick the wrong penalty, your
predictions could be way worse than OLS!
While we've only looked at one dataset, these are fairly typical results; in
many application domains, penalized regression models tend to outperform OLS
when using small to moderate-sized samples. When samples get very large relative
to the number of features, the performance difference usually disappears.
Outside of domains with extremely structured and high-dimensional data, this
latter principle tends to be true more globally -— i.e., "more data usually beats
better algorithms."
(ml-selection_ex2)=
#### Exercise
Use a similar approach to compare ridge regression with OLS.
## Beyond linear regression
We've focused our attention here on lasso and ridge regression because they're
regularized versions of the OLS estimator most social scientists are already
familiar with. But just as there's nothing special about OLS, there's also
nothing special about (penalized) regression in general. There are many
different estimators we could use to generate predictions. Each has its own
inductive biases and will perform better on some datasets and worse on others.
The science (or perhaps art) of machine learning lies in understanding the
available estimators, data, and validation methods well enough to know how to
tailor optimal (or at least good) machine learning workflows to the specific
problem at hand. Let's forge ahead with one more widely-used class of algorithms
that differs quite a bit from the regression algorithms you've seen here so far.
### Random forests
```{eval-rst}
.. index::
single: Decision trees
```
```{eval-rst}
.. index::
single: Random forests
```
To illustrate how easy it is to try out completely different types of estimators
in Scikit-learn -- and also draw attention to the fact that estimators based on
linear regression represent only a small part of the space--let's repeat some of
our earlier regression analyses using a very different estimator: the
`RandomForestRegressor`. Random forests are essentially collections (or
*ensembles*) of decision trees. A decision tree is a structure for generating
classification (classification trees) or regression (regression trees)
predictions. Here's an example of what a regression tree might look like if used
in our ABIDE-II dataset to try and predict age from structural brain features:
![decision tree example](figures/tree.png)
At each node, we evaluate a particular conditional. For example, the first
question we ask is: is the subject's observed value on the fsCT_218 variable
less than or equal to 2.49, or is it greater? We follow the appropriate branch
to the next node and answer the question we find there. We repeat this process
until we reach a terminal node (or "leaf"), whereupon the tree produces the
predicted value for that observation.
Random forests extend this idea by "bagging" multiple decision trees and
aggregating their outputs. Decision trees are very flexible and have the
propensity to overfit. By averaging over a lot of different trees (each one
generated by, e.g., resampling the data), we hopefully stabilize our predictions
and reduce overfitting. Random forests are popular because they're extremely
powerful, their constituent trees are highly interpretable, and they tend to
perform well when we have a lot of data. The downsides are that being very
flexible, they have a tendency to overfit, and careful tuning may be required to
achieve good performance. They can also be quite slow, as performance tends to
improve with the number of decision trees. In the example below we will reduce
tree complexity and mitigate overfitting by performing an explicit feature
selection step before fitting the model. This is done using the `SelectKBest`
object and an `f_regression` criterion, which selects the $k$ features that have
the strongest correlation with the target (based on an F-test -- hence the
name). Ideally, we'd want to use a larger number of trees, but this would be very
slow, so we stick with a small number (10) for demonstration purposes (but
please do experiment with changing these numbers!)
```{code-cell} ipython3
from sklearn.ensemble import RandomForestRegressor
from sklearn.feature_selection import SelectKBest, f_regression
from ndslib.viz import plot_learning_curves
number_of_features = 50
alpha = 10
number_of_trees_in_forest = 10
selector = SelectKBest(f_regression, k=number_of_features)
X = selector.fit_transform(features, y)
estimators = [Ridge(alpha), RandomForestRegressor(number_of_trees_in_forest)]
labels = ["Ridge regression", "Random forests"]
plot_learning_curves(estimators, [X, X], y, [100, 200, 500, 900], labels)
```
In this case, the random forest estimator slightly outperforms ridge regression,
particularly as the sample size gets larger. But it's worth reiterating our
earlier point: this isn't some kind of law, or even a particularly good
heuristic; if you play with various parameters above (e.g.,
`number_of_features`, `alpha`, etc.), you can easily find regimes under which
ridge regression performs better.
#### Interpreting random forests
Unlike regression-based methods, decision trees (and hence random forests) don't
have linear coefficients. Instead, we can obtain *feature importances*, which
tell us how important each feature is to the overall prediction (more important
features are better at reducing uncertainty, and occur closer to the root of the
tree). Let's randomly sample 30 features and work with those. After fitting, we
will print out the names of the 10 most important features.
```{code-cell} ipython3
rf = RandomForestRegressor(100)
X = features.sample(30, axis=1)
rf.fit(X, y)
import pandas as pd
print(pd.Series(
rf.feature_importances_, index=X.columns).sort_values(ascending=False).head(10))
```
Be aware that there are some important subtleties related to the
interpretation of feature importances, not least of which is that, as with any
model, feature importances are configural: a feature may seem important when
modeled alongside one set of other features, but unimportant when included with
a different set of features.
#### Visualizing trees
We can also plot the individual decision trees. Note that different trees can
often exhibit very different structures despite performing equally well
predictively—an observation that should lead us to exercise caution when
interpreting decision trees, even if they *seem* straightforward. In this
particular case, as we saw above, a handful of features dominate the list of
feature importances, so we might expect at least the first few nodes to look
relatively similar. Let's take a look at the first two trees in the forest.
```{code-cell} ipython3
from ndslib.viz import plot_graphviz_tree
plot_graphviz_tree(rf.estimators_[0], X.columns)
```
```{code-cell} ipython3
plot_graphviz_tree(rf.estimators_[2], X.columns)
```
### Additional resources
To learn more about the properties of the lasso and ridge algorithms, we'd recommend reading an excellent tutorial by Ryan Tibshirani that Larry Wasserman amended and posted on [his website](https://www.stat.cmu.edu/~larry/=sml/sparsity.pdf). If you want to think a bit more about the implications of the $\alpha$ parameter, particularly in ridge regression, you can read a paper that one of us wrote on the topic {cite}`Rokem2020-hz`.