---
title: "Lec 15 - scikit-learn<br/>Cross-validation"
subtitle: "<br/> Statistical Computing and Computation"
author: "Sta 663 | Spring 2022"
date: "<br/> Dr. Colin Rundel"
output:
  xaringan::moon_reader:
    css: ["slides.css"]
    lib_dir: libs
    nature:
      highlightLines: true
      highlightStyle: solarized-light
      countIncrementalSlides: false
      ratio: "16:9"
---
exclude: true

```{python setup}
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns

import sklearn

plt.rcParams['figure.dpi'] = 200

np.set_printoptions(
  edgeitems=30, linewidth=200,
  precision = 5, suppress=True
  #formatter=dict(float=lambda x: "%.5g" % x)
)

books = pd.read_csv("data/daag_books.csv")


from sklearn.metrics import mean_squared_error
```

```{r r_setup}
knitr::opts_chunk$set(
  fig.align="center",
  cache=FALSE
)
```

```{r hooks}
local({
  hook_err_old <- knitr::knit_hooks$get("error")  # save the old hook
  knitr::knit_hooks$set(error = function(x, options) {
    # now do whatever you want to do with x, and pass
    # the new x to the old hook
    x = sub("## \n## Detailed traceback:\n.*$", "", x)
    x = sub("Error in py_call_impl\\(.*?\\)\\: ", "", x)
    hook_err_old(x, options)
  })
  
  hook_warn_old <- knitr::knit_hooks$get("warning")  # save the old hook
  knitr::knit_hooks$set(warning = function(x, options) {
    x = sub("<string>:1: ", "", x)
    hook_warn_old(x, options)
  })
})
```

---

## Column Transformers

Are a tool for selectively applying transformer(s) to the columns of an array or DataFrame, they function in a way that is similar to a pipeline and similarly have a helper function `make_column_transformer()`.

.small[
```{python}
from sklearn.compose import make_column_transformer
from sklearn.preprocessing import StandardScaler, OneHotEncoder

ct = make_column_transformer(
  (StandardScaler(), ["volume"]),
  (OneHotEncoder(), ["cover"]),
).fit(
  books
)
```

```{python}
ct.get_feature_names_out()
ct.transform(books)
```
]

---

## Keeping or dropping other columns

One addition important argument is `remainder` which determines what happens to not specified columns. The default is `"drop"` which is why `weight` was removed, the alternative is `"passthrough"` which then retains untransformed columns.

.small[
```{python}
ct = make_column_transformer(
  (StandardScaler(), ["volume"]),
  (OneHotEncoder(), ["cover"]),
  remainder = "passthrough"
).fit(
  books
)
```

```{python}
ct.get_feature_names_out()
ct.transform(books)
```
]

---

## Column selection

One lingering issue with the above approach is that we've had to hard code the column names (can also use indexes). Often we want to select columns based on their dtype (e.g. categorical vs numerical) this can be done via pandas or sklearn,

```{python}
from sklearn.compose import make_column_selector
```

.pull-left[ .small[
```{python}
ct = make_column_transformer(
  ( StandardScaler(), 
    make_column_selector(dtype_include=np.number)),
  ( OneHotEncoder(), 
    make_column_selector(dtype_include=[object, bool]))
)

ct.fit_transform(books)
```
] ]


.pull-right[.small[
```{python}
ct = make_column_transformer(
  ( StandardScaler(), 
    books.select_dtypes(include=['number']).columns ),
  ( OneHotEncoder(), 
    books.select_dtypes(include=['object']).columns )
)

ct.fit_transform(books)
```
] ]

.footnote[`make_column_selector()` also supports selecting via `pattern` or excluding via `dtype_exclude`]

---
class: center, middle

## Demo 1 - Putting it together <br/> Interaction model


---
class: center, middle

## Cross validation &<br/>hyper parameter tuning

---

## hw2 ridge regression data

.small[
```{python}
d = pd.read_csv("data/ridge.csv")
d
```
]

--

.small[
```{python}
d = pd.get_dummies(d)
d
```
]

---

## Fitting a ridge regession model

The `linear_model` submodule also contains the `Ridge` model which can be used to fit a ridge regression, usage is identical other than `Ridge()` takes the parameter `alpha` to specify the regularization strength.


```{python}
from sklearn.linear_model import Ridge, LinearRegression

X, y = d.drop(["y"], axis=1), d.y

rg = Ridge(fit_intercept=False, alpha=10).fit(X, y)
lm = LinearRegression(fit_intercept=False).fit(X, y)
```

```{python}
rg.coef_
lm.coef_
```

--

```{python}
mean_squared_error(y, rg.predict(X))
mean_squared_error(y, lm.predict(X))
```

.footnote[Generally for a Ridge (or Lasso) model it is important to scale the features before fitting - in this case this is not necessary as $x\_1,\ldots,x\_4$ all have mean of ~0 and std dev of ~1 ]

---

## Test-Train split

The most basic form of CV is to split the data into a testing and training set, this can be achieved using `train_test_split` from the `model_selection` submodule.

```{python}
from sklearn.model_selection import train_test_split

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

--

.pull-left[
```{python}
X.shape
X_train.shape
X_test.shape
```
]

.pull-right[
```{python}
y.shape
y_train.shape
y_test.shape
```
]

---

.pull-left[ .small[
```{python}
X_train
```
] ]

.pull-right[ .small[
```{python}
y_train
```
] ]

---

## Train vs Test rmse

```{python}
alpha = np.logspace(-2,1, 100)
train_rmse = []
test_rmse = []

for a in alpha:
    rg = Ridge(alpha=a).fit(X_train, y_train)
    
    train_rmse.append( 
      mean_squared_error(y_train, rg.predict(X_train), squared=False) 
    )
    test_rmse.append( 
      mean_squared_error(y_test, rg.predict(X_test), squared=False) 
    )

res = pd.DataFrame(data = {"alpha": alpha, "train_rmse": train_rmse, "test_rmse": test_rmse})
res
```

---

```{python out.width="66%"}
g = sns.relplot(x="alpha", y="value", hue="variable", data = pd.melt(res, id_vars=["alpha"]))
g.set(xscale="log")
```

---

## Best alpha?

.pull-left[
```{python}
min_i = np.argmin(res.train_rmse)
min_i

res.iloc[[min_i],:]
```
]

.pull-right[
```{python}
min_i = np.argmin(res.test_rmse)
min_i

res.iloc[[min_i],:]
```
]

---

## k-fold cross validation

The previous approach was relatively straight forward, but it required a fair bit of book keeping code to implement and we only examined a single test train split. If we would like to perform k-fold cross validation we can use `cross_val_score` from the `model_selection` submodule. 

```{python}
from sklearn.model_selection import cross_val_score

cross_val_score(
  Ridge(alpha=0.59, fit_intercept=False), 
  X, y, 
  cv=5, 
  scoring="neg_root_mean_squared_error"
)
```


.footnote[
🚩🚩🚩 Note that the default k-fold cross validation used here does not shuffle your data which can be massively problematic if your data is ordered 🚩🚩🚩 
]

---

## Controling k-fold behavior

Rather than providing `cv` as an integer, it is better to specify a cross-validation scheme directly (with additional options). Here we will use the `KFold` class from the `model_selection` submodule. 

```{python}
from sklearn.model_selection import KFold

cross_val_score(
  Ridge(alpha=0.59, fit_intercept=False), 
  X, y, 
  cv = KFold(n_splits=5, shuffle=True, random_state=1234), 
  scoring="neg_root_mean_squared_error"
)
```


---

## KFold object

`KFold()` returns a class object which provides the method `split()` which in turn is a generator that returns a tuple with the indexes of the training and testing selects for each fold given a model matrix `X`,

```{python}
ex = pd.DataFrame(data = list(range(10)), columns=["x"])

cv = KFold(5)
for train, test in cv.split(ex):
  print(f'Train: {train} | test: {test}')
```

--

```{python}
cv = KFold(5, shuffle=True, random_state=1234)
for train, test in cv.split(ex):
  print(f'Train: {train} | test: {test}')
```




---

## Train vs Test rmse (again)

```{python}
alpha = np.logspace(-2,1, 30)
test_mean_rmse = []
test_rmse = []
cv = KFold(n_splits=5, shuffle=True, random_state=1234)

for a in alpha:
    rg = Ridge(fit_intercept=False, alpha=a).fit(X_train, y_train)
    
    scores = -1 * cross_val_score(
      rg, X, y, 
      cv = cv, 
      scoring="neg_root_mean_squared_error"
    )
    test_mean_rmse.append(np.mean(scores))
    test_rmse.append(scores)

res = pd.DataFrame(
    data = np.c_[alpha, test_mean_rmse, test_rmse],
    columns = ["alpha", "mean_rmse"] + ["fold" + str(i) for i in range(1,6) ]
)
res
```

---

```{python out.width="66%"}
g = sns.relplot(x="alpha", y="value", hue="variable", data=res.melt(id_vars=["alpha"]), marker="o", kind="line")
g.set(xscale="log")
```

---

## Best alpha? (again)

```{python}
i = res.drop(
  ["alpha"], axis=1
).agg(
  np.argmin
).to_numpy()

i = np.sort(np.unique(i))

res.iloc[ i, : ]
```

---

## Aside - Available metrics 

For most of the cross validation functions we pass in a string instead of a scoring function from the metrics submodule - if you are interested in seeing the names of the possible metrics, these are available via the `sklearn.metrics.SCORERS` dictionary,

```{python}
np.array( sorted(
  sklearn.metrics.SCORERS.keys()
) )
```

---

## Grid Search

We can further reduce the amount of code needed if there is a specific set of parameter values we would like to explore using cross validation. This is done using the `GridSearchCV` function from the `model_selection` submodule.

```{python}
from sklearn.model_selection import GridSearchCV

gs = GridSearchCV(
  Ridge(fit_intercept=False),
  {"alpha": np.logspace(-2, 1, 30)},
  cv = KFold(5, shuffle=True, random_state=1234),
  scoring = "neg_root_mean_squared_error"
).fit(
  X, y
)
```

```{python}
gs.best_index_
gs.best_params_
gs.best_score_
```

---

## `best_estimator_` attribute

If `refit = True` (the default) with `GridSearchCV()` then the `best_estimator_` attribute will be available which gives direct access to the "best" model or pipeline object. This model is constructed by using the parameter(s) that achieved the maximum score and refitting the model to the complete data set.

```{python}
gs.best_estimator_

gs.best_estimator_.coef_

gs.best_estimator_.predict(X)
```



---

## `cv_results_` attribute

Other useful details about the grid search process are stored in the dictionary `cv_results_` attribute which includes things like average test scores, fold level test scores, test ranks, test runtimes, etc.

```{python}
gs.cv_results_.keys()
```

```{python}
gs.cv_results_["param_alpha"]
gs.cv_results_["mean_test_score"]
gs.cv_results_["mean_fit_time"]
```

---

.pull-left[ .small[
```{python ridge_se, eval=FALSE}
alpha = np.array(gs.cv_results_["param_alpha"],dtype="float64")
score = -gs.cv_results_["mean_test_score"]
score_std = gs.cv_results_["std_test_score"]
n_folds = gs.cv.get_n_splits()

plt.figure(layout="constrained")

ax = sns.lineplot(x=alpha, y=score)
ax.set_xscale("log")

plt.fill_between(
  x = alpha,
  y1 = score + 1.96*score_std / np.sqrt(n_folds),
  y2 = score - 1.96*score_std / np.sqrt(n_folds),
  alpha = 0.2
)

plt.show()
```
] ]

.pull-right[
```{python ref.label="ridge_se", echo=FALSE}
```
]

---

## Ridge traceplot

```{python}
alpha = np.logspace(-2,3, 100)
betas = []

for a in alpha:
    rg = Ridge(alpha=a).fit(X, y)
    
    betas.append(rg.coef_)

res = pd.DataFrame(
  data = betas, columns = rg.feature_names_in_
).assign(
  alpha = alpha  
)

res
```

---

```{python out.width="66%"}
g = sns.relplot(
  data = res.melt(id_vars="alpha", value_name="coef values", var_name="feature"),
  x = "alpha", y = "coef values", hue = "feature",
  kind = "line", aspect=2
)
g.set(xscale="log")
```

---

## Exercise 1

Obtain the [diabetes dataset](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.load_diabetes.html) from sklearn using the following code,

```{python}
from sklearn import datasets
X, y = datasets.load_diabetes(return_X_y=True)
```

Our goal is to fit a Lasso model to these data and determine an optimal value of `alpha` using cross validation. Make sure to perform each of the following:

* Verify whether scaling is necessary for these data

* Even if scaling is not necessary, implement a pipeline that integrates `StandardScaler()` and `Lasso()`

* Find the "optimal" value of `alpha` using `GridSearchCV()` and an appropriate metric, how robust does this result appear to be?

* Time permitting, construct a traceplot of coefficients from the lasso models as a function of `alpha`

---

## Dataset details

.small[
```{python}
datasets.load_diabetes()["feature_names"]
print(datasets.load_diabetes()["DESCR"])
```
]