# Linear Regression and Regularization
Loosely following Chapter 10 of Python Machine Learning 3rd Edition, Raschka

>Disclaimer: Regression is a huge field. It is impossible to cover it all in one class (or even two).



[Image Source](https://scikit-learn.org/stable/auto_examples/linear_model/plot_ols.html)

>For a sense of the depth and potential complexity of regression models, see [Mostly Harmless Econometrics](https://www.mostlyharmlesseconometrics.com)

# Ordinary Least Squares (OLS) Linear Regression

>[All models are wrong, but some are useful.](https://en.wikipedia.org/wiki/All_models_are_wrong)
<br><br>George Box

Linear regression is one of the most popular, widely used, and foundational concepts in statistics, econometrics, and machine learning.  

Boils down having a numeric target value ($y$) we want to either predict or understand the variance drivers. We use data ($X$) we think impacts our target to understand the underlying **linear** relationship.

### Big Assumption: 
- The regression function $E(Y|X)$ is **linear** in the inputs $X_{1},\dots,X_{p}$
- Transformations can be applied to satisfy that requirment.

Typically will see it expressed as $y = \beta X$, or formally:  

$$
f(X)=\beta_{0}+\sum{X{j}\beta_{j}}
$$

- Linear model assumes the function is linear or reasonably linear.
- The true $\beta_{j}$'s are unknown parameters (coefficients/weights). 
- The features must be able to be represented within a non-null numeric matrix.


### Goal - Minimize the mean-squared error. Why?
Residuals will be positive and negative, need to penalize negative and positive equally.

Sum of errors: $\epsilon_{1} + \epsilon_{2} + \dots + \epsilon_{n}$  
RSS: $\epsilon_{1}^2 + \epsilon_{2}^2 + \dots + \epsilon_{n}^2$  
MSE: $\frac{1}{N}\sum{\epsilon_{i}^2}$


Most statistical programs solve for the $\beta$ values using plain old linear alegbra, in what is called the closed-form:  

$\hat\beta = (X^TX)^{-1}X^{T}y$

### Closed Form Derivation

$$RSS(\beta)=\sum{(y_{i}-f(x_{i}))^2}$$

$$=\sum(y_{i}-\beta_{0}-\sum{x_{ij}\beta{j}})^2$$

$(x_i,y_i)$ should be independent from other $(x_i,y_i)$'s  
- Time series models violate this without special treatment

We are seeking a $f(X)$ that minimizes the sum of squared residuals from $Y$:




[Image source: Elements of Statistical Learning, Figure 3.1](https://www.statlearning.com)

$$
RSS(\beta)=(y-X\beta)^T(y-X\beta)
$$

#### Differentiating:
$$
\frac{dRSS}{d\beta}=-2X^T(y-X\beta)
$$

#### And again:
$$
\frac{d^2RSS}{d\beta d \beta^T}=2X^TX
$$

#### Setting the first derivative to zero:
$$ 
X^T(y-X\beta)=0
$$

#### And we get:
$$
\hat{\beta}=(X^TX)^{-1}X^Ty
$$

#### And our predicted values are:
$$
\hat{y}=X\hat{\beta}
$$

#### And relates to $y$ by:
$$
y = \hat{y} + \epsilon =X\hat{\beta}+\epsilon 
$$

> Unique solution means we can derive with pure linear algebra, i.e., no need to use convergence algorithms.

## Slope and Intercept
### Remember in its simple form: $$y=mx+b$$



[Image source](https://en.wikipedia.org/wiki/Linear_regression#/media/File:Linear_least_squares_example2.png)

> Since we need an estimate for the intercept, we'll need to manually add the constant. Not all implementations do this automatically, e.g., statsmodels.

- Intercept: where the line go through the y-axis.  
- Slope: for a 1-unit change in x, y will increase by $\beta$


# Example - Credit Data
[Data from Elements of Statistical Learning](https://www.statlearning.com)


```python
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

credit = pd.read_csv('data/islr-credit.csv')
credit = credit.iloc[:, 1:]
credit.head()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Income</th>
      <th>Limit</th>
      <th>Rating</th>
      <th>Cards</th>
      <th>Age</th>
      <th>Education</th>
      <th>Gender</th>
      <th>Student</th>
      <th>Married</th>
      <th>Ethnicity</th>
      <th>Balance</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>14.891</td>
      <td>3606</td>
      <td>283</td>
      <td>2</td>
      <td>34</td>
      <td>11</td>
      <td>Male</td>
      <td>No</td>
      <td>Yes</td>
      <td>Caucasian</td>
      <td>333</td>
    </tr>
    <tr>
      <th>1</th>
      <td>106.025</td>
      <td>6645</td>
      <td>483</td>
      <td>3</td>
      <td>82</td>
      <td>15</td>
      <td>Female</td>
      <td>Yes</td>
      <td>Yes</td>
      <td>Asian</td>
      <td>903</td>
    </tr>
    <tr>
      <th>2</th>
      <td>104.593</td>
      <td>7075</td>
      <td>514</td>
      <td>4</td>
      <td>71</td>
      <td>11</td>
      <td>Male</td>
      <td>No</td>
      <td>No</td>
      <td>Asian</td>
      <td>580</td>
    </tr>
    <tr>
      <th>3</th>
      <td>148.924</td>
      <td>9504</td>
      <td>681</td>
      <td>3</td>
      <td>36</td>
      <td>11</td>
      <td>Female</td>
      <td>No</td>
      <td>No</td>
      <td>Asian</td>
      <td>964</td>
    </tr>
    <tr>
      <th>4</th>
      <td>55.882</td>
      <td>4897</td>
      <td>357</td>
      <td>2</td>
      <td>68</td>
      <td>16</td>
      <td>Male</td>
      <td>No</td>
      <td>Yes</td>
      <td>Caucasian</td>
      <td>331</td>
    </tr>
  </tbody>
</table>
</div>



## Find Function so $Rating=f(Limit)$
We'll need to convert the pandas objects to numpy arrays.


```python
credit.plot.scatter(x='Limit', y='Rating')
plt.show()
```

## Regression using closed-form:


```python
X = np.array(credit['Limit']).reshape(-1,1)
y = np.array(credit['Rating']).reshape(-1,1)

X.shape, y.shape
```




    ((400, 1), (400, 1))



And since we are going to implement a version, we'll need to manually add the constant for the intercept. Why?

$y=\beta_{0}(1)+\beta_{i}x_{i}$


```python
from numpy.linalg import inv

'''
- Manually adding the constant
- Sometimes this is done via the API (check the docs)
''' 
const = np.ones(shape=y.shape)
mat = np.concatenate( (const, X), axis=1)

# first 5 examples
mat[:5,:]
```




    array([[1.000e+00, 3.606e+03],
           [1.000e+00, 6.645e+03],
           [1.000e+00, 7.075e+03],
           [1.000e+00, 9.504e+03],
           [1.000e+00, 4.897e+03]])



## Betas
We have a feature matrix that has 2 columns, so we'll get estimate for the constant ($\beta_{0}$) and the credit limit ($\beta_{1}$).

### Calculate the coefficient estimates
Recall $\hat\beta = (X^TX)^{-1}X^{T}y$


```python
betas = inv(mat.transpose().dot(mat)).dot(mat.transpose()).dot(y)
b0, b1 = betas

print(f'Beta 0: {np.round(b0[0],3)}')
print(f'Beta 1: {np.round(b1[0],3)}')
```

    Beta 0: 38.492
    Beta 1: 0.067


### Predict $\hat{y}$ and plot the fit

$$
\begin{equation}
\hat{y}=\hat{\beta}X=
\hat{\beta_{0}}
\begin{pmatrix}
1 \\
\dots \\
1
\end{pmatrix}+
\hat{\beta_{1}}
\begin{pmatrix}
3606 \\
\dots \\
5524
\end{pmatrix}
\end{equation}
$$


```python
yhat = mat.dot(betas)

plt.plot(X, y, 'bo')
plt.plot(X, yhat, 'r')
plt.xlabel('Credit Limit')
plt.ylabel('Credit Rating')
plt.title('$Rating=f(Limit)$', loc='left')
plt.show()
```

# Quantifying fit with metrics
Common metrics:

### $R^2$ [Wikipedia](https://en.wikipedia.org/wiki/Coefficient_of_determination)
$$1 - \frac{\sum (\hat {y}-y)^{2}}{\sum ({\bar y}-y)^{2}}$$

### Mean squared error (MSE) [Wikipedia](https://en.wikipedia.org/wiki/Mean_squared_error)
$$\frac{\sum (\hat {y}-y)^{2}}{n}$$

### Mean Absolute Error (MAE) [Wikipedia](https://en.wikipedia.org/wiki/Mean_absolute_error)
$$1/n\sum |\hat {y}-y|$$

### Root mean squared error (RMSE) [Wikipedia](https://en.wikipedia.org/wiki/Root_mean_square)
$$\sqrt \frac{\sum (\hat {y}-y)^{2}}{n}$$

#### Notes:
- $r^2$ expresses the percent of variance explained (bound between 0% and 100%).  
- RMSE expresses the variance in unit terms.  
- MSE/MAE are heavily influenced by outliers.  
- Usually RMSE is chosen for optimizing since it's an unbiased estimator.
- If there are a lot of outliers, MAE may be a better choice.

### Further reading:
[$r^2$ vs. RMSE](https://www.statology.org/rmse-vs-r-squared/)

# Intrepretability
A nice property of linear regression is the relatively simple intrepretations that can be drawn. The implementation in statsmodels includes all the output needed to evaluate and interpret model output.

[statsmodels OLS regression](https://www.statsmodels.org/stable/regression.html)


```python
# betas we calculated:
betas.reshape(1,-1)[0]
```




    array([38.49178871,  0.06682326])



statsmodels output:


```python
import statsmodels.api as smf

simpleModel = smf.OLS(y, mat).fit()
print(simpleModel.summary())
```

                                OLS Regression Results                            
    ==============================================================================
    Dep. Variable:                      y   R-squared:                       0.994
    Model:                            OLS   Adj. R-squared:                  0.994
    Method:                 Least Squares   F-statistic:                 6.348e+04
    Date:                Fri, 31 Dec 2021   Prob (F-statistic):               0.00
    Time:                        09:51:08   Log-Likelihood:                -1568.1
    No. Observations:                 400   AIC:                             3140.
    Df Residuals:                     398   BIC:                             3148.
    Df Model:                           1                                         
    Covariance Type:            nonrobust                                         
    ==============================================================================
                     coef    std err          t      P>|t|      [0.025      0.975]
    ------------------------------------------------------------------------------
    const         38.4918      1.397     27.555      0.000      35.746      41.238
    x1             0.0668      0.000    251.949      0.000       0.066       0.067
    ==============================================================================
    Omnibus:                        7.074   Durbin-Watson:                   2.077
    Prob(Omnibus):                  0.029   Jarque-Bera (JB):                5.177
    Skew:                           0.155   Prob(JB):                       0.0751
    Kurtosis:                       2.537   Cond. No.                     1.20e+04
    ==============================================================================
    
    Notes:
    [1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
    [2] The condition number is large, 1.2e+04. This might indicate that there are
    strong multicollinearity or other numerical problems.


## Results Discussion:
- Overall $r^2$ of 99% - very strong linear relationship as we saw in the initial plot we created.  
- You can ignore the t-statistic on the constant - it isn't meaningful.  
- The t-statistic for the credit limit (x1) is very high (252), with a [p-value](https://en.wikipedia.org/wiki/P-value) of essentially 0 - recall a p-value is the probably that the effect we are seeing is by random chance.  
- We can conclude there is a very strong linear relationship.

### Further reading:
[Rethinking p-values](https://www.vox.com/2016/3/15/11225162/p-value-simple-definition-hacking)  
[Econometrics](https://en.wikipedia.org/wiki/Econometrics)

# A Note on p-Hacking



>Taken from "Having Too Little or Too Much Time Is Linked to Lower Subjective Well-Being", Sharif et al.

# scikit-learn's LinearRegression
[scikit-learn's LinearRegression](https://scikit-learn.org/stable/modules/linear_model.html#ordinary-least-squares)

This will probably be the implementation you'd want to use for comparing various regression models for prediction problems since the API will be similar.

scikit-learn is geared more towards prediction and won't have some of the friendly output that statsmodels has - if you are going for an intrepretation, you may be better off using statsmodels.


```python
from sklearn.linear_model import LinearRegression

# we added the intercept manually, so turn that option off
skOLS = LinearRegression(fit_intercept=False).fit(mat,y)
skOLS.coef_
```




    array([[38.49178871,  0.06682326]])



#### And they all provide the same coefficient/model estimates - assuming the same set-up.


```python
print(betas.reshape(1,-1)[0])
print(simpleModel.params)
print(skOLS.coef_[0])
```

    [38.49178871  0.06682326]
    [38.49178871  0.06682326]
    [38.49178871  0.06682326]


# Weakness - Outlier Sensitivity
There aren't any outliers in the data, I'm going to introduce one.


```python
matCopy = mat.copy()
matCopy[4, 1] = matCopy[4, 1] + 150000

outModel = smf.OLS(y, matCopy).fit()
print(outModel.summary())
```

                                OLS Regression Results                            
    ==============================================================================
    Dep. Variable:                      y   R-squared:                       0.086
    Model:                            OLS   Adj. R-squared:                  0.084
    Method:                 Least Squares   F-statistic:                     37.54
    Date:                Wed, 29 Dec 2021   Prob (F-statistic):           2.16e-09
    Time:                        17:41:20   Log-Likelihood:                -2565.7
    No. Observations:                 400   AIC:                             5135.
    Df Residuals:                     398   BIC:                             5143.
    Df Model:                           1                                         
    Covariance Type:            nonrobust                                         
    ==============================================================================
                     coef    std err          t      P>|t|      [0.025      0.975]
    ------------------------------------------------------------------------------
    const        325.3864      8.837     36.820      0.000     308.013     342.760
    x1             0.0058      0.001      6.127      0.000       0.004       0.008
    ==============================================================================
    Omnibus:                       42.701   Durbin-Watson:                   2.058
    Prob(Omnibus):                  0.000   Jarque-Bera (JB):              192.528
    Skew:                           0.295   Prob(JB):                     1.56e-42
    Kurtosis:                       6.347   Cond. No.                     1.12e+04
    ==============================================================================
    
    Notes:
    [1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
    [2] The condition number is large, 1.12e+04. This might indicate that there are
    strong multicollinearity or other numerical problems.



```python
yhat_out = outModel.predict(matCopy)

plt.plot(matCopy[:,1], y, 'bo')
plt.plot(matCopy[:,1], yhat_out, 'r')
plt.xlabel('Credit Limit')
plt.ylabel('Actual / Expected Credit Rating')
plt.show()
```

## Why the sensitivity?
Recall we are minimizing the sum of squared residuals. That really big outlier is going to a lot of influence.  

## How to combat?
- [RANdom SAmple Consensus (RANSAC)](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.RANSACRegressor.html)  
- Replace or remove the outliers. 


## RANSAC
- Select random samples. 
- Tests non-sample points and creates a inlier list (opposite of outlier).  
- Refits models with all inliers.  
- Evaluates error.  
- Stops or repeats until iterations/threshold met.
- **Not guaranteed to get same answer each time* - why?**

> [More details](https://scikit-learn.org/stable/modules/linear_model.html#ransac-regression)



[Image Source](https://scikit-learn.org/stable/auto_examples/linear_model/plot_ransac.html)


```python
from sklearn.linear_model import RANSACRegressor

ransac = RANSACRegressor().fit(matCopy, y)
yhat_ransac = ransac.predict(matCopy)

plt.plot(matCopy[:,1], y, 'bo')
plt.plot(matCopy[:,1], yhat_ransac, 'r')
plt.xlabel('Credit Limit')
plt.ylabel('Actual / Expected Credit Rating')
plt.show()
```

# Strength: Robust to Overfitting
>Simple is better than complex.

### No Overfitting
 

### Severe Overfitting


# Multiple Regression
Instead of an $mx1$ input matrix, we'll have $mxn$.

### $y = w_{0} + w_{1}x_{1} + \dots + w_{m}x_{m} = \sum{w_{i}x_{i}}=w^{T}x$

Coefficients still reference the effect on $y$ of a 1-unit change to a $x$ - all else held constant.

Example with the California housing dataset:

## California Housing


```python
import pandas as pd
import numpy as np
from sklearn.datasets import fetch_california_housing
from matplotlib import pyplot as plt
import warnings
%matplotlib inline

with warnings.catch_warnings():
    warnings.filterwarnings("ignore")
    X, y = fetch_california_housing(return_X_y=True, as_frame=True)
    x0, x1 = X.shape

print(f'Rows: {x0:,}\nFeatures: {x1}')
```

    Rows: 20,640
    Features: 8



```python
california_df = pd.concat([X,y], axis=1)
california_df.head()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>MedInc</th>
      <th>HouseAge</th>
      <th>AveRooms</th>
      <th>AveBedrms</th>
      <th>Population</th>
      <th>AveOccup</th>
      <th>Latitude</th>
      <th>Longitude</th>
      <th>MedHouseVal</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>8.3252</td>
      <td>41.0</td>
      <td>6.984127</td>
      <td>1.023810</td>
      <td>322.0</td>
      <td>2.555556</td>
      <td>37.88</td>
      <td>-122.23</td>
      <td>4.526</td>
    </tr>
    <tr>
      <th>1</th>
      <td>8.3014</td>
      <td>21.0</td>
      <td>6.238137</td>
      <td>0.971880</td>
      <td>2401.0</td>
      <td>2.109842</td>
      <td>37.86</td>
      <td>-122.22</td>
      <td>3.585</td>
    </tr>
    <tr>
      <th>2</th>
      <td>7.2574</td>
      <td>52.0</td>
      <td>8.288136</td>
      <td>1.073446</td>
      <td>496.0</td>
      <td>2.802260</td>
      <td>37.85</td>
      <td>-122.24</td>
      <td>3.521</td>
    </tr>
    <tr>
      <th>3</th>
      <td>5.6431</td>
      <td>52.0</td>
      <td>5.817352</td>
      <td>1.073059</td>
      <td>558.0</td>
      <td>2.547945</td>
      <td>37.85</td>
      <td>-122.25</td>
      <td>3.413</td>
    </tr>
    <tr>
      <th>4</th>
      <td>3.8462</td>
      <td>52.0</td>
      <td>6.281853</td>
      <td>1.081081</td>
      <td>565.0</td>
      <td>2.181467</td>
      <td>37.85</td>
      <td>-122.25</td>
      <td>3.422</td>
    </tr>
  </tbody>
</table>
</div>



## Run the regression

Need to add a constant/intercept term manually as statsmodels doesn't handle this automatically.


```python
california_df['const'] = 1

featureNames = [x for x in california_df.columns if x != 'MedHouseVal']
featureNames = ['const'] + list(featureNames)

print(featureNames)
```

    ['const', 'MedInc', 'HouseAge', 'AveRooms', 'AveBedrms', 'Population', 'AveOccup', 'Latitude', 'Longitude', 'const']


> Double and triple check that you properly separated your target variable from your features!


```python
import statsmodels.api as smf
statsModelsCoefs = smf.OLS(california_df['MedHouseVal'], california_df[featureNames]).fit()

print(statsModelsCoefs.summary())
```

                                OLS Regression Results                            
    ==============================================================================
    Dep. Variable:            MedHouseVal   R-squared:                       0.606
    Model:                            OLS   Adj. R-squared:                  0.606
    Method:                 Least Squares   F-statistic:                     3970.
    Date:                Wed, 29 Dec 2021   Prob (F-statistic):               0.00
    Time:                        17:41:20   Log-Likelihood:                -22624.
    No. Observations:               20640   AIC:                         4.527e+04
    Df Residuals:                   20631   BIC:                         4.534e+04
    Df Model:                           8                                         
    Covariance Type:            nonrobust                                         
    ==============================================================================
                     coef    std err          t      P>|t|      [0.025      0.975]
    ------------------------------------------------------------------------------
    const        -18.4710      0.329    -56.067      0.000     -19.117     -17.825
    MedInc         0.4367      0.004    104.054      0.000       0.428       0.445
    HouseAge       0.0094      0.000     21.143      0.000       0.009       0.010
    AveRooms      -0.1073      0.006    -18.235      0.000      -0.119      -0.096
    AveBedrms      0.6451      0.028     22.928      0.000       0.590       0.700
    Population -3.976e-06   4.75e-06     -0.837      0.402   -1.33e-05    5.33e-06
    AveOccup      -0.0038      0.000     -7.769      0.000      -0.005      -0.003
    Latitude      -0.4213      0.007    -58.541      0.000      -0.435      -0.407
    Longitude     -0.4345      0.008    -57.682      0.000      -0.449      -0.420
    const        -18.4710      0.329    -56.067      0.000     -19.117     -17.825
    ==============================================================================
    Omnibus:                     4393.650   Durbin-Watson:                   0.885
    Prob(Omnibus):                  0.000   Jarque-Bera (JB):            14087.596
    Skew:                           1.082   Prob(JB):                         0.00
    Kurtosis:                       6.420   Cond. No.                     1.94e+18
    ==============================================================================
    
    Notes:
    [1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
    [2] The smallest eigenvalue is 1.82e-26. This might indicate that there are
    strong multicollinearity problems or that the design matrix is singular.


## Check Residuals


```python
yhat = statsModelsCoefs.predict(california_df[featureNames])
resid = y - yhat

plt.figure(figsize=(5,5))
plt.plot(y, yhat, 'ro')
plt.xlabel('Actual')
plt.ylabel('Predicted')
plt.title('Actual vs. Predicted')
plt.show()
```


```python
plt.figure(figsize=(5,5))
plt.hist(resid)
plt.title('Residual Distribution')
plt.show()
```

### Patterns in residuals are not ideal. You'll want these to look like normally distributed white noise (ideally). We might be able to make those residuals behave with feature transformations and other techniques we'll talk about later.

>In certain cases, it may help to log-transform your target variable, which will compress some of the variance.

[Log-linear models](https://en.wikipedia.org/wiki/Log-linear_model)


```python
import statsmodels.api as smf
logMV = np.log(california_df['MedHouseVal'])

statsModelsCoefs = smf.OLS(logMV, california_df[featureNames]).fit()
print(statsModelsCoefs.summary())
```

                                OLS Regression Results                            
    ==============================================================================
    Dep. Variable:            MedHouseVal   R-squared:                       0.614
    Model:                            OLS   Adj. R-squared:                  0.614
    Method:                 Least Squares   F-statistic:                     4109.
    Date:                Wed, 29 Dec 2021   Prob (F-statistic):               0.00
    Time:                        17:41:20   Log-Likelihood:                -7819.3
    No. Observations:               20640   AIC:                         1.566e+04
    Df Residuals:                   20631   BIC:                         1.573e+04
    Df Model:                           8                                         
    Covariance Type:            nonrobust                                         
    ==============================================================================
                     coef    std err          t      P>|t|      [0.025      0.975]
    ------------------------------------------------------------------------------
    const        -12.1230      0.161    -75.393      0.000     -12.438     -11.808
    MedInc         0.1870      0.002     91.271      0.000       0.183       0.191
    HouseAge       0.0023      0.000     10.726      0.000       0.002       0.003
    AveRooms      -0.0333      0.003    -11.594      0.000      -0.039      -0.028
    AveBedrms      0.2344      0.014     17.068      0.000       0.207       0.261
    Population  7.843e-06   2.32e-06      3.384      0.001     3.3e-06    1.24e-05
    AveOccup      -0.0017      0.000     -7.039      0.000      -0.002      -0.001
    Latitude      -0.2845      0.004    -80.984      0.000      -0.291      -0.278
    Longitude     -0.2850      0.004    -77.525      0.000      -0.292      -0.278
    const        -12.1230      0.161    -75.393      0.000     -12.438     -11.808
    ==============================================================================
    Omnibus:                      803.061   Durbin-Watson:                   0.809
    Prob(Omnibus):                  0.000   Jarque-Bera (JB):             2517.584
    Skew:                          -0.036   Prob(JB):                         0.00
    Kurtosis:                       4.709   Cond. No.                     1.94e+18
    ==============================================================================
    
    Notes:
    [1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
    [2] The smallest eigenvalue is 1.82e-26. This might indicate that there are
    strong multicollinearity problems or that the design matrix is singular.



```python
logyhat = statsModelsCoefs.predict(california_df[featureNames])
logresid = logMV - logyhat

plt.figure(figsize=(5,5))
plt.plot(np.exp(logMV), np.exp(logyhat), 'ro', alpha=0.1)
plt.xlabel('Actual')
plt.ylabel('Predicted')
plt.title('Actual vs. Predicted')
plt.show()
```

> We'll go over models later that can catch the tail of this data better.

# What if we had categorical variables?

## Role of [Dummy Variables](https://en.wikipedia.org/wiki/Dummy_variable_(statistics))
A way to incorporate categorical data into modeling, since models require numerical matrices. Essentially acts to change the intercept.

### Simple Example


```python
dummy = pd.DataFrame([[2,4,1],
                     [3,6,1],
                     [4,8,1],
                     [6,12,1],
                     [7,14,1],
                     [2,6,0],
                     [4,10,0],
                     [6,14,0],
                     [7,16,0],
                     [3,8,0]], columns=['hgt', 'wgt', 'label'])


dummy
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>hgt</th>
      <th>wgt</th>
      <th>label</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>2</td>
      <td>4</td>
      <td>1</td>
    </tr>
    <tr>
      <th>1</th>
      <td>3</td>
      <td>6</td>
      <td>1</td>
    </tr>
    <tr>
      <th>2</th>
      <td>4</td>
      <td>8</td>
      <td>1</td>
    </tr>
    <tr>
      <th>3</th>
      <td>6</td>
      <td>12</td>
      <td>1</td>
    </tr>
    <tr>
      <th>4</th>
      <td>7</td>
      <td>14</td>
      <td>1</td>
    </tr>
    <tr>
      <th>5</th>
      <td>2</td>
      <td>6</td>
      <td>0</td>
    </tr>
    <tr>
      <th>6</th>
      <td>4</td>
      <td>10</td>
      <td>0</td>
    </tr>
    <tr>
      <th>7</th>
      <td>6</td>
      <td>14</td>
      <td>0</td>
    </tr>
    <tr>
      <th>8</th>
      <td>7</td>
      <td>16</td>
      <td>0</td>
    </tr>
    <tr>
      <th>9</th>
      <td>3</td>
      <td>8</td>
      <td>0</td>
    </tr>
  </tbody>
</table>
</div>



### Persistent differences between these lines. Track parallel to one another.


```python
class1 = dummy.query('label==1')
class2 = dummy.query('label==0')

plt.figure(figsize=(6,6))
plt.plot(class1['hgt'], class1['wgt'], 'bo')
plt.plot(class2['hgt'], class2['wgt'], 'ro')
plt.legend(['Label = 1', 'Label = 0'])
plt.show()
```

### Let's look at the means for the data by group


```python
dummy.groupby('label')['wgt'].mean().diff()[1]
```




    -2.0



### Compare running a model with and without a dummy variable


```python
from sklearn.linear_model import LinearRegression

Xa = np.array(dummy['hgt']).reshape(-1,1)
Xb = np.array(dummy[['hgt','label']])

y = np.array(dummy['wgt']).reshape(-1,1)

bothOLS = LinearRegression().fit(Xa,y)
yhat_both = bothOLS.predict(Xa)

sepOLS = LinearRegression().fit(Xb, y)
yhat_sep = sepOLS.predict(Xb)

print('No Dummy:\n')
print(f' Intercept: {np.round(bothOLS.intercept_[0],2)}')
print(f' Slope: {np.round(bothOLS.coef_[0][0],2)}\n')

print('w/ Dummy:\n')
print(f' Intercept: {np.round(sepOLS.intercept_[0], 2)}')
print(f' Slope: {np.round(sepOLS.coef_[0][0],2)}')
print(f' Dummy: {np.round(sepOLS.coef_[0][1], 0)}')
```

    No Dummy:
    
     Intercept: 1.0
     Slope: 2.0
    
    w/ Dummy:
    
     Intercept: 2.0
     Slope: 2.0
     Dummy: -2.0


### The dummy captures the mean difference! Otherwise the slope is identical.

If the dummy was present:
$$
y=1.0 + (2)(x_1) + (2)(1) 
$$

If the dummy was not present:
$$
y=2.0 + (2)(x_1) + (2)(0)
$$



# Incorporating Categorical Variables into a Model w/ Pipelines
[Example from "Hands-on Machine Learning with Scikit-Learn, Keras & Tensorflow"](https://github.com/ageron/handson-ml2)

This is an expanded, rawer form of the california housing data that is available in scikit-learn.


```python
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

housing = pd.read_csv('data/housing.csv')
housing.info()
```

    <class 'pandas.core.frame.DataFrame'>
    RangeIndex: 20640 entries, 0 to 20639
    Data columns (total 10 columns):
     #   Column              Non-Null Count  Dtype  
    ---  ------              --------------  -----  
     0   longitude           20640 non-null  float64
     1   latitude            20640 non-null  float64
     2   housing_median_age  20640 non-null  int64  
     3   total_rooms         20640 non-null  int64  
     4   total_bedrooms      20433 non-null  float64
     5   population          20640 non-null  int64  
     6   households          20640 non-null  int64  
     7   median_income       20640 non-null  float64
     8   ocean_proximity     20640 non-null  object 
     9   median_house_value  20640 non-null  int64  
    dtypes: float64(4), int64(5), object(1)
    memory usage: 1.6+ MB


## Goal - Predict Median House Value

#### Things we need to consider:
- Ocean Proximity is categorical. 
- Missing values in Total Bedrooms.  
- Data of significantly different scales.


```python
housing.median_house_value.hist()
plt.title('Distribution of Median Home Values')
plt.show()
```


```python
housing.hist(bins=50, figsize=(8,8))
plt.show()
```

Interesting items of note:
- Wide variances in scales. 
- Median house value truncated at $500,000  
- Outliers are present


```python
housing.groupby('ocean_proximity')['median_house_value'].median().plot.barh()
plt.xlabel('Median of Median Home Value')
plt.ylabel('')
plt.title('Price differences by Ocean Proximity')
plt.show()
```

Inland homes have significantly lower prices than homes closer to the water.


```python
housing.plot(kind='scatter', x='longitude', y='latitude', alpha=0.4, s=housing['population']/100,
            label='population', figsize=(10,7), c='median_house_value', cmap=plt.get_cmap('jet'),
            colorbar=True)
plt.legend('')
plt.show()
```

Higher values are largely clustered around the coast.

[Example from "Hands-on Machine Learning with Scikit-Learn, Keras & TensorFlow](https://github.com/ageron/handson-ml2)

# Missing Values

#### Options:
- Drop rows (can be bad, what if incoming data has missing values?). 
- Drop columns (could be bad, what if there value in that feature?).  
- Fill in the missing values (ding ding).  
- If categorical, might want to add a dummy to indicate it was missing (ding ding).  
- Best strategy will be situationally dependent. This can be treated as a hyperparameter - no perfect answer.

#### Strategies:
- Simple inputers with median, median, mode, random values.  
- Estimate the missing value with another machine learning model (increasing overal complexity).  

>Not all strategies will work for all data types, so you may need to split it up, e.g., one method for the numerical variables and another for the categorical variables.


```python
housing.isna().sum()
```




    longitude               0
    latitude                0
    housing_median_age      0
    total_rooms             0
    total_bedrooms        207
    population              0
    households              0
    median_income           0
    ocean_proximity         0
    median_house_value      0
    dtype: int64



### Filling with median


```python
from sklearn.impute import SimpleImputer

example_imputer = SimpleImputer(strategy='median')
example_imputer.fit_transform(np.array(housing.total_bedrooms).reshape(-1,1))

example_imputer = pd.Series(example_imputer)
example_imputer.isna().sum()
```




    0



### If you wanted to add an indictor for the missing value
> Probably more useful for categorical variables


```python
from sklearn.impute import SimpleImputer

example_imputer = SimpleImputer(strategy='median', add_indicator=True)
example_imputer.fit_transform(np.array(housing.total_bedrooms).reshape(-1,1))
```




    array([[ 129.,    0.],
           [1106.,    0.],
           [ 190.,    0.],
           ...,
           [ 485.,    0.],
           [ 409.,    0.],
           [ 616.,    0.]])



### Note on One-Hot Encoding

From "Hands-on Machine Learning with Scikit-Learn, Keras & TensorFlow":
> If a categorical attribute has a large number of possible categories, then one-hot encoding will result in a large number of input features. This may slow down training and degrade performance.

#### Possible alternatives in that situation:
- Recode to a numerical feature, e.g., distance to ocean.  
- Only use the most frequent $N$ categories.  
- Convert to embeddings.  

#### A risk for us:
May not have any islands in the training data, what would happen if we encountered that in our test/evaluation data?


```python
housing.ocean_proximity.value_counts()
```




    <1H OCEAN     9136
    INLAND        6551
    NEAR OCEAN    2658
    NEAR BAY      2290
    ISLAND           5
    Name: ocean_proximity, dtype: int64



### Islands are really rare

- Adding a dummy for this won't do much - it'll basically be zero.  
- Replace to nearest category?


```python
from sklearn.preprocessing import OneHotEncoder

example_ohe = OneHotEncoder()
example_ohe = example_ohe.fit_transform(np.array(housing['ocean_proximity']).reshape(-1,1))
example_ohe
```




    <20640x5 sparse matrix of type '<class 'numpy.float64'>'
    	with 20640 stored elements in Compressed Sparse Row format>



> Sparse matrix returns lists of coordinates in the matrix with a non-zero. It's a more efficient structure:


```python
print(example_ohe[:5,])
```

      (0, 3)	1.0
      (1, 3)	1.0
      (2, 3)	1.0
      (3, 3)	1.0
      (4, 3)	1.0


>Can be converted back to a dense format:


```python
example_ohe.toarray()[:5,:]
```




    array([[0., 0., 0., 1., 0.],
           [0., 0., 0., 1., 0.],
           [0., 0., 0., 1., 0.],
           [0., 0., 0., 1., 0.],
           [0., 0., 0., 1., 0.]])



### Should we use this for modeling?

>In statistics, multicollinearity (also collinearity) is a phenomenon in which one predictor variable in a multiple regression model can be linearly predicted from the others with a substantial degree of accuracy. In this situation, the coefficient estimates of the multiple regression may change erratically in response to small changes in the model or the data. Multicollinearity does not reduce the predictive power or reliability of the model as a whole, at least within the sample data set; it only affects calculations regarding individual predictors. That is, a multivariate regression model with collinear predictors can indicate how well the entire bundle of predictors predicts the outcome variable, but it may not give valid results about any individual predictor, or about which predictors are redundant with respect to others.
<br><br>[Wikipedia](https://en.wikipedia.org/wiki/Multicollinearity)

### Could argue that this could be represented as an ordinal variable (island > near ocean > near bay > ...)
See [OrdinalEncoder for as example](https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.OrdinalEncoder.html)

You could try this using both methods to see if one is better.

# A Brief Rant on LabelEncoder
Should not be used on your feature set! This is commonly done on Kaggle (incorrectly!). Do not use this for your feature processing!

[LabelEncoder Documentation](https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.LabelEncoder.html?highlight=labelencoder)

# Scaling Numerical Variables
Don't skip this - most models don't perform well when variables are on different scales.


```python
housing.select_dtypes(['float','integer']).describe().round(0)
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>longitude</th>
      <th>latitude</th>
      <th>housing_median_age</th>
      <th>total_rooms</th>
      <th>total_bedrooms</th>
      <th>population</th>
      <th>households</th>
      <th>median_income</th>
      <th>median_house_value</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>count</th>
      <td>20640.0</td>
      <td>20640.0</td>
      <td>20640.0</td>
      <td>20640.0</td>
      <td>20433.0</td>
      <td>20640.0</td>
      <td>20640.0</td>
      <td>20640.0</td>
      <td>20640.0</td>
    </tr>
    <tr>
      <th>mean</th>
      <td>-120.0</td>
      <td>36.0</td>
      <td>29.0</td>
      <td>2636.0</td>
      <td>538.0</td>
      <td>1425.0</td>
      <td>500.0</td>
      <td>4.0</td>
      <td>206856.0</td>
    </tr>
    <tr>
      <th>std</th>
      <td>2.0</td>
      <td>2.0</td>
      <td>13.0</td>
      <td>2182.0</td>
      <td>421.0</td>
      <td>1132.0</td>
      <td>382.0</td>
      <td>2.0</td>
      <td>115396.0</td>
    </tr>
    <tr>
      <th>min</th>
      <td>-124.0</td>
      <td>33.0</td>
      <td>1.0</td>
      <td>2.0</td>
      <td>1.0</td>
      <td>3.0</td>
      <td>1.0</td>
      <td>0.0</td>
      <td>14999.0</td>
    </tr>
    <tr>
      <th>25%</th>
      <td>-122.0</td>
      <td>34.0</td>
      <td>18.0</td>
      <td>1448.0</td>
      <td>296.0</td>
      <td>787.0</td>
      <td>280.0</td>
      <td>3.0</td>
      <td>119600.0</td>
    </tr>
    <tr>
      <th>50%</th>
      <td>-118.0</td>
      <td>34.0</td>
      <td>29.0</td>
      <td>2127.0</td>
      <td>435.0</td>
      <td>1166.0</td>
      <td>409.0</td>
      <td>4.0</td>
      <td>179700.0</td>
    </tr>
    <tr>
      <th>75%</th>
      <td>-118.0</td>
      <td>38.0</td>
      <td>37.0</td>
      <td>3148.0</td>
      <td>647.0</td>
      <td>1725.0</td>
      <td>605.0</td>
      <td>5.0</td>
      <td>264725.0</td>
    </tr>
    <tr>
      <th>max</th>
      <td>-114.0</td>
      <td>42.0</td>
      <td>52.0</td>
      <td>39320.0</td>
      <td>6445.0</td>
      <td>35682.0</td>
      <td>6082.0</td>
      <td>15.0</td>
      <td>500001.0</td>
    </tr>
  </tbody>
</table>
</div>



## Two main methods:
> Only fit these to the training data, no leaking information from the test set!


### Min-max scaling
- Simple  
- Bound between 0 and 1 - a lot of algorithms like that, especially neural networks
- scikit-learn gives you some additional flexibility in terms of the range  
- Very susceptible to outliers

$$
x_{scaled} = \frac{x - x_{min}}{x_{max}-x_{min}}
$$

### Standardization
- Little more involved. 
- More robust to outliers. 
- No specific range boundary. 

$$
x_{scaled} = \frac{x - \hat{x}}{\sigma_{x}}
$$

Since we are doing regression and don't have a scaling boundary requirement and there are probably outliers, we'll use standardization.


```python
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt

example_SS = StandardScaler()
example_SS = example_SS.fit_transform(np.array(housing.total_rooms).reshape(-1,1))

plt.hist(housing.total_rooms, bins=100)
plt.title('Total Rooms', loc='left')
plt.show()

plt.hist(example_SS, bins=100)
plt.title('Total Rooms - Standardized', loc='left')
plt.show()
```

# Training/Test Splits

>Put the test data aside and never look at it again.

- All of the feature transformations and model training should be on the training data.  
- In production, you wouldn't exactly know what the incoming data would look like ahead of time.  
- If you use the test data to inform **any** of the feature transformations or modeling, then you are letting that test data leak into the training data and that will (may) bias your evaluations. This is called **leakage** or **data snooping** - both are not good.

### Simpliest form is splitting your data into two parts:
- Training: will base feature transforms and modeling on this. 
- Test: evaluate the models on this data. 

>There are more robust methods that we'll talk about later. You can think of this simple splitting as a quick and dirty way to evaluate performance, but it isn't a methodology you'd want to use to estimate what your performance is truly likely to be.

### Split off the features and the target variable


```python
y = housing.median_house_value

features = ['housing_median_age', 'total_rooms', 'total_bedrooms', 
            'population', 'households', 'median_income', 'ocean_proximity', 'longitude', 'latitude'
           ]

X = housing[features]
X.info()
```

    <class 'pandas.core.frame.DataFrame'>
    RangeIndex: 20640 entries, 0 to 20639
    Data columns (total 9 columns):
     #   Column              Non-Null Count  Dtype  
    ---  ------              --------------  -----  
     0   housing_median_age  20640 non-null  int64  
     1   total_rooms         20640 non-null  int64  
     2   total_bedrooms      20433 non-null  float64
     3   population          20640 non-null  int64  
     4   households          20640 non-null  int64  
     5   median_income       20640 non-null  float64
     6   ocean_proximity     20640 non-null  object 
     7   longitude           20640 non-null  float64
     8   latitude            20640 non-null  float64
    dtypes: float64(4), int64(4), object(1)
    memory usage: 1.4+ MB


### Split into training and test sets
>80/20 split is pretty standard, but not universal. For very large datasets, I've heard of 99/1 splits.


```python
from sklearn.model_selection import train_test_split

X_training, X_test, y_training, y_test = train_test_split(X, y, test_size=0.20)

print(f'Training samples: {X_training.shape[0]:,}')
print(f'Test samples: {X_test.shape[0]:,}')
```

    Training samples: 16,512
    Test samples: 4,128


>Remember, the test data is only for evaluation.


```python
X_training.info()
```

    <class 'pandas.core.frame.DataFrame'>
    Int64Index: 16512 entries, 9505 to 15571
    Data columns (total 9 columns):
     #   Column              Non-Null Count  Dtype  
    ---  ------              --------------  -----  
     0   housing_median_age  16512 non-null  int64  
     1   total_rooms         16512 non-null  int64  
     2   total_bedrooms      16354 non-null  float64
     3   population          16512 non-null  int64  
     4   households          16512 non-null  int64  
     5   median_income       16512 non-null  float64
     6   ocean_proximity     16512 non-null  object 
     7   longitude           16512 non-null  float64
     8   latitude            16512 non-null  float64
    dtypes: float64(4), int64(4), object(1)
    memory usage: 1.3+ MB


# Pipelines

- Fill missing values. 
- Create dummies for categorical.  
- Standardize numerical variables.  
- Fit the model.  

### Pipelines can be made of collections of pipelines


```python
from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import OneHotEncoder
from sklearn.preprocessing import StandardScaler

cat_vars = ['ocean_proximity']
num_vars = ['housing_median_age', 'total_rooms', 'total_bedrooms', 'population',
            'households', 'median_income', 'longitude', 'latitude']

num_pipeline = Pipeline([('impute_missing', SimpleImputer(strategy='median')),
                           ('standardize_num', StandardScaler())
                        ])

cat_pipeline = Pipeline([('impute_missing_cats', SimpleImputer(strategy='most_frequent')),
                          ('create_dummies_cats', OneHotEncoder(handle_unknown='ignore', drop='first'))])

processing_pipeline = ColumnTransformer(transformers=[('proc_numeric', num_pipeline, num_vars),
                                                      ('create_dummies', cat_pipeline, cat_vars)])

print(processing_pipeline)
```

    ColumnTransformer(transformers=[('proc_numeric',
                                     Pipeline(steps=[('impute_missing',
                                                      SimpleImputer(strategy='median')),
                                                     ('standardize_num',
                                                      StandardScaler())]),
                                     ['housing_median_age', 'total_rooms',
                                      'total_bedrooms', 'population', 'households',
                                      'median_income', 'longitude', 'latitude']),
                                    ('create_dummies',
                                     Pipeline(steps=[('impute_missing_cats',
                                                      SimpleImputer(strategy='most_frequent')),
                                                     ('create_dummies_cats',
                                                      OneHotEncoder(drop='first',
                                                                    handle_unknown='ignore'))]),
                                     ['ocean_proximity'])])



```python
from sklearn.linear_model import LinearRegression

modeling_pipeline = Pipeline([('data_processing', processing_pipeline), ('lm', LinearRegression())])
modeling_pipeline.fit(X_training, y_training)
```




    Pipeline(steps=[('data_processing',
                     ColumnTransformer(transformers=[('proc_numeric',
                                                      Pipeline(steps=[('impute_missing',
                                                                       SimpleImputer(strategy='median')),
                                                                      ('standardize_num',
                                                                       StandardScaler())]),
                                                      ['housing_median_age',
                                                       'total_rooms',
                                                       'total_bedrooms',
                                                       'population', 'households',
                                                       'median_income', 'longitude',
                                                       'latitude']),
                                                     ('create_dummies',
                                                      Pipeline(steps=[('impute_missing_cats',
                                                                       SimpleImputer(strategy='most_frequent')),
                                                                      ('create_dummies_cats',
                                                                       OneHotEncoder(drop='first',
                                                                                     handle_unknown='ignore'))]),
                                                      ['ocean_proximity'])])),
                    ('lm', LinearRegression())])



# Evaluating the model
>Really evaluatign the entire preprocessing process and the model itself.


```python
housing_predictions = modeling_pipeline.predict(X_test)
```

### Get the mean squared error, root mean squared error, and $R^2$


```python
from sklearn.metrics import mean_squared_error
mse = mean_squared_error(y_test, housing_predictions)
mse
```




    4590349642.816633




```python
rmse = np.sqrt(mse)
rmse
```




    67752.11910203719




```python
from sklearn.metrics import r2_score
r2 = r2_score(y_test, housing_predictions)
r2
```




    0.6458661159771042



### Plot the test and predictions


```python
import matplotlib.pyplot as plt

plt.plot(y_test, housing_predictions, 'ro')
plt.xlabel('Actual')
plt.ylabel('Predicted')
plt.show()
```

# Observations
- If this was a perfect model, the residuals would be 0 for all actual/predicted values.  
- Residuals should look like white noise across all values - seeing some patterns.  
- Some information is leaking into the residuals that the model is capturing.  
    - Could be a feature we don't have access to.  
    - Could be noise in the data.  
    - Could be the underlying relationships are linear.  
    - Insert any number of additional explanations.  
    
There may be some overfitting - the training data fits better than the test data. We can explore other models to see if they are able to reduce the overfitting.

# Bias-variance Tradeoff



[Image source](http://scott.fortmann-roe.com/docs/BiasVariance.html)

>At its root, dealing with bias and variance is really about dealing with **over- and under-fitting**. Bias is reduced and variance is increased in relation to model complexity. As more and more parameters are added to a model, the complexity of the model rises and variance becomes our primary concern while bias steadily falls
<br><br>Understanding the Bias-Variance Tradeoff, Fortmann-Roe


From Raschka (paraphrased):
>Variance measures the consistency (or variability) of the model prediction. If we retrain the model on different subsets of the training get and observe difference results, we say it is subject to high variance.

>Bias measures how far off the predictions are from the correct values. This will be error that isn't due to differences in the training datasets.

[Bias and variance from Raschka's Evaluation Lecture Notes](https://sebastianraschka.com/pdf/lecture-notes/stat479fs18/08_eval-intro_notes.pdf)







# Simple Usually Triumphs Over the Complex
It's a balancing act though. You'll need a minimum level of complexity to capture the relationships in the data.



[Image source](http://scott.fortmann-roe.com/docs/BiasVariance.html)

## Potential Options
- Include less features.  
- Shrinkage methods.  
- Data over models.  

# Low Variance, Forward, and Backward Selection

## [Low Variance](https://scikit-learn.org/stable/modules/feature_selection.html)
>You can automatically omit features with zero-variance (i.e., constants).


```python
from sklearn.feature_selection import VarianceThreshold
X = np.array([[0, 0, 1], [0, 1, 0], [1, 0, 0], [0, 1, 1], [0, 1, 0], [0, 1, 1]])
print(X)
```

    [[0 0 1]
     [0 1 0]
     [1 0 0]
     [0 1 1]
     [0 1 0]
     [0 1 1]]



```python
sel = VarianceThreshold(threshold=(.8 * (1 - .8)))
sel.fit_transform(X)
```




    array([[0, 1],
           [1, 0],
           [0, 0],
           [1, 1],
           [1, 0],
           [1, 1]])



## [Forward or Backward Selection](https://scikit-learn.org/stable/modules/generated/sklearn.feature_selection.SequentialFeatureSelector.html#sklearn.feature_selection.SequentialFeatureSelector)

From scikit-learn:
>Forward-SFS is a greedy procedure that iteratively finds the best new feature to add to the set of selected features. Concretely, we initially start with zero feature and find the one feature that maximizes a cross-validated score when an estimator is trained on this single feature. Once that first feature is selected, we repeat the procedure by adding a new feature to the set of selected features. The procedure stops when the desired number of selected features is reached, as determined by the n_features_to_select parameter.

>Backward-SFS follows the same idea but works in the opposite direction: instead of starting with no feature and greedily adding features, we start with all the features and greedily remove features from the set. The direction parameter controls whether forward or backward SFS is used.
<br><br>In general, forward and backward selection do not yield equivalent results. Also, one may be much faster than the other depending on the requested number of selected features: if we have 10 features and ask for 7 selected features, forward selection would need to perform 7 iterations while backward selection would only need to perform 3.

>SFS differs from RFE and SelectFromModel in that it does not require the underlying model to expose a coef_ or feature_importances_ attribute. It may however be slower considering that more models need to be evaluated, compared to the other approaches. For example in backward selection, the iteration going from m features to m - 1 features using k-fold cross-validation requires fitting m * k models, while RFE would require only a single fit, and SelectFromModel always just does a single fit and requires no iterations.

[See a *Comparative Study of Techniques for Large-Scale Feature Selection* for more discussion.](https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.24.4369&rep=rep1&type=pdf)

### Basic Idea
- Start with all the features. 
- Determine the feature provided the least added benefit.    
- Remove the above feature.  
- Continue until reach the desired number of features or hit a threshold.  

### Rational
> Automatically select the most relevant subset of features.  

> Really only necessary if your models don't support regularization.

> Can help with the *Curse of Dimensionality* since it'll generally select a more parsimonious model with dense features.

From Machine Learning with Python, SBS (backward) showed a model with 3 features would have achieved 100% accuracy on the validation data set. See pages 137-139.


[Image source - Raschka's GitHub; Python for Machine Learning 3rd Edition, Figure 4.8](https://github.com/rasbt/python-machine-learning-book-3rd-edition/blob/master/ch04/images/04_08.png)


- I've never used these in practice - there are other ways to guard against overfitting and selecting a more parsimonious model.  
- Can add a lot of computational overhead.  

[See an example from scikit-learn for code example.](https://scikit-learn.org/stable/auto_examples/feature_selection/plot_select_from_model_diabetes.html#sphx-glr-auto-examples-feature-selection-plot-select-from-model-diabetes-py)

# Regularization
>Helps solve overfitting (high variance) - having too many parameters (i.e., too complex).  

>Also helps with multicolinearity (which we saw) and filtering out noise.  

### Types of regularization:
- $L1$ (lasso - least absolute shrinkage and selection operator): penalizing the sum of $\lvert \beta_{j} \rvert$  
- $L2$ (ridge): penalizing the sum of the $\beta_{j}$'s

# Shrinkage with Ridge Regression ($l_2$ Regularization)
## We can modify our loss function to penalize complexity
Shrinking parameter values to penalize for increased complexity (i.e., parameters with meaningful weights).
>Shrinks, but does not eliminate. Coefficients will be non-zero. Will that help when there any many many features?

## Reduces the influence of parameters that carry less weight

Recall $\hat\beta = (X^TX)^{-1}X^{T}y$  
Goal: $\frac{1}{N}\sum{\epsilon_{i}^2}$

#### We can force parameters to shrink towards zero, effectively reducing their influence by adding a penality to the loss function:
$argmin\frac{1}{N}\sum{\epsilon_{i}^2} + \sum{\lambda\beta_{j}^2}$
- What happens with $\lambda=0$?

New closed form: $\hat\beta = (X^TX+\lambda{I})^{-1}X^{T}y$ 

> Ridge has a unique solution, similar to OLS. You can modify the closed form implementation from last class as an independent proof. scikit-learn's implementation uses various optimization methods to solve the loss optimization problem, so it won't be 100% comparable to OLS with $\alpha=0$, it'll be close though.

## Introduces $\lambda$ - our first (official) hyperparameter!
$\lambda$ controls the amount of the penality, it is bounded between 0 and $\infty$.
>When $\lambda=0$, ridge will provide the same coefficients as OLS, since the $\lambda{I}$ will become zero.


```python
from sklearn.linear_model import Ridge

modeling_pipeline_ridge = Pipeline([('data_processing', processing_pipeline), ('ridge', Ridge(alpha=0))])
modeling_pipeline_ridge.fit(X_training, y_training)
modeling_pipeline_ridge['ridge'].coef_
```




    array([ 13124.73824305, -11243.30743505,  29618.1648861 , -44863.09838405,
            31106.85393638,  73800.61689507, -55108.18089897, -56417.73637031,
           -38005.75860477, 100213.5262696 ,  -2754.34307915,   5249.19152698])



> Compared to the below from the original model:


```python
modeling_pipeline['lm'].coef_
```




    array([ 13124.73824305, -11243.30743505,  29618.1648861 , -44863.09838405,
            31106.85393638,  73800.61689507, -55108.18089897, -56417.73637031,
           -38005.75860477, 100213.5262696 ,  -2754.34307915,   5249.19152698])



# Now to evaluate different $\lambda$ values
>We'll need to evaluate a gradient of $\lambda$ values to determine the best one to use.


```python
from collections import defaultdict
alphas = [0, 1, 2, 5, 10, 50]

ridge_results = defaultdict(dict)

for alph in alphas:
    modeling_pipeline_ridge = Pipeline([('data_processing', processing_pipeline), ('ridge', Ridge(alpha=alph))])
    modeling_pipeline_ridge.fit(X_training, y_training)
    ridge_results['coefficients'][alph] = modeling_pipeline_ridge['ridge'].coef_
    ridge_results['training score'][alph] = modeling_pipeline_ridge.score(X_training, y_training)
    ridge_results['test score'][alph] = modeling_pipeline_ridge.score(X_test, y_test)

print('Done')
```

    Done



```python
coefficients_ridge = pd.DataFrame.from_dict(ridge_results['coefficients'])
coefficients_ridge = coefficients_ridge.reset_index()
coefficients_ridge = coefficients_ridge.rename(columns={'index':'coefficient_nbr'})
coefficients_ridge = coefficients_ridge.melt(id_vars='coefficient_nbr', var_name='alpha', value_name='coefficient')

(
coefficients_ridge.pivot_table(index='alpha', columns='coefficient_nbr', values='coefficient')
    .plot(figsize=(8,4),legend=False)
)

plt.title('Ridge Coefficients', loc='left')
plt.xlabel('Alpha (Regularization Amount)')
plt.ylabel('Coefficient')
plt.show()
```

### Changes in $R^2$


```python
ridge_training_r2 = pd.Series(ridge_results['training score'])
ridge_test_r2 = pd.Series(ridge_results['test score'])

ridge_training_r2.plot()
ridge_test_r2.plot()
plt.title('$R^2$ for Ridge Regression')
plt.legend(['Training','Test'])
plt.xlabel('Alpha (Regularization Level)')
plt.ylabel('Percent of Variance Explained')
plt.ylim(0.4, 1)
plt.show()
```

# Another Option is Lasso - Requires Gradient Descent (GD)
OLS and Ridge regression have unique solutions (even though scikit-learn uses optimization). In order to talk about some of the other variants, we need to talk about an optimization technique called gradient descent.

Gradient descent is an optimization technique that allows us to "learn" what the coefficients should be by iteration and continuous improvment. Traditional statisticans don't love this technique since it's not too far away from guessing a bunch of times until you can't guess much better.  



[Image source](https://www.google.com/url?sa=i&url=https%3A%2F%2Fwww.datasciencecentral.com%2Fprofiles%2Fblogs%2Falternatives-to-the-gradient-descent-algorithm&psig=AOvVaw1ki8gWYTrWRy-NKpu7RFgo&ust=1631036623607000&source=images&cd=vfe&ved=0CAsQjRxqFwoTCJiSq6nz6vICFQAAAAAdAAAAABAD)

Gradient descent essentially is looking for the local minimums of loss functions that are differentiable.  

While least-squares does not a closed-form solution, you can also approximate it using gradient descent. Gradient descent will reappear with other algorithms.

GD requires a loss function, which for OLS regression is the sum of squared errors:
$$J(w)=\frac{1}{2}\sum(y^{(i)} - \hat{y}^{(i)})^2$$

This also be used in logistic regression and neural networks.

## Updating weights
All of the weights are set simultaneously.

1. Initialize weights to 0 or small random numbers. 
2. For each training example, $x^{i}$:  
    a. Compute the output value, $\hat{y}$.  
    b. Update the weights.


```python
def gradientDescent(x, y, theta, alpha, m, numIterations):
    thetaHistory = list()
    
    xTrans = x.transpose()
    costList = list()
    
    for i in range(0, numIterations):
        # data x feature weights = y_hat
        hypothesis = np.dot(x, theta)
        # how far we are off
        loss = hypothesis - y 
        # mse
        cost = np.sum(loss ** 2) / (2 * m)
        costList.append(cost)

        # avg gradient per example
        gradient = np.dot(xTrans, loss) / m 

        # update
        theta = theta - alpha * gradient
        thetaHistory.append(theta)
        
    return thetaHistory, costList
```

### Create training data


```python
data_pipeline = Pipeline([('data_processing', processing_pipeline)])
data_pipeline.fit(X_training)

gs_training_data = data_pipeline.fit_transform(X_training)
```

### Run the model


```python
import datetime

start_ts = datetime.datetime.now()
betaHistory, costList = gradientDescent(gs_training_data,y_training,
                                        theta=np.zeros(gs_training_data.shape[1]),
                                        alpha=0.01,
                                        m=gs_training_data.shape[0], numIterations=5000)
                                                                 
end_ts = datetime.datetime.now()

print(f'Completed in {end_ts-start_ts}')

plt.plot(costList)
plt.title(f'Final cost: {costList[-1]:,.2f}', loc='left')
plt.show()
```

## Show changes in $\beta$ throughout the iterations


```python
from collections import defaultdict
thetas = defaultdict(list)

for i in range(len(betaHistory)):
    for j in range(len(betaHistory[i])):
        thetas[j].append(betaHistory[i][j])
        
thetasD = pd.DataFrame.from_dict(thetas)
thetasD.plot(legend=False)
plt.title('Beta Estimates')
plt.ylabel('Coefficient')
plt.xlabel('Iteration')
plt.show()
```

### Predictions


```python
gs_betas = betaHistory[4999]
gs_predictions = np.dot(gs_training_data, gs_betas)

plt.plot(y_training, gs_predictions, 'bo', alpha=0.4)
plt.xlabel('Actual')
plt.ylabel('Predicted')
plt.title('Gradient Descent Regression Fit on Training Data')
plt.show()
```

# Where things can go wrong:
- The learning rate (alpha generally) is really important. If you pick a rate that is too large, you may hop over the minimum and the models will either be very poor or never converge.  
- Other thresholds may take some trial and error.  
- You can get stuck at a local minima and never find the local maxima.  
- You need to know when to stop. It'll keep adjusting coefficients until it reaches an iteration limit or a derivative threshold. Too long and overfitting could occur.

# Lasso Regression ($l_1$ Regularization)
No closed form solution - will need to use optimization techniques regardless.  

$$
J(w)_{lasso}=\sum{(y^{(i)}-\hat{y}^{(i)})^2+\lambda \lvert \lvert w \rvert \rvert_1}
$$

$$
L1:\lambda \lvert \lvert w \rvert \rvert_1 = \lambda \sum{\lvert w_j \rvert}
$$

Could modify the gradient descent with the above, but we'll let sci-kit learn handle it.

## Differences from Ridge



[Introduction to Statistical Learning, Figure 3.11](https://www.statlearning.com)
>Estimation picture for the lasso (left) and ridge regression (right). Shown are contours of the error and constraint functions. The solid blue areas are the constraint regions $\lvert \beta_1 \rvert \leq t$ and $\beta_{1}^2 + \beta_{2}^2 \leq t^2$, respectivity, while the red ellipses are the contours of the least squares error function.

#### Explanation from Raschka (Python Machine Learning 3rd Edition, Chapter 4, pages 129-131):
> We can think of regularization as adding a penalty term to the cost function to encourage smaller weighs; in other words, we penalize large weights. Thus, by increasing the regularization strength via the regularization parameter, we shrink the weights toward zero and decrease the dependence of our model on the training data.

> The shaded regions of represent the regularization "budget" - the combination of the weights cannot exceed those limits. As the regularization term increases, so does the area of that shaded region.

[See Elements of Statistical Learning Section 3.4 for a more thorough discussion.](https://web.stanford.edu/~hastie/ElemStatLearn/)

### TL:DR - Lasso can create sparse models, Ridge cannot.
Ridge will have non-zero estimates for its $\beta$ values, and lasso can result in some $\beta$ values equal to zero (i.e., sparse).

- Lasso should provide better protection to overfitting than Ridge and OLS.  
- Can also be a technique by itself for feature selection.  


```python
from sklearn.linear_model import Lasso
from collections import defaultdict

alphas = [1, 2, 5, 10, 50]

lasso_results = defaultdict(dict)

for alph in alphas:
    modeling_pipeline_lasso = Pipeline([('data_processing', processing_pipeline), ('lasso', Lasso(alpha=alph))])
    modeling_pipeline_lasso.fit(X_training, y_training)
    lasso_results['coefficients'][alph] = modeling_pipeline_lasso['lasso'].coef_
    lasso_results['training score'][alph] = modeling_pipeline_lasso.score(X_training, y_training)
    lasso_results['test score'][alph] = modeling_pipeline_lasso.score(X_test, y_test)
    
coefficients_lasso = pd.DataFrame.from_dict(lasso_results['coefficients'])
coefficients_lasso = coefficients_lasso.reset_index()
coefficients_lasso = coefficients_lasso.rename(columns={'index':'coefficient_nbr'})
coefficients_lasso = coefficients_lasso.melt(id_vars='coefficient_nbr', var_name='alpha', value_name='coefficient')

coefficients_lasso.pivot_table(index='alpha', columns='coefficient_nbr', values='coefficient').plot(figsize=(8,4))
plt.title('Lasso Coefficients', loc='left')
plt.xlabel('Alpha (Regularization Amount)')
plt.ylabel('Coefficient')
plt.legend('')
plt.show()
```


```python
coefficients_ridge.query('coefficient_nbr == 9')
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>coefficient_nbr</th>
      <th>alpha</th>
      <th>coefficient</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>9</th>
      <td>9</td>
      <td>0</td>
      <td>100213.526270</td>
    </tr>
    <tr>
      <th>21</th>
      <td>9</td>
      <td>1</td>
      <td>75158.290077</td>
    </tr>
    <tr>
      <th>33</th>
      <td>9</td>
      <td>2</td>
      <td>60135.335203</td>
    </tr>
    <tr>
      <th>45</th>
      <td>9</td>
      <td>5</td>
      <td>37614.934183</td>
    </tr>
    <tr>
      <th>57</th>
      <td>9</td>
      <td>10</td>
      <td>23187.012620</td>
    </tr>
    <tr>
      <th>69</th>
      <td>9</td>
      <td>50</td>
      <td>5765.391206</td>
    </tr>
  </tbody>
</table>
</div>




```python
coefficients_lasso.query('coefficient_nbr == 9')
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>coefficient_nbr</th>
      <th>alpha</th>
      <th>coefficient</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>9</th>
      <td>9</td>
      <td>1</td>
      <td>94710.893674</td>
    </tr>
    <tr>
      <th>21</th>
      <td>9</td>
      <td>2</td>
      <td>89208.276458</td>
    </tr>
    <tr>
      <th>33</th>
      <td>9</td>
      <td>5</td>
      <td>72700.508688</td>
    </tr>
    <tr>
      <th>45</th>
      <td>9</td>
      <td>10</td>
      <td>45187.627961</td>
    </tr>
    <tr>
      <th>57</th>
      <td>9</td>
      <td>50</td>
      <td>0.000000</td>
    </tr>
  </tbody>
</table>
</div>



> Lasso shrunk to $0$ and ridge did not!

# Reading

[Spurious correlations (High $r^2$ useless relationships)](https://www.tylervigen.com/spurious-correlations)
<br>[Mostly Harmless Econometrics](https://www.mostlyharmlesseconometrics.com)
<br>[Additional examples from statsmodels](https://www.statsmodels.org/stable/examples/notebooks/generated/ols.html)
<br> [Common pitfalls](https://scikit-learn.org/stable/auto_examples/inspection/plot_linear_model_coefficient_interpretation.html)
<br>[Andrew Ng on Regularization](https://www.youtube.com/watch?v=NyG-7nRpsW8&list=PLkDaE6sCZn6Hn0vK8co82zjQtt3T2Nkqc&index=6)
<br>[Great notes from Goker Erdogan](https://github.com/gokererdogan/JaverianaMLCourse/blob/master/Lectures/05.pdf).
<br>[Raschka on Regularization Performance](https://sebastianraschka.com/faq/docs/regularized-logistic-regression-performance.html)
<br>[Raschka's Linear Regression Notes](http://rasbt.github.io/mlxtend/user_guide/regressor/LinearRegression/)

From Raschka's post:
>Regularization does NOT improve the performance on the data set that the algorithm used to learn the model parameters (feature weights). However, it can improve the generalization performance, i.e., the performance on new, unseen data, which is exactly what we want.


```python

```

107cd72711479c23b596084aabd48c0011b2803a
ipynb
Jupyter Notebook
regression/linear-regression-and-regularization.ipynb
appliedecon/data602-lectures
ecc394095270aff9487a68462ad79aa3ddae458b
MIT
2021-09-28T00:27:00.000Z
2022-03-19T17:30:27.000Z
regression/linear-regression-and-regularization.ipynb
appliedecon/data602-lectures
ecc394095270aff9487a68462ad79aa3ddae458b
MIT
regression/linear-regression-and-regularization.ipynb
appliedecon/data602-lectures
ecc394095270aff9487a68462ad79aa3ddae458b
MIT
2021-09-27T16:16:49.000Z
2022-03-23T21:29:25.000Z
Qwen/Qwen-72B
1. YES
2. YES