---
title: "scikit-learn"
subtitle: "Lecture 14"
author: "Dr. Colin Rundel"
footer: "Sta 663 - Spring 2023"
format:
  revealjs:
    theme: slides.scss
    transition: fade
    slide-number: true
    self-contained: true
execute: 
  echo: true
---

```{python setup}
#| include: false
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns

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

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

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


## scikit-learn

> Scikit-learn is an open source machine learning library that supports supervised and unsupervised learning. It also provides various tools for model fitting, data preprocessing, model selection, model evaluation, and many other utilities.
>
> <br/>
>
> * Simple and efficient tools for predictive data analysis
> * Accessible to everybody, and reusable in various contexts
> * Built on NumPy, SciPy, and matplotlib
> * Open source, commercially usable - BSD license


::: {.aside}
This is one of several other "scikits" (e.g. scikit-image) which are scientific toolboxes built on top of scipy. For a more complete list see [here](https://pypi.org/search/?q=scikit).
:::

## Submodules

The `sklearn` package contains a large number of submodules which are specialized for different tasks / models,

:::: {.columns .tiny}
::: {.column width='50%'}
- `sklearn.base` - Base classes and utility functions
- `sklearn.calibration` - Probability Calibration
- `sklearn.cluster` - Clustering
- `sklearn.compose` - Composite Estimators
- `sklearn.covariance` - Covariance Estimators
- `sklearn.cross_decomposition` - Cross decomposition
- `sklearn.datasets` - Datasets
- `sklearn.decomposition` - Matrix Decomposition
- `sklearn.discriminant_analysis` - Discriminant Analysis
- `sklearn.ensemble` - Ensemble Methods
- `sklearn.exceptions` - Exceptions and warnings
- `sklearn.experimental` - Experimental
- `sklearn.feature_extraction` - Feature Extraction
- `sklearn.feature_selection` - Feature Selection
- `sklearn.gaussian_process` - Gaussian Processes
- `sklearn.impute` - Impute
- `sklearn.inspection` - Inspection
- `sklearn.isotonic` - Isotonic regression
- `sklearn.kernel_approximation` - Kernel Approximation
:::

::: {.column width='50%'}
- `sklearn.kernel_ridge` - Kernel Ridge Regression
- `sklearn.linear_model` - Linear Models
- `sklearn.manifold` - Manifold Learning
- `sklearn.metrics` - Metrics
- `sklearn.mixture` - Gaussian Mixture Models
- `sklearn.model_selection` - Model Selection
- `sklearn.multiclass` - Multiclass classification
- `sklearn.multioutput` - Multioutput regression and classification
- `sklearn.naive_bayes` - Naive Bayes
- `sklearn.neighbors` - Nearest Neighbors
- `sklearn.neural_network` - Neural network models
- `sklearn.pipeline` - Pipeline
- `sklearn.preprocessing` - Preprocessing and Normalization
- `sklearn.random_projection` - Random projection
- `sklearn.semi_supervised` - Semi-Supervised Learning
- `sklearn.svm` - Support Vector Machines
- `sklearn.tree` - Decision Trees
- `sklearn.utils` - Utilities
:::
::::


# Model Fitting

## Sample data

To begin, we will examine a simple data set on the size and weight of a number of books. The goal is to model the weight of a book using some combination of the other features in the data. 

:::: {.columns .small}
::: {.column width='50%'}
The included columns are:

* `volume` - book volumes in cubic centimeters

* `weight` - book weights in grams

* `cover` - a categorical variable with levels `"hb"` hardback, `"pb"` paperback
:::

::: {.column width='50%'}
```{python}
books = pd.read_csv("data/daag_books.csv"); books
```
:::
::::

::: {.aside}
These data come from the `allbacks` data set from the `DAAG` package in R
:::

##

```{python}
#| out-width: 50%
sns.relplot(data=books, x="volume", y="weight", hue="cover")
```

## Linear regression

scikit-learn uses an object oriented system for implementing the various modeling approaches, the class for `LinearRegression` is part of the `linear_model` submodule.

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

. . .

Each modeling class needs to be constructed (potentially with options) and then the resulting object will provide attributes and methods. 

. . .

:::: {.columns .small}
::: {.column width='50%'}
```{python}
lm = LinearRegression()

m = lm.fit(
  X = books[["volume"]],
  y = books.weight
)

m.coef_
m.intercept_
```
:::

::: {.column width='50%' .fragment}
Note `lm` and `m` are labels for the same object,

```{python}
lm.coef_
lm.intercept_
```
:::
::::


## A couple of considerations

When fitting a model, scikit-learn expects `X` to be a 2d array-like object (e.g. a `np.array` or `pd.DataFrame`) and so it will not accept a `pd.Series` or 1d `np.array`.

:::: {.columns .small}
::: {.column width='50%'}
```{python error=TRUE}
lm.fit(
  X = books.volume,
  y = books.weight
)
```
:::

::: {.column width='50%'}
```{python error=TRUE}
lm.fit(
  X = np.array(books.volume),
  y = books.weight
)
```
:::
::::

. . .

:::: {.columns .small}
::: {.column width='50%'}
```{python error=TRUE}
lm.fit(
  X = np.array(books.volume).reshape(-1,1),
  y = books.weight
)
```
:::

::: {.column width='50%'}
```{python error=TRUE}
lm.fit(
  X = books.drop(["weight", "cover"], axis=1),
  y = books.weight
)
```
:::
::::


## Model parameters

Depending on the model being used, there will be a number of parameters that can be configured when creating the model object or via the `set_params()` method.

::: {.small}
```{python}
lm.get_params()
```
:::

. . .

::: {.small}
```{python}
lm.set_params(fit_intercept = False)
```
:::

. . .

::: {.small}
```{python}
lm = lm.fit(X = books[["volume"]], y = books.weight)
lm.intercept_
lm.coef_
```
:::

## Model prediction

Once the model coefficients have been fit, it is possible to predict using the model via the `predict()` method, this method requires a matrix-like `X` as input and in the case of `LinearRegression` returns an array of predicted y values. 

::: {.small}
```{python}
lm.predict(X = books[["volume"]])
```

```{python}
books["weight_lm_pred"] = lm.predict(X = books[["volume"]])
books
```
:::

##

::: {.small}
```{python}
plt.figure()
sns.scatterplot(data=books, x="volume", y="weight", hue="cover")
sns.lineplot(data=books, x="volume", y="weight_lm_pred", color="c")
plt.show()
```
:::


## Residuals?

There is no built in functionality for calculating residuals, so this needs to be done by hand.

::: {.small}
```{python}
#| out-width: 66%
books["resid_lm_pred"] = books["weight"] - books["weight_lm_pred"]

plt.figure(layout="constrained")
ax = sns.scatterplot(data=books, x="volume", y="resid_lm_pred", hue="cover")
ax.axhline(c="k", ls="--", lw=1)
plt.show()
```
:::


## Categorical variables?

Scikit-learn expects that the model matrix be numeric before fitting,

::: {.small}
```{python error=TRUE}
lm = lm.fit(
  X = books[["volume", "cover"]],
  y = books.weight
)
```
:::

. . .

the solution here is to dummy code the categorical variables - this can be done with pandas via `pd.get_dummies()` or with a scikit-learn preprocessor.

::: {.small}
```{python}
pd.get_dummies(books[["volume", "cover"]])
```
:::


## What went wrong?

Do the following results look reasonable? What went wrong?

::: {.small}
```{python}
lm = LinearRegression().fit(
  X = pd.get_dummies(books[["volume", "cover"]]),
  y = books.weight
)

lm.intercept_
lm.coef_
```
:::



## Quick comparison with R

::: {.small}
```{r}
d = read.csv('data/daag_books.csv')
d['cover_hb'] = ifelse(d$cover == "hb", 1, 0)
d['cover_pb'] = ifelse(d$cover == "pb", 1, 0)
lm = lm(weight~volume+cover_hb+cover_pb, data=d)
summary(lm)
```
:::


## Avoiding co-linearity

:::: {.columns .small}
::: {.column width='50%'}
```{python}
lm = LinearRegression(
  fit_intercept = False
).fit(
  X = pd.get_dummies(books[["volume", "cover"]]),
  y = books.weight
)

lm.intercept_
lm.coef_
lm.feature_names_in_
```
:::

::: {.column width='50%'}
```{python}
lm = LinearRegression(
).fit(
  X = pd.get_dummies(
    books[["volume", "cover"]], 
    drop_first=True
  ),
  y = books.weight
)

lm.intercept_
lm.coef_
lm.feature_names_in_
```
:::
::::


# Preprocessors

## Preprocessors

These are a set of transformer classes present in the `sklearn.preprocessing` submodule that are designed to help with the preparation of raw feature data into quantities more suitable for downstream modeling tools.

Like the modeling classes, they have an object oriented design that shares a common interface (methods and attributes) for bringing in data, transforming it, and returning it.


## OneHotEncoder

For dummy coding we can use the `OneHotEncoder` preprocessor, the default is to use one hot encoding but standard dummy coding can be achieved via the `drop` parameter.

::: {.small}
```{python}
from sklearn.preprocessing import OneHotEncoder
```
:::

:::: {.columns .small}
::: {.column width='50%'}
```{python}
enc = OneHotEncoder(sparse_output=False)
enc.fit(X = books[["cover"]])
enc.transform(X = books[["cover"]])
```
:::

::: {.column width='50%' .fragment}
```{python}
enc = OneHotEncoder(sparse_output=False, drop="first")
enc.fit_transform(X = books[["cover"]])
```
:::
::::


## Other useful bits

```{python, include=FALSE}
enc = OneHotEncoder(sparse_output=False)
enc.fit(X = books[["cover"]])
```

:::: {.columns .small}
::: {.column width='50%'}
```{python}
enc.get_feature_names_out()
f = enc.transform(X = books[["cover"]])
f
```
:::

::: {.column width='50%' .fragment}
```{python}
enc.inverse_transform(f)
```
:::
::::



## A cautionary note

Unlike `pd.get_dummies()` it is not safe to use `OneHotEncoder` with both numerical and categorical features, as the former will also be transformed.

::: {.small}
```{python}
enc = OneHotEncoder(sparse_output=False)
X = enc.fit_transform(X = books[["volume", "cover"]])
pd.DataFrame(data=X, columns = enc.get_feature_names_out())
```
:::


## Putting it together

:::: {.columns .small}
::: {.column width='50%'}
```{python}
cover = OneHotEncoder(
  sparse_output=False
).fit_transform(
  books[["cover"]]
)
X = np.c_[books.volume, cover]

lm2 = LinearRegression(
  fit_intercept=False
).fit(
  X = X,
  y = books.weight
)

lm2.coef_
```
:::

::: {.column width='50%' .fragment}
```{python}
books["weight_lm2_pred"] = lm2.predict(X=X)
books.drop(["weight_lm_pred", "resid_lm_pred"], axis=1)
```
:::
::::

::: {.aside}
We'll see a more elegant way of doing this in the near future
:::

## Model fit

```{python}
#| echo: false
plt.figure()
sns.scatterplot(data=books, x="volume", y="weight", hue="cover")
sns.lineplot(data=books, x="volume", y="weight_lm2_pred", hue="cover")
plt.show()
```

## Model residuals

```{python}
#| echo: false
books["resid_lm2_pred"] = books["weight"] - books["weight_lm2_pred"]

plt.figure(layout="constrained")
ax = sns.scatterplot(data=books, x="volume", y="resid_lm2_pred", hue="cover")
ax.axhline(c="k", ls="--", lw=1)
plt.show()
```


## Model performance

Scikit-learn comes with a number of builtin functions for measuring model performance in the `sklearn.metrics` submodule - these are generally just functions that take the vectors `y_true` and `y_pred` and return a scalar score.

::: {.small}
```{python}
from sklearn.metrics import mean_squared_error, r2_score
```
:::

:::: {.columns .small}
::: {.column width='50%'}
```{python}
r2_score(books.weight, books.weight_lm_pred)
mean_squared_error(
  books.weight, books.weight_lm_pred
)
mean_squared_error(
  books.weight, books.weight_lm_pred,
  squared=False
)
```
:::

::: {.column width='50%'}
```{python}
r2_score(books.weight, books.weight_lm2_pred)
mean_squared_error(
  books.weight, books.weight_lm2_pred
) 
mean_squared_error(
  books.weight, books.weight_lm2_pred, 
  squared=False
)
```
:::
::::

::: {.aside}
See [API Docs](https://scikit-learn.org/stable/modules/classes.html#module-sklearn.metrics) for a list of available metrics
:::


## Exercise 1

Create and fit a model for the `books` data that includes an interaction effect between `volume` and `cover`. 

You will need to do this manually with `pd.getdummies()` and some additional data munging.

The data can be read into pandas with,
```{python}
#| eval: false
books = pd.read_csv(
  "https://sta663-sp23.github.io/slides/data/daag_books.csv"
)
```



# Other transformers

## Polynomial regression

We will now look at another flavor of regression model, that involves preprocessing and a hyperparameter - namely polynomial regression.

```{python, out.width="40%"}
df = pd.read_csv("data/gp.csv")
sns.relplot(data=df, x="x", y="y")
```


## By hand

It is certainly possible to construct the necessary model matrix by hand (or even use a function to automate the process), but this is less then desirable generally - particularly if we want to do anything fancy (e.g. cross validation)

:::: {.columns .small}
::: {.column width='50%'}
```{python}
X = np.c_[
    np.ones(df.shape[0]),
    df.x,
    df.x**2,
    df.x**3
]

plm = LinearRegression(
  fit_intercept = False
).fit(
  X=X, y=df.y
)

plm.coef_
```
:::

::: {.column width='50%' .fragment}
```{python}
df["y_pred"] = plm.predict(X=X)

plt.figure(layout="constrained")
sns.scatterplot(data=df, x="x", y="y")
sns.lineplot(data=df, x="x", y="y_pred", color="k")
plt.show()
```
:::
::::


## PolynomialFeatures

This is another transformer class from `sklearn.preprocessing` that simplifies the process of constructing polynormial features for your model matrix. Usage is similar to that of `OneHotEncoder`.

::: {.small}
```{python}
from sklearn.preprocessing import PolynomialFeatures
X = np.array(range(6)).reshape(-1,1)
```
:::

:::: {.columns .small}
::: {.column width='50%'}
```{python}
pf = PolynomialFeatures(degree=3)
pf.fit(X)
pf.transform(X)
pf.get_feature_names_out()
```
:::

::: {.column width='50%' .fragment}
```{python}
pf = PolynomialFeatures(
  degree=2, include_bias=False
)
pf.fit_transform(X)
pf.get_feature_names_out()
```
:::
::::


## Interactions

If the feature matrix `X` has more than one column then `PolynomialFeatures` transformer will include interaction terms with total degree up to `degree`.

:::: {.columns .small}
::: {.column width='50%'}
```{python}
X.reshape(-1, 2)

pf = PolynomialFeatures(
  degree=3, include_bias=False
)
pf.fit_transform(
  X.reshape(-1, 2)
)
pf.get_feature_names_out()
```
:::

::: {.column width='50%' .fragment}
```{python}
X.reshape(-1, 3)

pf = PolynomialFeatures(
  degree=2, include_bias=False
)
pf.fit_transform(
  X.reshape(-1, 3)
)
pf.get_feature_names_out()
```
:::
::::


## Modeling with PolynomialFeatures

:::: {.columns .small}
::: {.column width='50%'}
```{python}
def poly_model(X, y, degree):
  X  = PolynomialFeatures(
    degree=degree, include_bias=False
  ).fit_transform(
    X=X
  )
  y_pred = LinearRegression(
  ).fit(
    X=X, y=y
  ).predict(
    X
  )
  return mean_squared_error(y, y_pred, squared=False)
```

```{python}
poly_model(X = df[["x"]], y = df.y, degree = 2)
poly_model(X = df[["x"]], y = df.y, degree = 3)
```
:::

::: {.column width='50%' .fragment}
```{python}
#| out-width: 66%
degrees = range(1,10)
rmses = [
  poly_model(X=df[["x"]], y=df.y, degree=d) 
  for d in degrees
]
sns.relplot(x=degrees, y=rmses)
```
:::
::::

##

```{python}
#| echo: false
res = df.copy().drop("y_pred", axis=1)
for d in range(1,10):
  X  = PolynomialFeatures(
    degree=d, include_bias=False
  ).fit_transform(
    X=res[["x"]]
  )
  res[str(d)] = LinearRegression().fit(X=X, y=res.y).predict(X)

g = sns.relplot(
  data = res.melt(id_vars=["x","y"], var_name="degree"),
  x = "x", y="value", col = "degree",
  col_wrap=3, kind="line", color="k"
)

[ ax.scatter(res.x, res.y, alpha=0.3)  for ax in g.axes ]
```

