<jupyter_start><jupyter_text>Credit Balance Prediction _Motivated by_: https://github.com/Mashimo/datascience/blob/master/01-Regression/Regularisation.ipynbThe Credit data set records Balance (average credit card debt for a number of individuals) as well as several quantitative predictors: age, cards (number of credit cards), education (years of education), income (in thousands of dollars), limit (credit limit), and rating (credit rating). In addition to these quantitative variables, we also have four qualitative variables: Gender, Student (student status), Married (marital status), and Ethnicity (Caucasian, African American or Asian).We will need pandas as before. This time, we will use [scikit-learn](http://scikit-learn.org/stable/) directly to peform the linear regression, rather than `statsmodel`, which uses libraries such as `scikit-learn` in the background.<jupyter_code>import pandas as pd
from sklearn import linear_model
import os
if not os.path.exists('res'):
    os.makedirs('res')

dataFile = "data/Credit.csv"
data = pd.read_csv(dataFile, index_col=0) 

data.shape
data.columns
data.head()
data.describe()
data.info()<jupyter_output><empty_output><jupyter_text>Preprocessing The `linear_model` prefers to see a separate `X` matrix of predictors and `y` vector of dependent variable. Therefore, we take the dataframe `data` and split it accordingly.<jupyter_code>y = data['Balance'].copy() # this is the Credit Balance vector - our 'y' vector
X = data.copy()
X.drop(['Balance'], axis=1, inplace=True) # X is a copy of the dataframe but with the 'Balance' vector removed, leaving just the matrix of predictors `X`.<jupyter_output><empty_output><jupyter_text>Recoding the `Gender` variable to _indicator_ formAs presented in the data, the `Gender` variable is categorical (nominal) with two values `Male` and `Female`. Such variables are not suitable predictors for linear regression unless they are recoded as _indicator_ variables, with values either `0` or `1`.Arbitrarily, we can map ` Male` to `0` and `Female` to `1`. According to this mapping, the `Gender` variable could be interpreted as `IsFemale`. However, we could swap the mappings without changing the structure of the model, although the interpretation of the `Gender` coefficient would change to `IsMale`.Note the space character before "Male" in the mapping statement below as this is present in the data, and the strings need to match exactly.<jupyter_code>X.Gender = X.Gender.map({'Female':1, ' Male':0})<jupyter_output><empty_output><jupyter_text>Recoding the `Student` and `Married` variables in indicator formThese variables are more obviously binary-values, so they can be mapped using the same approach.<jupyter_code>X.Student = X.Student.map({'Yes':1, 'No':0})
X.Married = X.Married.map({'Yes':1, 'No':0})<jupyter_output><empty_output><jupyter_text>Recoding the 3-valued `Ethnicity` variable as _three_ indicator variablesThe `Ethnicity` variable has 3 values and so cannot be mapped to a single (1,1)-valued variable. Pandas provides a utility function to create indicator variables for each of those three values. These variables are added to the `X` matrix and the redundant `Ethnicity` variable should be removed from `X`, as its information is captured in the 3 new indicator variables.Optionally, these 3 indicator variables can be combined into two indicator variables. For example, let the `Asian` variable be defined implictly as the case wehn both `Caucasian` and `AfricanAmerican` take the value `0`. From a modelling perspective, this is fine. However, this can make interpretation slightly more difficult, so we choose not to combine them here.<jupyter_code>ethnicityIndicators = pd.get_dummies(X['Ethnicity'])
X = X.join(ethnicityIndicators)
X.drop(['Ethnicity'], axis=1, inplace=True)<jupyter_output><empty_output><jupyter_text>Single predictor modelsWith `scikit-learn`, we first create a template model, which we later configure for specific fits.<jupyter_code>model = linear_model.LinearRegression()<jupyter_output><empty_output><jupyter_text>Credit Balance by `Gender` With `scikit-learn`, there is no model formula as there was with `statsmodels`. Instead you need to use columns in the dataframe. Also note that the constant term is assumed (so there is no need to `add_constant`.`sklearn` objects come with a `fit` function. In this case, you need to provide $X$ (the dataframe of features) and $y$ the values to fit.<jupyter_code>model.fit(X[['Gender']], y)
[model.intercept_, model.coef_]<jupyter_output><empty_output><jupyter_text>We can interpret this as follows. The average credit balance is \$509.80\$ for males and \$509.80 + \$19.73 = \$529.53\$ for females.However this is probably a gross simplification, as we are ignoring all the other variables.  Exercise: How does the Credit Balance vary with each of the other categorical variables?  Forward selection Previously, in the case of the Diamonds and Advertising data, we used forward selection to find the best set of predictors. We used a manual search that worked well because the number of candidate predictors was quite low.For a larger set, as here, it makes sense to leave much of the work to a python program.We identify two operations:1. Finding the $R^2$ metric for each of a candidate set of predictors, and returning the predictor with the largest metric2. Accumulate the predictors found to date, and look for the next best predictor.The following Python function `findNextBestPredictor()` provides operation 1., and the next block of python provides operation 2.<jupyter_code>def findNextBestPredictor(X,foundPredictors):
    nP = X.shape[1] # number of columns in X
    allPredictors = list(X) # See https://stackoverflow.com/a/19483025
    predictorsToSearch = set(allPredictors) - set(foundPredictors)
    maxScore = 0 # can usually do better than this!
    for predictor in predictorsToSearch: # loop over all remaining columns (predictors) in X
        trialPredictors = set(foundPredictors)
        trialPredictors.add(predictor) # Add this predictor to the existing predictors
        XcolSubset = X.loc[:,trialPredictors] # all rows and just the trial predictors
        model.fit(XcolSubset, y) # fit the model to y
        score = model.score(XcolSubset, y)
        if score > maxScore: # identify the largest score and its associated predictor
            maxScore = score
            bestPredictorFound = predictor

    return (maxScore, bestPredictorFound)

findNextBestPredictor(X, list())
nP = X.shape[1]
scores = [0]
foundPredictors = list()

for i in range(nP): # loop over all columns (predictors) in X
    (score, bestPredictorFound) = findNextBestPredictor(X, foundPredictors)
    foundPredictors.append(bestPredictorFound)
    scores.append(score)

print(foundPredictors)
print(scores)<jupyter_output><empty_output><jupyter_text>The first features are Rating, Income and Student.  With these 3 features, the $R^2$ score is already at 0.95 and it increases only marginally thereafter, see plot below.<jupyter_code>import matplotlib.pyplot as plt

plt.title("Score versus predictors")
plt.xlabel('Number of predictors') 
plt.ylabel('R squared')
plt.plot(scores)
resFig = "res/rSqByPred.pdf"
plt.savefig(resFig)<jupyter_output><empty_output><jupyter_text>The first 3 features are the best predictors, the others do not add much value. As we add more and more parameters to our model its complexity increases. Generally, this results in increasing variance and decreasing bias, i.e. overfitting, so the model does not generalise as well from the training data to new data. We have two options:1. Use model selection techniques, e.g., use `Rating`, `Income` and `Student` as predictors, and ignore the remaining predictors as offering little explanatory power. We have described this approach before, in other notebooks.2. Use regularisation techniques: _ridge regression_ and the _lasso_  RegularisationIf we attempt to use too many predictors, we can run into multicollinearity issues, where several predictors compete to explain the data. The regression procedure will choose one of the set, and the remaining predictors in that multicollinear set have extremely high variance.This is not good news. We want the set of predictors to have roughly equal, small variance.If you imagine the function we wish to minimise (the sum of the squares of the residuals) as a function of the $\beta$ parameters, it will look like an inverted ridge, with more variance along the spine of the ridge than across it. We wish to constrain the OLS search to home in on the "best" point, and not to get distracted by running off along the ridge.This is the motivation for both ridge and lasso regression. Standardised scalingBefore we compare $\beta$ for different types of regression, we should note that the linear model looks after scaling in a natural way. Recall that $\hat{y}_i = \sum_j X(i,j) \beta_j $ can be written as $ \hat{y}_i = \sum_j (cX(i,j)) \frac{\beta_j}{c}$ for any nonzero $c$. Therefore, if we wish to ensure that the Euclidean norm of each column of $X$ is 1, or each column has the same variance, we can scale the $X$ matrix appropriately.In that regard, `scikit-learn` provides a `StandardScaler()` to derive $c$ and `scaler.fit_transform()` to apply this scaling to $X$.<jupyter_code>from sklearn import preprocessing

scaler = preprocessing.StandardScaler()
scaledXarray = scaler.fit_transform(X) # Note that X is no longer a dataframe after scaling :-(
scaledX = pd.DataFrame(scaledXarray, index=X.index, columns=X.columns) # Convert it back<jupyter_output><empty_output><jupyter_text>We first have a look at the size of the regression coefficients $\beta$:<jupyter_code>model.fit(scaledX, y)
predictors = scaledX.columns
coef = pd.Series(model.coef_,predictors).sort_values()
coef<jupyter_output><empty_output><jupyter_text>Note that there is a huge difference in the sensitivity of the regression across all the predictors, which is something we often see when there are too many predictors. `scikit-learn` helpfully offers a plot of the coefficients $\beta$ that shows this even more clearly, see below:<jupyter_code>fig = coef.plot(kind='bar', title="Model coefficients")
resFig = "res/coeffs.pdf"
plt.savefig(resFig)<jupyter_output><empty_output><jupyter_text>Interestingly, the equivalent plot for unscaled data showed that the $|\beta_{\rm Cards}|$ coefficient was much larger than the others and so it dominates the regression model. Yet we also know that, on its own, it adds little to the regression. Therefore it was trying to compensate for other predictors, rather than adding much of value in its own right. The scaling helped to reduce this problem but there are still significant differences in the size of different coefficients. Ridge RegressionThe ridge regression coefficients $\beta_{\rm (ridge)}$ minimize$\sum_i (y_i - X_i\beta_{\rm (ridge)})^T(y_i - X_i\beta_{\rm (ridge)}) + \lambda ||\beta_{\rm (ridge)}||_2$The first term is the ordinary least squares (OLS) estimate. The second term is a measure of the size of the $\beta$ coefficients. If we minimise both together, we are looking for the $\beta$ coefficients that (simultaneously) make the regression fit close to the data _and_ have small values (measured as their Euclidean distance from the origin in $\beta$-space).The $\lambda$ multiplier weights the penalty term (the one that shrinks the coefficients) relative to the normal OLS term. Small values of $\lambda$ favour the OLS objective; larger ones favour the shrinkage objective. The easiest way to find a "good" $\lambda$ and hence a good overall fit is by cross validation.![Ridge Regression: geometrical interpretation](pic/ridgeRegression.png)In the visualisation above, the OLS component of the objective function has the characteristic contours of a ridge with its minimum indicated with a red dot. The regularisation component is represented by the disk, with its minimum represented by the red dot at the origin. The overall objective depends on the scaling of the disk and where it intersects with the ridge, represented by the third red dot above. As can be seen,* the ridge estimate does not usually coincide with the OLS estimate, and* it generally has as many dimensions (2 nonzero values in this case) as the OLS estimate.Note that the combined objective function used in ridge regression is _not_ scale-invariant. Hence $X$ should be scaled appropriately so that different fits can be compared across $\lambda$, without the $\lambda$ values themselves affecting the scaling of the $\beta$ coefficients.Here is an example using `scikit-learn`'s `Ridge()` function. For this example, we set $\lambda$ (specified as `alpha` in the call to `Ridge`) to 50.<jupyter_code>from sklearn.linear_model import Ridge

# assign lambda (or alpha as is used in Ridge())
ridgeReg = Ridge(alpha=50)

# Now compute the ridge regression using X, the full matrix of predictors, and y
ridgeReg.fit(scaledX, y)

ridgeReg.score(scaledX, y)
coef = pd.Series(ridgeReg.coef_,predictors).sort_values()
coef<jupyter_output><empty_output><jupyter_text>As you can see, the $R^2$ score is reduced, reflecting the fact that the OLS criterion is no longer the only criterion and so the bias is greater, but it is also clear that the size of the coefficients has also shrunk and the predictors with the largest $|\beta|$ tend to be those we know from previous analysis are the more important predictors anyway. How does the penalty weight $\lambda$ affect the coefficients $\mathbf{\beta}$?We can collect the coefficients $\beta$ for a selection of $\lambda$ values.<jupyter_code>import numpy as np

lambdas = [0.01, 1, 100, 10000]
coefficients = np.zeros(shape=(nP, len(lambdas)))
i=0
for l in lambdas:
    ridgeReg = Ridge(alpha=l)
    ridgeReg.fit(scaledX, y)
    coefficients[:,i] = ridgeReg.coef_
    i += 1<jupyter_output><empty_output><jupyter_text>Now we can plot the coefficients as $\lambda$ increases, highlighting those of interest:<jupyter_code>fig, ax = plt.subplots()

plt.title("Ridge regularisation")
plt.xlabel("lambda")
plt.ylabel("standardised coefficients")
styles=['-','--','-.',':']

labels = ['0.01','','1','', '100','', '1e4','']
ax.set_xticklabels(labels) # custom x labels

chosenPredictors = {"Income", "Rating", "Student"}
for i in range(0,nP):
    s = styles[i % len(styles)]
    if predictors[i] in chosenPredictors:
        plt.plot(coefficients[i], label=predictors[i], linestyle=s)
    else:
        plt.plot(coefficients[i])

plt.legend(loc='best')
resFig = "res/ridge.pdf"
plt.savefig(resFig)<jupyter_output><empty_output><jupyter_text>Each curve corresponds to the ridge regression coefficient estimate for one of the variables, plotted as a function of $\lambda$.At the left side of the plot, $\lambda=0.01$ is essentially zero, and so the corresponding ridge coefficient estimates are the same as the usual least squares estimates. But as $\lambda$ increases, the ridge coefficient estimates shrink towards zero. When $\lambda$ is extremely large, all of the ridge coefficient estimates are basically zero; this corresponds to the null regression model that contains no predictors and just the regularisation component of the objective function.In this plot, the income, rating, and student variables are highlighted, since these variables are the ones to have the largest coefficient estimates. The orange trace is interesting - it appears to represent the intercept term $\beta_0$.As noted in the visualisation of ridge regression earlier, ridge regression generally does not reduce the number of active predictors - it just drives many of them near to zero as $\lambda$ increases.  The lassoThe lasso regression coefficients $\beta_{\rm (lasso)}$ minimize$\sum_i (y_i - X_i\beta_{\rm (lasso)})^T(y_i - X_i\beta_{\rm (lasso)}) + \lambda ||\beta_{\rm (lasso)}||_1$The first term is the ordinary least squares (OLS) estimate. The second term is a measure of the size of the $\beta$ coefficients. Here the Manhattan norm is used, where the Euclidean norm was used for ridge regression. If we minimise both together, we are looking for the $\beta$ coefficients that (simultaneously) make the regression fit close to the data _and_ have small values (measured as their Manhattan distance from the origin in $\beta$-space).![Lasso Regression: geometrical interpretation](pic/lassoRegression.png)In the visualisation above, the OLS component of the objective function has the characteristic contours of a ridge with its minimum indicated with a red dot. The regularisation component is represented by the rotated square, with its minimum represented by the red dot at the origin. The overall objective depends on the scaling of the square and where it intersects with the ridge, represented by the third red dot above. As can be seen,* the lasso estimate does not usually coincide with the OLS estimate, and* it generally has fewer dimensions (1 in this case, since $\beta_0 = 0$) as the OLS estimate (which has 2 in this case), because the intersection point is generally at a corner of the hypercube (square in this case) bounding the penalty component of the objective function.Using lasso regression from `sklearn` is very similar to how ridge regression is used. Again, we apply lasso regression to a range of $\lambda$ values and review how the $\beta$ coefficients tend to 0 as $\lambda$ increases.<jupyter_code>from sklearn.linear_model import Lasso

lambdas = [0.01, 1, 10, 50, 100, 200, 500, 1000]
coefficients = np.zeros(shape=(nP, len(lambdas)))
i=0
for l in lambdas:
    lassoReg = Lasso(alpha=l)
    lassoReg.fit(scaledX, y)
    coefficients[:,i] = lassoReg.coef_
    i += 1
fig, ax = plt.subplots()

plt.title("Lasso Regularisation")
plt.xlabel("lambda")
plt.ylabel("standardised coefficients")
styles=['-','--','-.',':']

labels = ['0.01','1','10','50', '100', '200','500','1000']

# See https://stackoverflow.com/a/63755285/1988855

# Outdated code, generates warning but works
ax.set_xticklabels(labels) # custom x labels

chosenPredictors = {"Income", "Rating", "Student"}
for i in range(0,12):
    s = styles[i % len(styles)]
    if predictors[i] in chosenPredictors:
        plt.plot(coefficients[i], label=predictors[i], linestyle=s)
    else:
        plt.plot(coefficients[i])

plt.legend(loc='best')
resFig = "res/lasso.pdf"
plt.savefig(resFig)<jupyter_output><empty_output><jupyter_text>$\lambda$ has a similar job in both ridge and lasso regression, but the behaviour is different. Look at how the optimum `Rating` actually _increases_ over the range $10 < \lambda < 500$ while all the other $\beta$ coefficients were driven to zero.The challenge is to choose this tuning parameter ($\lambda$). Selecting the Tuning Parameter lambdaCross-validation provides a simple way to tackle this problem.1. Choose a grid of $\lambda$ values2. Compute the cross-validation error for each value of $\lambda$.3. Select the tuning parameter $\lambda$ for which the cross-validation error is smallest.4. Re-fit using all the observations and the selected value of $\lambda$. Cross-Validation with $n=10$ splits is relatively straightforward to implement, although it requires significant computation. The data is randomly partitioned into $n$ subsets. $n-1$ of these subsets are used for training, and the remaining subset is used for validation. The prediction error (root-mean-square residual error) is computed for this validation subset. The procedure is repeated $n=10$ times in all, so that each of the $n$ subsets has a turn as the validation set. The average prediction error is computed and returned as the CV score. The python code below shows how this can be achieved using `sklearn`.<jupyter_code>from sklearn.model_selection import cross_val_score, KFold

lambdas = np.linspace(500,0.01,num=50) # Step 1
scoresCV = []
for l in lambdas:
    lassoReg = Lasso(alpha=l)
    lassoReg.fit(scaledX, y)    
    
    scoreCV = cross_val_score(lassoReg, scaledX, y, scoring='neg_mean_squared_error', 
                             cv=KFold(n_splits=10, shuffle=True,
                                            random_state=1))
    scoresCV.append(np.mean(scoreCV))
# Step 2 complete
plt.title("Lambda tuning for Lasso Regularisation")
plt.xlabel("lambda")
plt.ylabel("Cross-Validation error (MSE)")

plt.plot(lambdas,scoresCV)
resFig = "res/lambdaLassoCV.pdf"
plt.savefig(resFig)<jupyter_output><empty_output><jupyter_text>The plot shows how the CV error changes with $\lambda$. It looks like the CV error is minimised near $\lambda = 0$.Using `LassoCV`, which applies cross-validation to the Lasso, we can find the $\lambda$ that minimises the CV error.<jupyter_code>from sklearn.linear_model import LassoCV

lassoCV = LassoCV()
lassoCV.fit(scaledX, y)

lassoCV.alpha_
lassoReg = Lasso(lassoCV.alpha_)
lassoReg.fit(scaledX, y)
lassoReg.score(scaledX, y)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08441931774626293, 'fraction_numerical': 0.009486393254120352, 'mean_word_length': 5.058925979680697, 'pattern_counts': {'":': 0, '<': 75, '<?xml version=': 0, '>': 74, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '119054', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 6040, 'n_tokens_neox': 5631, 'n_words': 2605, 'path': 'topics/04-Regression/02-Practical/files/Case4_CreditBalance.ipynb', 'repo_name': 'DataMining-2021S2/live', 'size': '20872'}