<jupyter_start><jupyter_text>Overview**This tutorial will focus on _Generalized Linear Models_**We will include some additional topics, including: * Attribute normalization to improve interpretability* L1 regularization* L2 regularization* (Cross)-validation**The lecture relies on the the following proficiencies and tools:*** [Python](https://www.python.org/) programming* Vector computations via [the numpy module](https://numpy.org/)* Data handling using [the pandas module](https://pandas.pydata.org/)* Plotting using [matplotlib](https://matplotlib.org/)* Training and using Machine Learning model via [scikit-learn](https://scikit-learn.org/)You will need them only if you plan to handle these tasks yourself  Normalization  Impact of Individual Attributes**We may be interested in understanding the impact of individual attributes*** E.g. is the price more impacted by age, distance from MRT stations, or what?   A Straightforward Approach**Let's start by repeating our basic linear regression approach**<jupyter_code>import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score

data = pd.read_csv('data/real_estate.csv', sep=',')
cols = data.columns
X = data[cols[:-1]]
y = np.log(data[cols[-1]])
X_tr, X_ts, y_tr, y_ts = train_test_split(X, y, test_size=0.34, random_state=42)
m = LinearRegression()
m.fit(X_tr, y_tr);

y_pred_tr, y_pred_ts = m.predict(X_tr), m.predict(X_ts)
print(f'R2: {r2_score(y_tr, y_pred_tr):.3} (training), {r2_score(y_ts, y_pred_ts):.3} (test)')<jupyter_output>R2: 0.691 (training), 0.645 (test)<jupyter_text>* We load the dataset, separate input & output, split the training and test set* The we fit the model and evaluate the results  A Straightforward Approach**We can think of looking at the coefficients**<jupyter_code>from matplotlib import pyplot as plt
idx = np.argsort(-np.abs(m.coef_)) # Sort indexes by decreasing absolute coefficient
scores = m.coef_[idx] # Coefficients, in the order defined by idx
plt.figure(figsize=(9, 3))
plt.bar(0.5 + np.arange(len(idx)), scores, tick_label=X.columns[idx]);<jupyter_output><empty_output><jupyter_text>* Intuition: the higher the (absolute) weight, the more important the attribute* Positive coefficient = price increase, negative coefficient = price decrease  A Straightforward Approach**There is a (huge) pitfall with this approach*** Each coefficient depends on the importance of an attribute* ...But also on its scale!An attribute with _larger values_ may end up with a _smaller coefficient_**We can preprocess the data so that all attributes have _similar scales_**A couple of common approaches:* Normalization (all values in the 0-1 range)$$x^\prime = \frac{x - \min(x)}{\max(x) - \min(x)}$$* Standardization (the mean becomes 0, the st. dev. becomes 1):$$x^\prime = \frac{x - \text{mean}(x)}{\text{std}(x)}$$  Normalization/Standardization**Scikit learn provides simple tools to take care of this step**The are called `MinMaxScaler` and `StandardScaler`:<jupyter_code>from sklearn.preprocessing import StandardScaler # for standardization

scaler = StandardScaler()

X_tr_s = scaler.fit_transform(X_tr)<jupyter_output><empty_output><jupyter_text>**On the training set:*** We calibrate the scaler (i.e. compute the mean and st. dev.)* The we apply the transformation<jupyter_code>X_ts_s = scaler.transform(X_ts)<jupyter_output><empty_output><jupyter_text>**On the test set we simply apply the transformation*** We do not want information about the test data to "leak" into the model  The Correct Approach**We can then repeat the process using the scaled data**<jupyter_code>m = LinearRegression()
m.fit(X_tr_s, y_tr);

y_pred_ts = m.predict(X_ts_s)
print(f'R2 on the test set: {r2_score(y_ts, y_pred_ts):.3}')<jupyter_output>R2 on the test set: 0.645<jupyter_text>* The quality of the predictions is the same as before* ...Since we have applied only a linear transformation to the input  The Correct Approach**Now we can build and inspect the same plot as before**<jupyter_code>idx = np.argsort(-np.abs(m.coef_)) # Sort indexes by decreasing absolute coefficient
scores = m.coef_[idx] # Coefficients, in the order defined by idx
plt.figure(figsize=(9, 3))
plt.bar(0.5 + np.arange(len(idx)), scores, tick_label=X.columns[idx]);<jupyter_output><empty_output><jupyter_text>* The order of the attributes has changed!* Positive/negative values = positive/negative effect if above average  LASSO  Attribute Selection**Let's look again at the last plot**<jupyter_code>plt.figure(figsize=(9, 3))
plt.bar(0.5 + np.arange(len(idx)), scores, tick_label=X.columns[idx]);<jupyter_output><empty_output><jupyter_text>* The "longitude" attribute has very low importance* Given that, is it worth keeping the attribute at all?* Collecting useless attributes has still a cost  Attribute Selection**Choosing which attributes to keep is known as _attribute selection_**There are multiple solution approaches* An exact (but costly) one: consider all possible sets and repeat training* A cheap (but heuristic) one: remove attributes with small (abs.) coefficients**In the case of linear regression, we have _a third option_**We can _act on the structure of the loss (error) function_ used during training:$$L(w) = \mathit{MSE}(w) + \alpha \|w\|_1$$* Where $\|w\|_1 = \sum_{j=1}^n |w_j|$* The additional term $\|w\|$ is introduced to make mode model "well behaved"* ...And for this reason it is called a _reguralizer_* In particular, since we are using the L1 norm, we call it an _L1 regularizer_  LASSO**At training time, we will need to solve**$$\text{argmin}_w \mathit{MSE}(w) + \alpha \|w\|_1$$For reasons not discussed in detail:* If raising a coefficient does not provide enough error improvement* ...Then the attribute will not be used at all* What "enough" means in practice depends on the value of $\alpha$**As a result, we naturally end-up doing attribute selection**The approach is known as _Least Absolute Shrinkage and Selection Operator_* ...Or LASSO for short :-)  Using LASSO**Scikit-learn provides a convenient LASSO implementation**<jupyter_code>from sklearn.linear_model import Lasso

m2 = Lasso(alpha=1.0)<jupyter_output><empty_output><jupyter_text>* When building the model, we need to pick a value for $\alpha$* Higher $\alpha$ values promote using fewer attributes* ...Bu they also reduce the prediction quality**For now, we will fix $\alpha$ via trial-and-error approach:**<jupyter_code>m2.fit(X_tr_s, y_tr);
y_pred_tr2 = m2.predict(X_tr_s)
print(f'R2 on the training set: {r2_score(y_tr, y_pred_tr2):.3}')<jupyter_output>R2 on the training set: 0.0<jupyter_text>* This is too low, and therefore $\alpha = 1.0$ is too large  Choosing $\alpha$<jupyter_code>m2 = Lasso(alpha=0.1)
m2.fit(X_tr_s, y_tr);
y_pred_tr2 = m2.predict(X_tr_s)
print(f'R2 on the training set: {r2_score(y_tr, y_pred_tr2):.3}')<jupyter_output>R2 on the training set: 0.549<jupyter_text>* The prediction quality is still much lower than the original one<jupyter_code>m2 = Lasso(alpha=0.05)
m2.fit(X_tr_s, y_tr);
y_pred_tr2 = m2.predict(X_tr_s)
print(f'R2 on the training set: {r2_score(y_tr, y_pred_tr2):.3}')<jupyter_output>R2 on the training set: 0.643<jupyter_text>* Now we are back to the original value (more or less)!  Impact of Each Attribute**Once we are satisfied with the $\alpha$ values, we can inspect the coefficients**<jupyter_code>idx = np.argsort(-np.abs(m2.coef_))
scores = m2.coef_[idx]
plt.figure(figsize=(9, 3))
plt.bar(0.5 + np.arange(len(idx)), scores, tick_label=X.columns[idx]);<jupyter_output><empty_output><jupyter_text>* The "longitude" attribute has been dropped* When there are many attributes, the LASSO method can be extremely useful  Ridge Regression  Overfitting**We have already observed how our model has a modest degree of overfitting**<jupyter_code>print(f'R2 on the training set (linear regression): {r2_score(y_tr, y_pred_tr):.3}')
print(f'R2 on the test set (linear regression): {r2_score(y_ts, y_pred_ts):.3}')<jupyter_output>R2 on the training set (linear regression): 0.691
R2 on the test set (linear regression): 0.645<jupyter_text>**Overfitting arises when:*** The model learns patterns in the training set* ...That do not repeat on unseen data (e.g. since they due to random sampling)**This becomes more critical with many attributes and scarce data*** A linear regression model with $n$ coefficients* ...Can interpolate perfectly any collections of $n$ examples...But typically, that means we learning noise  First, Let's Make It Worse**First, let's amplify the problem, so as to make any change more significant**In particular, we will reduce the amount of available training data<jupyter_code>X_tr2, X_ts2, y_tr2, y_ts2 = train_test_split(X, y, test_size=0.938, random_state=42)

scaler2 = StandardScaler()
X_tr_s2 = scaler2.fit_transform(X_tr2)
X_ts_s2 = scaler2.transform(X_ts2)
print(f'Training set size: {len(X_tr_s2)}')

m3 = LinearRegression()
m3.fit(X_tr2, y_tr2)

y_pred_tr3, y_pred_ts3 = m3.predict(X_tr2), m3.predict(X_ts2)
print(f'R2: {r2_score(y_tr2, y_pred_tr3):.3} (training), {r2_score(y_ts2, y_pred_ts3):.3} (test)')<jupyter_output>Training set size: 25
R2: 0.628 (training), 0.509 (test)<jupyter_text>* Having access to little data is not uncommon* Raw measurements are easy to obtain, but targets may be expensive to collect  Ridge Regression**We will show how to reduce overfitting by using _ridge regression_**...Which is obtained from linear regression by adding a quadratic reguralizer:$$\text{argmin}_w \mathit{MSE}(w) + \alpha \|w\|_2^2$$where $\|w\|_2^2 = \sum_{j=1}^2 w_j^2$* Since the $\|w\|_2$ is an L2 norm, this is called an _L2 regularizer_* It looks similar to the L1 regularizer in the LASSO approach* ...But it has a very different effect**This regularizer _prevents weights from having a large absolute value_*** Smaller (absolute) weights result in smoother behavior* ...Which typically leads to better generalization  Using Ridge Regression**Scikit-learn provides an implementation for Ridge Regression**<jupyter_code>from sklearn.linear_model import Ridge
m4 = Ridge(alpha=1.0)
m4.fit(X_tr_s2, y_tr2);<jupyter_output><empty_output><jupyter_text>It works similarly to the `Lasso` class* When we build the model, we need to specify a value for $\alpha$* ...Then we train the model as usual<jupyter_code>y_pred_tr4 = m4.predict(X_tr_s2)
y_pred_ts4 = m4.predict(X_ts_s2)
print(f'R2 on the training set: {r2_score(y_tr2, y_pred_tr4):.3}')
print(f'R2 on the test set: {r2_score(y_ts2, y_pred_ts4):.3}')<jupyter_output>R2 on the training set: 0.627
R2 on the test set: 0.534<jupyter_text>* The default value for $\alpha$ is 1.0* ...And it is already helping a bit the model on unseen data  Cross-Validation and Hyper-Parameter Tuning**This time, we will follow a systematic approach to pick a value for $\alpha$*** Our goal is to improve the model results on _unseen examples_* So we can choose $\alpha$ to minimize the error on examples _not in the training set_**It would be tempting to optimize $\alpha$ for the test set...**...But it would also be _very wrong_* The test set is our ultimate approach to assess generalization* For this, its content _stay unseen_ whenever we calibrate parameters**We can however split _yet another set_ on this purpose*** This would not be used to train the model* ...But only for tuning the training parameters (also called _hyper-parameters_)**This is called a _validation set_*** But there is a drawback: we have _even less data available for training_  Cross-Validation and Hyper-Parameter Tuning**As an alternative, we can use _$k$-fold cross-validation_:*** We partition the training data in $k$ sets, called _folds_* For $i = 1..k$:  - We separate the $i$-th fold  - We train a model on the remaining data  - We evaluate the prediction quality on the separated fold  - Then we put the fold back  During each iteration, the examples in the $i$-th folder are effectively unseen**At the end of the process we have:*** $k$ models (that typically we _don't_ use)* $k$ prediction quality evaluations (of which we can take the mean)We can use cross validation to _assess how good a model is with unseen data_  Applying Cross-Validation and Hyper-Parameter Tuning**Scikit learn allows us to automate this process with a few lines of code*** First we specify which parameters we want to tune* ...And which range we want to investigateWe do so by building a dictionary:<jupyter_code>param_ranges = {'alpha': np.logspace(0, 2, 100)}<jupyter_output><empty_output><jupyter_text>* The dictionary keys are the parameters to be tuned* The entries are the values to be testedThen we build a `GridSearchCV` object:<jupyter_code>from sklearn.model_selection import GridSearchCV

gscv = GridSearchCV(estimator=Ridge(), param_grid=param_ranges, cv=5)<jupyter_output><empty_output><jupyter_text>* The `cv` parameter specifies the number of folds  Applying Cross-Validation and Hyper-Parameter Tuning**The we can "train" the `GridSearchCV` object*** This will try all combinations of the values in parameter ranges* ...And pick the parameters leading to the best cross-validation score<jupyter_code>gscv.fit(X_tr_s2, y_tr2);<jupyter_output><empty_output><jupyter_text>The best parameters can be accessed via the `best_params_` field:<jupyter_code>gscv.best_params_<jupyter_output><empty_output><jupyter_text>Applying Cross-Validation and Hyper-Parameter Tuning**By default, the `GridSearchCV` object:*** Trains an estimator with the best parameters _on all data_* This is the model that we are going to useCalling `predict` on `GridSearchCV` calls the same method on the best estimator:<jupyter_code>y_pred_tr4 = gscv.predict(X_tr_s2)
y_pred_ts4 = gscv.predict(X_ts_s2)
print(f'R2 on the training set: {r2_score(y_tr2, y_pred_tr4):.3}')
print(f'R2 on the test set: {r2_score(y_ts2, y_pred_ts4):.3}')<jupyter_output>R2 on the training set: 0.516
R2 on the test set: 0.638
{'fraction_non_alphanumeric': 0.1039426523297491, 'fraction_numerical': 0.017777777777777778, 'mean_word_length': 5.852161100196463, 'pattern_counts': {'":': 0, '<': 76, '<?xml version=': 0, '>': 76, 'https://': 5, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '603634', 'license': 'MIT', 'n_tokens_mistral': 4578, 'n_tokens_neox': 4388, 'n_words': 1643, 'path': '3. Generalized Linear Models.ipynb', 'repo_name': 'lompabo/linear-models-in-ml', 'size': '13950'}