# 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
# Read in the data
import pandas as pd 

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

Run the code below to perform best subset selection with linear regression. The code implements algorithm 6.1 from the book as shown above. For each number of features $k$, the code will print the number of feature combinations that is explored. We fit a model for each combination, so there is quite some computation involved.


```python
import scipy as scipy
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import KFold, cross_validate
from itertools import combinations 
from matplotlib import pyplot as plt
%matplotlib inline

# Specify features and labels
X = concretedata.drop(['compressiveStrength'], axis=1)
y = concretedata['compressiveStrength']

# Dicts to store the feature combinations and their scores
scores_perNfeat = {'mean':[], 'max':[], 'var':[]}
best_feature_perNfeat = []
# Define splitter object outside loop to ensure same splits over loop
splitter = KFold(n_splits=5, shuffle=True, random_state=None)

# For k in 1,...,p
for k in np.arange(1,X.shape[1]+1):
    print('Number of combinations  with {} features: {}'.format(k, scipy.special.comb(X.shape[1],k)))
    best_score_k = 0
    scores_k_comb = []
    # Create and loop over all the possible feature combinations given k
    for c in combinations(X.columns, k):
        # Select features from data
        X_c = X[list(c)]
        # Define model
        model = LinearRegression()
        # Fit model and compute score on training data
        scores = cross_validate(model, X_c.values, y.values, cv=splitter)
        mean_cv_score = np.mean(scores['test_score'])
        # Save scores
        scores_k_comb.append(mean_cv_score)
        if mean_cv_score > best_score_k:
            best_features = c
            best_score_k = mean_cv_score
    # Save scores based on # features

    scores_perNfeat['max'].append(np.max(scores_k_comb))
    scores_perNfeat['var'].append(np.var(scores_k_comb))
    scores_perNfeat['mean'].append(np.mean(scores_k_comb))
    # Save best feature combination
    best_feature_perNfeat.append(best_features)

#Get variance of scores within grouped scores 
for i, score in enumerate(scores_perNfeat['max']):
    print("Best CV score using {} features: {}".format(i+1, score))

fig, ax = plt.subplots(1,1,figsize=(10,8), sharey=True)

ax.errorbar(np.arange(len(scores_perNfeat['mean']))+1, scores_perNfeat['mean'], scores_perNfeat['var'], label='mean score')
ax.scatter(np.arange(len(scores_perNfeat['max']))+1, scores_perNfeat['max'], alpha=0.5, label='max score')
ax.set_xlabel('#features')
ax.set_ylabel('R^2')
fig.legend()
print('Best feature combination: {}'.format(best_feature_perNfeat[np.argmax(scores_perNfeat['max'])]))
```

We will do the same thing for 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.


```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()
```


```python
from sklearn.linear_model import LogisticRegression
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA

# Specify features and labels
X = bacterialdata.drop(['species'], axis=1)
y = bacterialdata['species']

# Dicts to store the feature combinations and their scores
scores_perNfeat = {'mean':[], 'max':[], 'var':[]}
best_feature_perNfeat = []
# Define splitter object outside loop to ensure same splits over loop
splitter = KFold(n_splits=5, shuffle=True, random_state=None)

# For k in 1,...,p
for k in np.arange(1,X.shape[1]+1):
    print('Number of combinations  with {} features: {}'.format(k, scipy.special.comb(X.shape[1],k)))
    best_score_k = 0
    scores_k_comb = []
    # Create and loop over all the possible feature combinations given k
    for c in combinations(X.columns, k):
        # Select features from data
        X_c = X[list(c)]
        # Define model
        model = LogisticRegression() # or model = LDA()
        # Fit model and compute score on training data
        scores = cross_validate(model, X_c.values, y.values, cv=splitter, n_jobs=4)
        mean_cv_score = np.mean(scores['test_score'])
        # Save scores
        scores_k_comb.append(mean_cv_score)
        if mean_cv_score > best_score_k:
            best_features = c
            best_score_k = mean_cv_score
    # Save scores based on # features

    scores_perNfeat['max'].append(np.max(scores_k_comb))
    scores_perNfeat['var'].append(np.var(scores_k_comb))
    scores_perNfeat['mean'].append(np.mean(scores_k_comb))
    # Save best feature combination
    best_feature_perNfeat.append(best_features)

#Get variance of scores within grouped scores 
for i, score in enumerate(scores_perNfeat['max']):
    print("Best CV score using {} features: {}".format(i+1, score))

fig, ax = plt.subplots(1,1,figsize=(10,8), sharey=True)

ax.errorbar(np.arange(len(scores_perNfeat['mean']))+1, scores_perNfeat['mean'], scores_perNfeat['var'], label='mean score')
ax.scatter(np.arange(len(scores_perNfeat['max']))+1, scores_perNfeat['max'], alpha=0.5, label='max score')
ax.set_xlabel('#features')
ax.set_ylabel('accuracy')
fig.legend()
print('Best feature combination: {}'.format(best_feature_perNfeat[np.argmax(scores_perNfeat['max'])]))
```

<div class="alert alert-success">
    <h2>Exercise</h2>
    <p>Implement stepwise selection for features as explained at the start of the notebook</p>
</div>


```python
from sklearn.linear_model import LogisticRegression

# Specify features and labels
X = bacterialdata.drop(['species'], axis=1)
y = bacterialdata['species']

# Lists to store the feature combinations and their scores
scores_perNfeat = {'mean':[], 'max':[], 'var':[]}
best_feature_perNfeat = []

features = [] # start with empty set of features

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

# For k in 1,...,p
for k in np.arange(1,X.shape[1]+1):
    best_score_k = 0
    scores_k_comb = []
    # Create and loop over all the possible features to add
    for feature_idx in np.arange(0,X.shape[1]):
        ... # make sure that you don't features are not added that are already present in 'features'
        temp_features = features + [feature_idx]
        # Select features from data
        X_c = X.iloc[:, temp_features]
        # Define model
        model = LogisticRegression()
        # Fit model and compute score on training data
        scores = cross_validate(model, X_c.values, y.values, cv=splitter, n_jobs=4)
        mean_cv_score = np.mean(scores['test_score'])
        # Save scores
        scores_k_comb.append(mean_cv_score)
        if mean_cv_score > best_score_k:
            best_feature = ...
            best_score_k = ...
    features.append(best_feature)
    # Save scores based on # features
    scores_perNfeat['max'].append(np.max(scores_k_comb))
    scores_perNfeat['var'].append(np.var(scores_k_comb))
    scores_perNfeat['mean'].append(np.mean(scores_k_comb))
    # Save best feature combination
    best_feature_perNfeat.append(features[:])

#Get variance of scores within grouped scores 
for i, score in enumerate(scores_perNfeat['max']):
    print("Best training score using {} features: {}".format(i+1, score))

fig, ax = plt.subplots(1,1,figsize=(10,8), sharey=True)
ax.errorbar(np.arange(len(scores_perNfeat['mean']))+1, scores_perNfeat['mean'], scores_perNfeat['var'], label='tune_data')
ax.scatter(np.arange(len(scores_perNfeat['max']))+1, scores_perNfeat['max'], alpha=0.5)
ax.set_xlabel('#features')
ax.set_ylabel('accuracy')
fig.legend()
print('final set of features: {}'.format(set(best_feature_perNfeat[np.argmax(scores_perNfeat['max'])])))
best_feature_perNfeat
```

## 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.

### 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)
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 dividing them by their respective standard deviation: (see book p. 232)**

\begin{equation}
\tilde{x}_{ij} = \frac{x_{ij}}{\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. 


```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):
    
    ... # extract train and test splits from train_index & test_index
    
    ... # call and fit scaler
    
    linregmodel = LinearRegression()
    ridgemodel = RidgeCV(alphas=np.logspace(-3,3,100))
    

    ... # fit linear regression model and append scores to linreg_scores
    
    ... # fit ridgemodel and append scores to ridge_scores
    
    print('Regularization parameter: {}'.format(ridgemodel.alpha_))
    

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)))

# ** solution
```

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):
    ... # Repeat from previous but with X_expanded
    # consider using a larger range of alphas for RidgeCV, think for yourself why?
    

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)))
# ** solution
```


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
#scaler = StandardScaler(with_mean=False)
#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(with_mean=False)
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=(40,30), 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);
```

### 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
data.head()
```


```python
data = pd.read_csv('./meatNIR1000.csv',header=None)
colnames = pd.read_csv('./meatNIR1000.colnames', header=None)
data.columns = colnames.values[0]

data.head()
```


```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)
data.head()
```

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


```python
data.shape
```


```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. Look at the model coefficients. What happens?
    </p>
</div>


```python
... # import, fit

... # look at the first twentycoefficients of the model
```

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


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

#** solution
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):
    ... # Define X_train, X_test, y_train, y_test
    
    ... # scaler fit & transform
    
    ... # Fit linear, ridge and lasso 
    
    ... # Append test scores to the lists

    

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)))
#**solution
```

98145f9a171d50faef2839ec5c7d0e1a947e1471
ipynb
Jupyter Notebook
predmod/lab6/PClab6_reg.ipynb
gdewael/teaching
a78155041918422a843f31c863dd11e8afc5646a
MIT
predmod/lab6/PClab6_reg.ipynb
gdewael/teaching
a78155041918422a843f31c863dd11e8afc5646a
MIT
predmod/lab6/PClab6_reg.ipynb
gdewael/teaching
a78155041918422a843f31c863dd11e8afc5646a
MIT
Qwen/Qwen-72B
1. YES
2. YES