---
title: "scikit-learn<br/>classification"
subtitle: "Lecture 16"
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}
#| include: false
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns

import sklearn

from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import OneHotEncoder, StandardScaler
from sklearn.model_selection import GridSearchCV, KFold, StratifiedKFold, train_test_split
from sklearn.metrics import classification_report

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

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

pd.set_option("display.width", 1000)
pd.set_option("display.max_columns", 10)
pd.set_option("display.precision", 6)
```

```{r}
#| include: false

options(
  width=55
)

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


## OpenIntro - Spam

We will start by looking at a data set on spam emails from the [OpenIntro project](https://www.openintro.org/). A full data dictionary can be found [here](https://www.openintro.org/data/index.php?data=email). To keep things simple this week we will restrict our exploration to including only the following columns: `spam`, `exclaim_mess`, `format`, `num_char`,  `line_breaks`, and `number`.

* `spam` - Indicator for whether the email was spam.
* `exclaim_mess` - The number of exclamation points in the email message.
* `format` - Indicates whether the email was written using HTML (e.g. may have included bolding or active links).
* `num_char` - The number of characters in the email, in thousands.
* `line_breaks` - The number of line breaks in the email (does not count text wrapping).
* `number` - Factor variable saying whether there was no number, a small number (under 1 million), or a big number.

##

:::: {.columns .small}
::: {.column width='50%'}
```{python}
email = pd.read_csv('data/email.csv')[ 
  ['spam', 'exclaim_mess', 'format', 'num_char', 'line_breaks', 'number'] 
]
email
```
:::

::: {.column width='50%'}
Given that `number` is categorical, we will take care of the necessary dummy coding via `pd.get_dummies()`,
```{python}
email_dc = pd.get_dummies(email)
email_dc
```
:::
::::


##

::: {.small}
```{python}
#| out-width: 75%
sns.pairplot(email, hue='spam', corner=True, aspect=1.25)
```
:::

## Model fitting

::: {.small}
```{python}
from sklearn.linear_model import LogisticRegression

y = email_dc.spam
X = email_dc.drop('spam', axis=1)

m = LogisticRegression(fit_intercept = False).fit(X, y)
```
:::

. . .

::: {.small}
```{python}
m.feature_names_in_
m.coef_
```
:::

## A quick comparison

```{r include=FALSE}
d = read.csv("data/email.csv")
d = dplyr::select(d, spam, exclaim_mess, format, num_char, line_breaks, number)
```

:::: {.columns .small}
::: {.column width='50%'}
*R output*

```{r}
glm(spam ~ . - 1, data = d, family=binomial) 
```
:::

::: {.column width='50%'}
*sklearn output*

```{python}
m.feature_names_in_
m.coef_
```
:::
::::


## `sklearn.linear_model.LogisticRegression`

From the documentations,

> This class implements regularized logistic regression using the ‘liblinear’ library, ‘newton-cg’, ‘sag’, ‘saga’ and ‘lbfgs’ solvers. **Note that regularization is applied by default.** It can handle both dense and sparse input. Use C-ordered arrays or CSR matrices containing 64-bit floats for optimal performance; any other input format will be converted (and copied).


## Penalty parameter

🚩🚩🚩 `LogisticRegression()` has a parameter called penalty that applies a `"l1"` (lasso), `"l2"` (ridge), `"elasticnet"` or `None` with `"l2"` being the default. To make matters worse, the regularization is controlled by the parameter `C` which defaults to 1 (not 0) - also `C` is the inverse regularization strength (e.g. different from `alpha` for ridge and lasso models). 🚩🚩🚩

$$
\min_{w, c} \frac{1 - \rho}{2}w^T w + \rho |w|_1 + C \sum_{i=1}^n \log(\exp(- y_i (X_i^T w + c)) + 1),
$$


## Another quick comparison

```{r include=FALSE}
d = read.csv("data/email.csv")
d = dplyr::select(d, spam, exclaim_mess, format, num_char, line_breaks, number)
```

:::: {.columns .small}
::: {.column width='50%'}
*R output*

```{r}
glm(spam ~ . - 1, data = d, family=binomial) 
```
:::

::: {.column width='50%'}
*sklearn output (penalty `None`)*

```{python}
m = LogisticRegression(
  fit_intercept = False, penalty=None
).fit(
  X, y
)
m.feature_names_in_
m.coef_
```
:::
::::


## Solver parameter

It is also possible specify the solver to use when fitting a logistic regression model, to complicate matters somewhat the choice of the algorithm depends on the penalty chosen: 

* `newton-cg` - [`"l2"`, `None`]
* `lbfgs` - [`"l2"`, `None`]
* `liblinear` - [`"l1"`, `"l2"`]
* `sag` - [`"l2"`, `None`]
* `saga` - [`"elasticnet"`, `"l1"`, `"l2"`, `None`]

Also the can be issues with feature scales for some of these solvers:

> **Note:** ‘sag’ and ‘saga’ fast convergence is only guaranteed on features with approximately the same scale. You can preprocess the data with a scaler from sklearn.preprocessing.


## Prediction

Classification models have multiple prediction methods depending on what type of output you would like,

::: {.small}
```{python}
m.predict(X)
```
:::

:::: {.columns .small}
::: {.column width='50%'}
```{python}
m.predict_proba(X)
```
:::

::: {.column width='50%'}
```{python}
m.predict_log_proba(X)
```
:::
::::

## Scoring

Classification models also include a `score()` method which returns the model's *accuracy*,

::: {.small}
```{python}
m.score(X, y)
```
:::

. . .

Other scoring options are available via the [metrics](https://scikit-learn.org/stable/modules/classes.html#classification-metrics) submodule

::: {.small}
```{python}
from sklearn.metrics import accuracy_score, roc_auc_score, f1_score, confusion_matrix
```
:::

:::: {.columns .small}
::: {.column width='50%'}
```{python}
accuracy_score(y, m.predict(X))
roc_auc_score(y, m.predict_proba(X)[:,1])
f1_score(y, m.predict(X))
```
:::

::: {.column width='50%'}
```{python}
confusion_matrix(y, m.predict(X), labels=m.classes_)
```
:::
::::

## Scoring visualizations - confusion matrix

::: {.small}
```{python}
from sklearn.metrics import ConfusionMatrixDisplay
cm = confusion_matrix(y, m.predict(X), labels=m.classes_)

disp = ConfusionMatrixDisplay(cm).plot()
plt.show()
```
:::


## Scoring visualizations - ROC curve {.smaller}

::: {.small}
```{python}
from sklearn.metrics import auc, roc_curve, RocCurveDisplay

fpr, tpr, thresholds = roc_curve(y, m.predict_proba(X)[:,1])
roc_auc = auc(fpr, tpr)
disp = RocCurveDisplay(fpr=fpr, tpr=tpr, roc_auc=roc_auc,
                       estimator_name='Logistic Regression').plot()
plt.show()
```
:::


## Scoring visualizations - Precision Recall {.smaller}

::: {.small}
```{python}
from sklearn.metrics import precision_recall_curve, PrecisionRecallDisplay

precision, recall, _ = precision_recall_curve(y, m.predict_proba(X)[:,1])
disp = PrecisionRecallDisplay(precision=precision, recall=recall).plot()
plt.show()
```
:::


## Another visualization

::: {.small}
```{python}
def confusion_plot(truth, probs, threshold=0.5):
    
    d = pd.DataFrame(
        data = {'spam': y, 'truth': truth, 'probs': probs}
    )
    
    # Create a column called outcome that contains the labeling outcome for the given threshold
    d['outcome'] = 'other'
    d.loc[(d.spam == 1) & (d.probs >= threshold), 'outcome'] = 'true positive'
    d.loc[(d.spam == 0) & (d.probs >= threshold), 'outcome'] = 'false positive'
    d.loc[(d.spam == 1) & (d.probs <  threshold), 'outcome'] = 'false negative'
    d.loc[(d.spam == 0) & (d.probs <  threshold), 'outcome'] = 'true negative'
    
    # Create plot and color according to outcome
    plt.figure(figsize=(12,4))
    plt.xlim((-0.05,1.05))
    sns.stripplot(y='truth', x='probs', hue='outcome', data=d, size=3, alpha=0.5)
    plt.axvline(x=threshold, linestyle='dashed', color='black', alpha=0.5)
    plt.title("threshold = %.2f" % threshold)
    plt.show()
```
:::

##

::: {.small}
```{python out.width='66%'}
#| out-width: 66%
truth = pd.Categorical.from_codes(y, categories = ('not spam','spam'))
probs = m.predict_proba(X)[:,1]
confusion_plot(truth, probs, 0.5)
confusion_plot(truth, probs, 0.25)
```
:::

# Example 1 - DecisionTreeClassifier

# Example 2 - SVC

# MNIST

## MNIST handwritten digits

::: {.small}
```{python}
from sklearn.datasets import load_digits

digits = load_digits(as_frame=True)
```
:::


:::: {.columns .small}
::: {.column width='50%'}
```{python}
X = digits.data
X
```
:::

::: {.column width='50%'} 
```{python}
y = digits.target
y
```
:::
::::



## digit description

::: {.small}
```{python echo=FALSE}
print(digits.DESCR)
```
:::

## Example digits

::: {.small}
```{python}
#| echo: false
fig, axes = plt.subplots(nrows=5, ncols=10, figsize=(10, 6), layout="constrained")
axes2 = [ax for row in axes for ax in row]

for ax, image, label in zip(axes2, digits.images, digits.target):
    ax.set_axis_off()
    img = ax.imshow(image, cmap=plt.cm.gray_r, interpolation="nearest")
    txt = ax.set_title(f"{label}")
    
plt.show()
```
:::


## Doing things properly - train/test split

To properly assess our modeling we will create a training and testing set of these data, only the training data will be used to learn model coefficients or hyperparameters, test data will only be used for final model scoring.

::: {.small}
```{python}
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.33, shuffle=True, random_state=1234
)
```
:::


## Multiclass logistic regression

Fitting a multiclass logistic regression model will involve selecting a value for the `multi_class` parameter, which can be either `multinomial` for multinomial regression or `ovr` for one-vs-rest where `k` binary models are fit.

::: {.small}
```{python}
#| code-line-numbers: |3
mc_log_cv = GridSearchCV(
  LogisticRegression(penalty=None, max_iter = 5000),
  param_grid = {"multi_class": ["multinomial", "ovr"]},
  cv = KFold(10, shuffle=True, random_state=12345)
).fit(
  X_train, y_train
)
```
:::

. . .

::: {.small}
```{python}
mc_log_cv.best_estimator_
mc_log_cv.best_score_
```
:::

. . .

::: {.small}
```{python}
for p, s in  zip(mc_log_cv.cv_results_["params"], mc_log_cv.cv_results_["mean_test_score"]):
  print(p,"Score:",s)
```
:::


## Model coefficients

::: {.small}
```{python}
pd.DataFrame(
  mc_log_cv.best_estimator_.coef_
)

mc_log_cv.best_estimator_.coef_.shape

mc_log_cv.best_estimator_.intercept_
```
:::

## Confusion Matrix

```{python}
#| include: false
np.set_printoptions(
  linewidth=60
)
```


:::: {.columns .small}
::: {.column width='50%'}
**Within sample**
```{python}
accuracy_score(
  y_train, 
  mc_log_cv.best_estimator_.predict(X_train)
)
confusion_matrix(
  y_train, 
  mc_log_cv.best_estimator_.predict(X_train)
)
```
:::

::: {.column width='50%'}
**Out of sample**
```{python}
accuracy_score(
  y_test, 
  mc_log_cv.best_estimator_.predict(X_test)
)
confusion_matrix(
  y_test, 
  mc_log_cv.best_estimator_.predict(X_test),
  labels = digits.target_names
)
```
:::
::::

## Report

::: {.small}
```{python}
print( classification_report(
  y_test, 
  mc_log_cv.best_estimator_.predict(X_test)
) )
```
:::

## ROC & AUC?

These metrics are slightly awkward to use in the case of multiclass problems since they depend on the probability predictions to calculate.

::: {.small}
```{python}
#| error: true
roc_auc_score(
  y_test, mc_log_cv.best_estimator_.predict_proba(X_test)
)
```
:::

. . .

:::: {.columns .small}
::: {.column width='50%'}
```{python error=TRUE}
roc_auc_score(
  y_test, mc_log_cv.best_estimator_.predict_proba(X_test),
  multi_class = "ovr"
)

roc_auc_score(
  y_test, mc_log_cv.best_estimator_.predict_proba(X_test),
  multi_class = "ovo"
)
```
:::

::: {.column width='50%'}
```{python error=TRUE}
roc_auc_score(
  y_test, mc_log_cv.best_estimator_.predict_proba(X_test),
  multi_class = "ovr", average = "weighted"
)

roc_auc_score(
  y_test, mc_log_cv.best_estimator_.predict_proba(X_test),
  multi_class = "ovo", average = "weighted"
)
```
:::
::::

## Prediction

:::: {.columns .small}
::: {.column width='50%'}
```{python}
mc_log_cv.best_estimator_.predict(X_test)
```
:::

::: {.column width='50%'}
```{python}
mc_log_cv.best_estimator_.predict_proba(X_test),
```
:::
::::

## Examining the coefs

::: {.small}
```{python}
#| out-width: 66%
coef_img = mc_log_cv.best_estimator_.coef_.reshape(10,8,8)

fig, axes = plt.subplots(nrows=2, ncols=5, figsize=(10, 5), layout="constrained")
axes2 = [ax for row in axes for ax in row]

for ax, image, label in zip(axes2, coef_img, range(10)):
    ax.set_axis_off()
    img = ax.imshow(image, cmap=plt.cm.gray_r, interpolation="nearest")
    txt = ax.set_title(f"{label}")
    
plt.show()
```
:::

## Example 3 - DecisionTreeClassifier

Using these data we will now fit a `DecisionTreeClassifier` to these data, we will employ `GridSearchCV` to tune some of the parameters (`max_depth` at a minimum) - see the full list [here](https://scikit-learn.org/stable/modules/generated/sklearn.tree.DecisionTreeClassifier.html).

::: {.small}
```{python}
from sklearn.datasets import load_digits
digits = load_digits(as_frame=True)


X, y = digits.data, digits.target
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.33, shuffle=True, random_state=1234
)
```
:::

# Example 4 - GridSearchCV w/ Multiple models <br/> (Trees vs Forests)