# PC lab 6: Linear model selection - Regularization

We have already seen linear regression to tackle regression problems. With linear regression, we model a continous outcome as a linear function of the features:

\begin{equation}
\hat{y} = w_{0}x_{0} + w_{1}x_{1} + ... + w_{p}x_{p} = \sum\limits_{i=0}^{p}w_{i}x_{i}
\end{equation}

This works well when there are a lot of training observations and when the number of features (the dimensionality of the problem) is not too large. However, there are a couple of situations where ordinary linear regression might give problems:

* When the number of features $p$ becomes large with respect to the number of observations $n$, the variance of the model weights estimated by linear regression increases, which might result in poor predictive performance. Futhermore, there is no longer an analytical solution provided by least squares when $p > n$. 
* It is possible that there are a lot of uninteresting or redundant features. If we want a sparse and interpretable model, we might want to do feature selection to reduce $p$. 

In this lab, we will cover two solutions to the problems above: subset selection and regularization methods.

## Subset selection methods

In subset selection, we only use a subset of the features that are available. The goal is to come up with a model that is sparse and that generalizes better to unseen data. There are two main strategies for subset selection: in *best subset selection*, we fit all the $p \choose k$ models for each $k = 1, 2.. p$ and retain the best model for each $k$. Finally, we select the model that performs best on some measure that controls for overfitting: 



This becomes quickly unfeasible for large values of $p$. Therefore, an alternative approach is to perform *stepwise selection*, which explores a much smaller set of feature combinations. Stepwise selection can be performed either backward or forward. For large $p$ they are the only computationally feasible subset selection methods. 



Finally, it is important to account for the fact that the MSE (or, equivalently, the RSS) will always go down on the training data as we add more and more features. To select the best model out of several candidates, it is important to have an estimate of the test error of each model. This can be done indirectly by using a metric that penalizes for model complexity such as the AIC or the adjusted $R^2$. Another option is to use cross-validation to get an estimate of the test error.

**Let's apply subset selection on two datasets.** The first dataset contains features of mixtures used to produce concrete. The goal is to predict the compressive strength of the concrete.


```python
import urllib.request
urllib.request.urlretrieve("https://raw.githubusercontent.com/gdewael/teaching/main/predmod/lab6/concreteComprStrength.txt", "concreteComprStrength.txt")
urllib.request.urlretrieve("https://raw.githubusercontent.com/gdewael/teaching/main/predmod/lab6/fc_data_new.csv", "fc_data_new.csv")
urllib.request.urlretrieve("https://raw.githubusercontent.com/gdewael/teaching/main/predmod/lab6/communities.csv", "communities.csv")
urllib.request.urlretrieve("https://raw.githubusercontent.com/gdewael/teaching/main/predmod/lab6/meatNIR1000.csv", "meatNIR1000.csv")
urllib.request.urlretrieve("https://raw.githubusercontent.com/gdewael/teaching/main/predmod/lab6/meatNIR1000.colnames", "meatNIR1000.colnames")
```


```python
# Read in the data
import pandas as pd

concretedata = pd.read_table('concreteComprStrength.txt', delim_whitespace=True, header=0, index_col=None)
print(concretedata.shape)
concretedata.head()
```

Let's build up an intuition for the code we need for doing best subset selection with linear regression. We will implement algorithm 6.1 from the book as shown above. First we split the data and look at how we can generate all possible combinations of features:


```python
X = concretedata.drop(['compressiveStrength'], axis=1)
y = concretedata['compressiveStrength']

from itertools import chain, combinations

list_of_things = [1, 2, 3]

for subset in chain.from_iterable(
    combinations(list_of_things, r) 
    for r in range(1, len(list_of_things)+1)
    ):

    print(subset)
```

Now we can use this code to iterate not through `list_of_things` but through the possible features. For every subset we fit a model and keep it in a dictionary that records the score for that specific subset of features:


```python
from sklearn.linear_model import LinearRegression

X = concretedata.drop(['compressiveStrength'], axis=1)
y = concretedata['compressiveStrength']

from itertools import chain, combinations

list_of_things = X.columns

scoring_dict = {}

for subset in chain.from_iterable(
    combinations(list_of_things, r) 
    for r in range(1, len(list_of_things)+1)
    ):

    X_subset = X[list(subset)]

    model = LinearRegression()

    model.fit(X_subset, y)

    scoring_dict[subset] = model.score(X_subset, y)
```

Taking five random samples from this dict to get an idea what's in it:


```python
import random
random.sample(list(scoring_dict.items()), 5)
```

Getting the best model:


```python
max(scoring_dict, key=scoring_dict.get), scoring_dict[max(scoring_dict, key=scoring_dict.get)]
```

We see that the best model (R² score-wise) is the one where all features are present. However, this is estimated on training data, so it might very well be that this is overfit. To do this in a better way, we will use cross-validation.
When doing this, it is important that for every model we train, we use exactly the same splits, otherwise, extra unwanted sources of variation are present in our model performance estimates. The following code block expands upon the previous code block with these things in mind: **This code may take a little while to run**


```python
from sklearn.model_selection import KFold, cross_validate
import numpy as np


X = concretedata.drop(['compressiveStrength'], axis=1)
y = concretedata['compressiveStrength']
list_of_things = X.columns
scoring_dict = {}

#NEW:
# Define splitter object outside loop to ensure same splits over loop
splitter = KFold(n_splits=5, shuffle=True, random_state=None)
####

for subset in chain.from_iterable(
    combinations(list_of_things, r) 
    for r in range(1, len(list_of_things)+1)
    ):

    X_subset = X[list(subset)]

    model = LinearRegression()

    # NEW: instead of simply fitting the model on all train data,
    # now we use cross validation to fit and evaluate model
    # cross_validate() returns a dict with key 'test_score' containing
    # a list of  CV performances of every split
    scores = cross_validate(model, X_subset.values, y.values, cv=splitter)
    scoring_dict[subset] = np.mean(scores['test_score'])
```

#### Exercise 1: how much models did we just train, exactly?


```python

```

Let's visualize how the test performance changes in function of how many features we had in the model:


```python
import matplotlib.pyplot as plt

plt.scatter([len(k) for k in scoring_dict.keys()], list(scoring_dict.values()), alpha = 0.25)
plt.xlabel('number of features')
plt.ylabel('R² score')
# taking the best model for every 'n' number of features and plotting them separately:
best_score_n = [
    max([v for k, v in scoring_dict.items() if len(k) == i])
    for i in range(1, X.shape[1]+1)
]
plt.scatter(range(1, X.shape[1]+1), best_score_n)
print('Best feature combination:', max(scoring_dict, key=scoring_dict.get))
```

We see that one certain combination of 5 features already results in a performance on +- the same level as the one with all features. This finding begs the question: could we have found this combination more efficiently, without testing out all possible feature combinations?

Now, we will introduce a second dataset. This time, the features are measurements from a flow cytometry experiment. The 'SC' features measure scatter, and say something about the morphologhy of the cells (FSC: forwad scatter, SSC: sideway scatter). The 'FL' features are fluorescence features from different parts of the spectrum. There are two possible bacterial species present in the dataset. The goal is to classify the correct species based on the measurements from the flow cytometer. Species number one corresponds to *Pseudomonas putida*, while species number 6 is *Brachybacterium faecium*. We will use logistic regression to do the classification.

This dataset has 10 features instead of the 8 features encountered in the previous dataset. This makes the number of possible combinations (power set): $(2^{10} -1 ) \cdot n_{folds} = 5115$, instead of $(2^{8} -1) \cdot n_{folds} = 1275$. We could repeat the same *best subset selection* method, but we will be smarter about it and implement *stepwise forward selection*.


```python
# Read in the data
bacterialdata = pd.read_csv('fc_data_new.csv', index_col=0)
bacterialdata = bacterialdata.drop(['Width', 'Time'], axis=1).reset_index(drop=True)
```


```python
bacterialdata.head()
```

#### Exercise: Considering we have 10 features and will use 5-fold cross validation: how much models will we train using stepwise forward selection?


```python

```

#### Exercise: Complete the code below to implement stepwise forward selection:


```python
from sklearn.linear_model import LogisticRegression

X = bacterialdata.drop(['species'], axis=1)
y = bacterialdata['species']
list_of_things = X.columns
scoring_dict = {}

splitter = KFold(n_splits=5, shuffle=True, random_state=None)


features_already_in_model = set()
features_not_in_model = set(X.columns)

# iterate through the numbers of features that we can add at every step
for k in np.arange(1, X.shape[1] + 1):
    # at every step: iterate through the features that we can still add
    # keep a separate scoring dictionary for features added at this step:
    scoring_dict_k = {} 
    for feature in features_not_in_model:
        # select the previously added features + a feature to (maybe) add
        X_subset = X[list(features_already_in_model) + [feature]] 

        ... # call the model

        scores = ... # do cross validation on the model

        scoring_dict_k[feature] = ... # record the mean cross validation score

    feature_to_add = max(scoring_dict_k, key=scoring_dict_k.get)
    features_already_in_model.add(feature_to_add)
    features_not_in_model.remove(feature_to_add)
    scoring_dict[tuple(features_already_in_model)] = max(scoring_dict_k.values())

import matplotlib.pyplot as plt

plt.scatter([len(k) for k in scoring_dict.keys()], list(scoring_dict.values()))
plt.xlabel('number of features')
plt.ylabel('accuracy')
```


```python
for k, v in scoring_dict.items():
    print(k, v)
```

#### Exercise: Let's say we performed this feature selection method with k-nearest neighbor classifiers instead of logistic regression, and let's say we wanted to additionally tune the number of nearest neighbors we consider. How would you go about combining both feature selection and hyperparameter tuning?

## Regularization methods: ridge regression and the lasso


Regularization is one of the most important concepts in machine learning to avoid overfitting. It comes in many forms. In linear regression, regularization techniques typically constrain the coefficient estimates. In return for a little extra bias, this reduces the variance of the coefficient esimates. The two main shrinkage techniques are **ridge regression** and the **lasso**.

**Ridge regression penalizes the sum of the squares of the model weights by adding a term to the MSE loss function**:


The lasso does a similar thing, but penalizes the absolute value of the model coefficients. The effect of both approaches is that the model coefficients are shrunk towards zero, resulting in less overfitting and less variance in the predictions (at the cost of a little more bias). We will apply both models on two datasets.

Note that we shown the equation for regularization pentalties here for linear regression with MSE loss, but the technique is equally applicable to classification loss functions (e.g. cross entropy)

### Linear models for high dimensional data

We will apply ridge regression on the Communities and Crime Data Set. The dataset contains 123 population statistics on 1994 communities. We would like to predict the number of violent crimes per 100000 inhabitants. This is the final column of the dataframe. Of the 123 features, a lot contain missing values, so we will drop these columns and use only 99 features.


```python
import pandas as pd
import numpy as np

data = pd.read_csv('./communities.csv')

# Drop columns with missing values
data = data.dropna(axis=1)
print(data.shape)
data.head()
```

<div class="alert alert-success">
    <h2>Exercise</h2>
    <p>Use linear regression and ridge regression to predict the number of violent crimes per 100,000 inhabitants. Use 5-fold cross-validation to evaluate both models. The scikit-learn implementation of RidgeCV automatically performs cross-validation to tune the hyperparameter that determines the amount of regularization, so you don't need to implement a second cross-validation loop. Which model performs best?</p>
</div>

**Before we apply ridge regression, it's important to make sure that all the features are on the same scale. If one of the features is on a completely different scale (let's say, income can be measured in dollars or in thousands of dollars), this might lead the ridge regression coefficient to change substantially because of the penalty term in the optimization problem. We can make sure that all the features are on the same scale by use of the standard scaling: (see book p. 232)**

\begin{equation}
\tilde{x}_{ij} = \frac{x_{ij}-\bar{x}_{j}}{\sqrt{\frac{1}{n}\sum_{i=1}^{n}(x_{ij} - \bar{x}_{j})^2}}
\end{equation}

We can do this with the ```StandardScaler``` from scikit-learn. We will do the scaling each time in the cross-validation loop using only training data statistics (*think for yourselves: do you still remember why we only use training statistics to do this?*). 


```python
from sklearn.linear_model import LinearRegression, RidgeCV
from sklearn.model_selection import KFold
from sklearn.preprocessing import StandardScaler
```


```python
# Select X and y
y = data['ViolentCrimesPerPop'].values
X = data.drop(['ViolentCrimesPerPop'], axis=1).values

kf = KFold(n_splits=5)

linreg_scores = []
ridge_scores = []

for train_index, test_index in kf.split(X):
    # select training and testing X & y from indices
    
    # scale X

    # call linear regression

    # fit linear regression

    # evaluate and append score to linreg_scores

    # do the same for ridgeCV

print('Average validation fold R² of linear regression: {}'.format(np.mean(linreg_scores)))
print('Average validation fold R² of ridge regression: {}'.format(np.mean(ridge_scores)))
```

Now suppose that we don't have 99 features, but four times as many features. And suppose that a lot of features are correlated. This situation is very common in lots of datasets. We will mimic this situation by adding correlated features to our original feature matrix:


```python
X_1 = X + np.random.normal(loc=0, scale=0.05, size=(X.shape))
X_2 = X + np.random.normal(loc=0, scale=0.1, size=(X.shape))
X_3 = X + np.random.normal(loc=0, scale=0.01, size=(X.shape))
X_expanded = np.concatenate((X, X_1, X_2, X_3), axis=1)

X_expanded.shape
```

<div class="alert alert-success">
    <h2>Exercise</h2>
    <p>Repeat the comparison between ridge regression and linear regression from above, but with the new feature matrix that contains correlated features. What happens with the performance of linear regression? What happens with the regularization parameter? Now, add even more correlated features and repeat the analysis.</p>
</div>


```python
kf = KFold(n_splits=5)

linreg_scores = []
ridge_scores = []

for train_index, test_index in kf.split(X_expanded):
    # select training and testing X & y from indices
    
    # scale X

    # call linear regression

    # fit linear regression

    # evaluate and append score to linreg_scores

    # do the same for ridgeCV
    

print('Average validation fold R² of linear regression: {}'.format(np.mean(linreg_scores)))
print('Average validation fold R² of ridge regression: {}'.format(np.mean(ridge_scores)))
```


Finally, it's interesting to look at how the amount of regularization in ridge regression and lasso regression affects the magnitudes of the fitted weights.


```python
from sklearn.linear_model import Ridge, Lasso
from matplotlib import pyplot as plt
%matplotlib inline

y = data['ViolentCrimesPerPop'].values
X = data.drop(['ViolentCrimesPerPop'], axis=1).values
 
# Scale X, note the "wrong" way in which we do it here. This code is only to make a nice plot!
scaler = StandardScaler()
x = scaler.fit_transform(X)
x = X
alphas = np.logspace(-5,2,100)

ridge_coefficients = np.ndarray(shape=(50, len(alphas)))
lasso_coefficients = np.ndarray(shape=(50, len(alphas)))

for i,a in enumerate(alphas):
    ridgemodel = Ridge(alpha=a)
    lassomodel = Lasso(alpha=a, max_iter=10000)
    
    ridgemodel.fit(X,y)
    lassomodel.fit(X,y)
    
    ridge_coefficients[:, i] = ridgemodel.coef_[:50]
    lasso_coefficients[:, i] = lassomodel.coef_[:50]
```


```python
fig, ((ax1, ax2)) = plt.subplots(figsize=(15,10), nrows=2, sharex=True)

for c in range(ridge_coefficients.shape[0]):
    pd.Series(ridge_coefficients[c,:]).plot(ax=ax1)
    pd.Series(lasso_coefficients[c,:]).plot(ax=ax2)

ax2.set_xlabel('Regularization parameter').set_fontsize(20)
ax1.set_ylabel('Model weights')
ax2.set_ylabel('Model weights')
ax1.set_title('Model weights for increasing regularization - Ridge regression').set_fontsize(20)
ax2.set_title('Model weights for increasing regularization - Lasso').set_fontsize(20)
ax2.get_xaxis().set_ticks(alphas);
```

<div class="alert alert-success">
    <h2>Exercise</h2>
    <p>Make sure you understand the plots above. What is the main difference between ridge regression and the lasso?
    </p>
</div>

### Feature selection with the lasso

Suppose that we are interested in selecting only a couple of features out of a high dimensional dataset. Let's fit ridge regression and lasso regression on the data and look at the model coefficients for both models:


```python
import matplotlib.pyplot as plt
%matplotlib inline
from sklearn.linear_model import LassoCV
from sklearn.model_selection import train_test_split

ridgemodel = RidgeCV(cv=5)
lassomodel = LassoCV(cv=5, max_iter=10000)

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)

# Scale X
scaler = StandardScaler()
scaler.fit(X_train)
X_train = scaler.transform(X_train)
X_test = scaler.transform(X_test)

ridgemodel.fit(X_train, y_train)
ridgescore = np.round(ridgemodel.score(X_test, y_test),2)
lassomodel.fit(X_train, y_train)
lassoscore = np.round(lassomodel.score(X_test, y_test),2)

# Plot of the coefficients for ridge regression
fig, ((ax1, ax2)) = plt.subplots(figsize=(15, 10), nrows=2)
pd.Series(ridgemodel.coef_).plot(kind='bar', ax=ax1)
ax1.set_title('Ridge regression coefficients. Test set R²: {}'.format(ridgescore)).set_fontsize(40)
ax1.get_xaxis().set_ticks([])
ax1.set_xlabel('Features')

# Plot of the coefficients for the lasso
pd.Series(lassomodel.coef_).plot(kind='bar', ax=ax2)
ax2.set_title('Lasso regression coefficients. Test set R²: {}'.format(lassoscore)).set_fontsize(40)
ax2.get_xaxis().set_ticks([])
ax2.set_xlabel('Features');
```

The lasso applies regularization by constraining the sum of the absolute values of the model coefficients (the L1-norm). A result of this is that a lot of model coefficients are set to zero: the lasso performs **feature selection**. This is not the case for ridge regression: the model weights are rarely set to zero. Feature selection is a nice property if we want an interpretable model. Let's list the features with non-zero coefficients in the lasso: 


```python
import seaborn as sns
sns.set_style('whitegrid')

fig, ax = plt.subplots(figsize=(12,15))
sns.barplot(x= pd.Series(lassomodel.coef_[lassomodel.coef_ != 0]),
            y=data.drop(['ViolentCrimesPerPop'], axis=1).columns[lassomodel.coef_ != 0],
            palette=['steelblue' if n > 0 else 'coral' for n in lassomodel.coef_[lassomodel.coef_ != 0]],
            ax=ax);
ax.set_title('Lasso model coefficients for the Community Crime dataset').set_fontsize(20);
```

#### Extra note: **machine learning bias**

Take a look at which features our model has identified as important for predicting whether a community will have a higher percentage of crime. You will perhaps be shocked (or unsurprised) to find that our model has quite a racial bias.

An in-depth discussion of this topic is out-of-scope for this PC-lab. However we note two important factors contributing to this bias here:
- Sampling bias: were communities selected uniformly, or more pertinently "fairly" (How do you define "fairness"?)
- Correlation between features: socio-economical factors causal for crime rates will be correlated with all of the features. Remember that our linear model simply captures patterns/correlations between features and outcomes. Without a view of the causal processes underlying the data-generation, we can never infer causalities from data in itself.

When a model trained on real-world (past) data is employed to influence future decisions, it is prudent to keep in mind that that model has learned to reflect the state of the world as it is/was, not how it should be.

### Extra parts: Regularization methods for $n < p$ data

In high dimensional data, we often have more features than observations. This is often called the $n < p$ scenario. In this situation, linear regression breaks down: the variance on the weight estimates blows up and the model will fail on unseen data. Both ridge regression and the lasso are valuable solutions here.

We will work with a dataset that contains spectral measurements on food samples. The target variables are the water, fat and protein content of the samples.


```python
# Read in the data
data = pd.read_csv('./meatNIR1000.csv', header=None)
colnames = pd.read_csv('./meatNIR1000.colnames', header=None)
data.columns = colnames.values[0]
W = data['water']
F = data['fat']
P = data['protein']

data = data.drop(['water', 'fat', 'protein'], axis=1)
print(data.shape)
data.head()
```

We have 240 observations, but 500 features, so we are in a $n < p$ setting:


```python
# let's plot some random data points

fig, ax = plt.subplots(figsize=(20,10))
for i in range(20):
    ax.plot(data.iloc[np.random.choice(range(len(data)))])
ax.set_xlabel('Wavelength').set_fontsize(20)
ax.set_ylabel('Absorbance').set_fontsize(20)
```

<div class="alert alert-success">
    <h2>Exercise</h2>
    <p>Try to fit a linear regression model to predict the fat content of a sample on the entire dataset (don't bother with cross validation / splitting / scaling, just fit a model). Look at the model coefficients (look up documentation to see how to) What happens?
    </p>
</div>


```python
from sklearn.linear_model import LinearRegression
X = data.values
y = F
```

<div class="alert alert-success">
    <h2>Exercise</h2>
    <p>Compare the performance of linear regression with ridge regression and with the lasso. Use 5-fold cross-validation to evaluate both models.</p>
</div>


```python
import warnings
warnings.filterwarnings('ignore')

from sklearn.linear_model import LassoCV, RidgeCV
# Select X and y
y = F
X = data.values

kf = KFold(n_splits=5)

linreg_scores = []
ridge_scores = []
lasso_scores = []

for train_index, test_index in kf.split(X):

print('Average validation fold R² of linear regression: {}'.format(np.mean(linreg_scores)))
print('Average validation fold R² of ridge regression: {}'.format(np.mean(ridge_scores)))
print('Average validation fold R² of lasso regression: {}'.format(np.mean(lasso_scores)))
```


```python

```

ee75bc473f746da57ddfc14e2cec5f8022e7744d
ipynb
Jupyter Notebook
predmod/lab6/PClab6_reg_2022.ipynb
gdewael/teaching
a78155041918422a843f31c863dd11e8afc5646a
MIT
predmod/lab6/PClab6_reg_2022.ipynb
gdewael/teaching
a78155041918422a843f31c863dd11e8afc5646a
MIT
predmod/lab6/PClab6_reg_2022.ipynb
gdewael/teaching
a78155041918422a843f31c863dd11e8afc5646a
MIT
Qwen/Qwen-72B
1. YES
2. YES