ml_engineering/4.3_evaluation_pt.2/evaluation_pt.2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 4.3: Evaluation Pt.2\n",
    "\n",
    "[**Lecture Slides**](https://docs.google.com/presentation/d/1HXO3G2Z9fAHOV7itzJ7jmUbggEDNZrg99rB7P6d_1OE/edit?usp=sharing)\n",
    "\n",
    "This lecture, we are going to develop machine learning models for a small imbalanced dataset.\n",
    "\n",
    "**Learning goals:**\n",
    "- stratified sampling of train, val, & test sets\n",
    "- compare accuracy and F-score metrics\n",
    "- apply cross-validation to estimate validation metrics\n",
    "- combine cross-validation and stratification\n",
    "- estimate confidence intervals with the bootstrap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Introduction\n",
    "\n",
    "The [International Sommelier Association](https://asi.info/) want to change the world of wines with new technologies. 🍷 They are interested to know if the quality assessment of wines can be automated: sommeliers are expensive, slow, and rare... if only robots could analyse new harvests! Their goal is to find wines with high potential.\n",
    "\n",
    "We are asked to investigate potential machine learning solutions. We are provided a dataset with 6500 wines, 11 features, a binary label `tasty`, and our task is to find the `tasty` wines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "df = pd.read_csv('wine_quality.csv')\n",
    "df.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(df['tasty'] == 1).sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Only 3% of wines are considered tasty, 198 in total, thus this is an imbalanced learning regime. ⚖️ In some ways, this is also a low-data regime, because there are few positives, and we are interested in finding the positive class."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Stratification\n",
    "\n",
    "We have already defined the task, the data, and the project seems feasible. So we are ready to randomly sample test and validation sets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "train_val_df, test_df = train_test_split(df, test_size=0.2, random_state=1337)\n",
    "train_df, val_df = train_test_split(train_val_df, test_size=0.2, random_state=1337)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since there are only 200 good quality wines, we'd like to see how these were distributed between the three splits:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def tasty_ratio(df):\n",
    "    n_tasty = df.loc[df['tasty'] == 1, 'tasty'].count()\n",
    "    return n_tasty / len(df)\n",
    "\n",
    "print(f'train set tasty ratio: {tasty_ratio(train_df):.3f}')\n",
    "print(f'val set tasty ratio: {tasty_ratio(val_df):.3f}')\n",
    "print(f'test set tasty ratio: {tasty_ratio(test_df):.3f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our validation set has almost 50% more good quality wines than our train set! 🙀 This goes against Andrew Ng's famous words:\n",
    "\n",
    "> Choose dev & test sets to reflect data you expect to get in the future, and want to do well on.\n",
    "\n",
    "If the validation set has a unnatural positive class ratio, then its evaluation metrics will be _biased_. We can fix this through _stratification_. sklearn knows this is a common usecase, and incorporates the `stratify` argument in the `train_test_split()` function. One can pass the variable defining the groups which are to be stratified: in our case, the labels `tasty`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "train_val_df, test_df = train_test_split(df, test_size=0.2, random_state=1337, stratify=df['tasty'])\n",
    "train_df, val_df = train_test_split(train_val_df, test_size=0.2, random_state=1337, stratify=train_val_df['tasty'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'train set tasty ratio: {tasty_ratio(train_df)}')\n",
    "print(f'val set tasty ratio: {tasty_ratio(val_df)}')\n",
    "print(f'test set tasty ratio: {tasty_ratio(test_df)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now our train, val, & test sets have the same class _prevalence_ , and this has avoided some of the bias in our evaluations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Data Munging\n",
    "\n",
    "Let's start by putting aside our test set, and saving it to disk:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_df.to_csv('wine_quality_test.csv', index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then convert the train and val sets to features:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def to_features(df):\n",
    "    X = df.drop(['tasty'], axis=1).values\n",
    "    y = df['tasty'].values\n",
    "    return X, y\n",
    "\n",
    "X_train, y_train = to_features(train_df)\n",
    "X_val, y_val = to_features(val_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And let's remember to scale our features. Feature scaling should **not** be fit on evaluation sets. This will leak information from the val and/or test set to the trained model, and potentially cause overfitting.\n",
    "\n",
    "i.e use the mean and standard deviation of the _training set_ , and apply the transformation to all sets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "scaler = StandardScaler().fit(X_train)\n",
    "X_train = scaler.transform(X_train)\n",
    "X_val = scaler.transform(X_val)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Training\n",
    "\n",
    "We'll start our experiments with three common models as baselines: a random forest classifier, a support vector machine with radial basis function kernel, and a logistic regression model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.svm import SVC\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "\n",
    "rf_clf = RandomForestClassifier(random_state=0).fit(X_train, y_train)\n",
    "svm_clf = SVC(kernel='rbf', random_state=0).fit(X_train, y_train)\n",
    "lr_clf = LogisticRegression(random_state=0).fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Evaluation\n",
    "\n",
    "### 5.1 Accuracy\n",
    "\n",
    "Once these are trained, we can call the `.score()` method of the model api to calculate the validation accuracy, like last lecture:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clfs = [rf_clf, svm_clf, lr_clf]\n",
    "\n",
    "for clf in clfs:\n",
    "    accuracy = clf.score(X_val, y_val)\n",
    "    print(f'classifier: {type(clf).__name__}, validation accuracy: {accuracy}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All these baselines look excellent! 🤩\n",
    "\n",
    "### 5.2 F-score\n",
    "\n",
    "We should be cautious however: we have identified this dataset as _imbalanced_ , and we have learned in the lecture slides that accuracy is a bad metric for the imbalanced learning regime. 🤨 Let's try another metric more apt for the task: F1 score. We'll also record precision and recall to give more insights into the model performance, but F1 score is our choice of single number metric which defines model quality.\n",
    "\n",
    "Just like the `accuracy_score()` function in the `metrics` module, sklearn provides `f1_score`, `precision_score`, and `recall_score` functions. We can wrap these in an `evaluate()` function to avoid explicitly calling all of them for every measurement. We'll also use a [`namedtuple`](https://docs.python.org/2/library/collections.html#collections.namedtuple) to hold the three metric results.\n",
    "\n",
    "Let's evaluate the logitic regression which had a promising 97% accuracy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import precision_score, recall_score, f1_score\n",
    "from collections import namedtuple\n",
    "\n",
    "FPR = namedtuple('FPR', ['f1', 'precision', 'recall'])\n",
    "\n",
    "def evaluate(clf, X_val, y_val):\n",
    "    y_pred = clf.predict(X_val)\n",
    "    p = precision_score(y_val, y_pred, zero_division=0)\n",
    "    r = recall_score(y_val, y_pred)\n",
    "    f1 = f1_score(y_val, y_pred)\n",
    "    return FPR(f1=f1, precision=p, recall=r)\n",
    "\n",
    "print(f'logistic regression validation metrics: {evaluate(lr_clf, X_val, y_val)}')\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "🤯 The precision and recall are 0, which means that the logistic regression model found _none_ of the `tasty` wines! The 97% accuracy was just a consequence of the fact that 97% of the wines are not `tasty`. Logistic regression literally ended up being our \"terrible rule-based model\" example from the lecture slides, and predicted _all_ wines as negatives. This can happen when the training set is heavily imbalanced, and there are [techniques](https://github.com/scikit-learn-contrib/imbalanced-learn) to deal with this which we will mention later in the course.\n",
    "\n",
    "Let's check the other two baselines:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'SVM validation metrics: {evaluate(svm_clf, X_val, y_val)}')\n",
    "print(f'random forest validation metrics: {evaluate(rf_clf, X_val, y_val)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our SVM is also a useless model, but the random forest has a non-zero F1-score.\n",
    "\n",
    "🧠 Define precision and recall. What do they represent here? What type of error do we expect this model to make most of the time?\n",
    "\n",
    "Our suspicion is that the SVM's regularization hyperparameters are causing it to fail. Let's search for a better value of `C` by measuring our single-number validation metric like last lecture: 🕵️"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c_values = [0.1, 1, 10, 100, 1000, 10000]\n",
    "svms = [SVC(kernel='rbf', C=c, random_state=0) for c in c_values]\n",
    "svms = [svm.fit(X_train, y_train) for svm in svms]\n",
    "f1_scores = [evaluate(svm, X_val, y_val).f1 for svm in svms]\n",
    "results = zip(c_values, f1_scores)\n",
    "\n",
    "for c, f1 in results:\n",
    "    print(f'SVM, C: {c}, f1-score: {f1:.3f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`C=1000` improves the performance of the SVM model, but it's still not as good as our random forest model. Let's try to improve that instead:\n",
    "\n",
    "\n",
    "💪💪 Find better hyperparameters for our random forest model. \n",
    "* this is an open-ended exercise, check the [documentation](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestClassifier.html), revisit lecture 3.10 for a fresher on random forest optimization, and be creative! \n",
    "* your results should be f1-scores comparable to the value obtained a few cells above. Remember that even bad scores are important information when tuning hyperparameters, and should be reported too\n",
    "* based on your experiments, conclude which random forest hyperparameters best solve this task"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# INSERT YOUR CODE HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "💪Use this [metric](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.classification_report.html) to create a classification report for our random forest classifier's predictions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# INSERT YOUR CODE HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "🧠🧠 What do the two rows `0` & `1` for precision, recall, and f1-score mean? what's the macro vs weighted accuracy averages?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "💪 Calculate the F2-score of the random forest classifier.\n",
    "* F2 is the same as F$\\beta$ with beta=2\n",
    "* Use the [official documentation](https://scikit-learn.org/stable/modules/classes.html#module-sklearn.metrics), "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# INSERT YOUR CODE HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.3 Crossvalidation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_val.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our validation set contains only 32 positive examples. A metric estimated on this set will probably have high _variance_. However, in the lecture slides we learned about a technique which leverages our dataset more efficiently: cross-validation. Cross-validation allows to validate on the _entire_ dataset, by splitting the examples into folds.\n",
    "\n",
    "\n",
    "For this, we need to return to our `train_val_df`, which we aim to cut into K folds of independent train & val set pairs. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train_val, y_train_val = to_features(train_val_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `sklearn.model_selection` module contains a variety of splitter classes, like [`KFold`](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.KFold.html#sklearn.model_selection.KFold). They can be used to iterate through the folds, as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import KFold\n",
    "\n",
    "kf = KFold(n_splits=5, random_state=1337, shuffle=True)\n",
    "\n",
    "for train_indices, val_indices in kf.split(X_train_val, y_train_val):\n",
    "    X_train, y_train = X_train_val[train_indices], y_train_val[train_indices]\n",
    "    X_val, y_val = X_train_val[val_indices], y_train_val[val_indices]\n",
    "    print(f'train set size: {len(X_train)}, val set size: {len(X_val)}, val positives count: {y_val.sum()}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`KFold` returns the rotating indices corresponding to each train/val split for each fold. These indices can then be used to access the correct examples from our feature matrix & label vector. As expected, the fold sizes are consistent, but once again we observe random variation in the class balance of our splits. To remedy this effect, we can use the `StratifiedKFold` splitter instead:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import StratifiedKFold\n",
    "\n",
    "skf = StratifiedKFold(n_splits=5, random_state=1337, shuffle=True)\n",
    "\n",
    "for train_indices, val_indices in skf.split(X_train_val, y_train_val):\n",
    "    X_train, y_train = X_train_val[train_indices], y_train_val[train_indices]\n",
    "    X_val, y_val = X_train_val[val_indices], y_train_val[val_indices]\n",
    "    print(f'train set size: {len(X_train)}, val set size: {len(X_val)}, val positives count: {y_val.sum()}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our folds are now class balanced! 🤹\n",
    "\n",
    "Using splitters gives a lot of control over the splitting loop. However, cross-validation is often used as part of model evaluation. sklearn provides the [`cross_val_score()`](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.cross_val_score.html) function for this exact purpose. It takes care of _almost_ for us, including splitting the folds, training the model on the train sets, evaluating them on the validation sets. 😎 However, it doesn't take care of feature scaling... \n",
    "\n",
    "We could fit a `StandardScaler` on the whole `X_train_val` feature matrix. But this would break the train/validation set independence, and risk _overfitting_. 😰 We need the feature scaling to be fitted on _each_ train split, and applied to each matching val split. To do so, we can add `StandardScaler` to our `RandomForestClassifier` with a [`Pipeline`](https://scikit-learn.org/stable/modules/compose.html):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.pipeline import make_pipeline\n",
    "\n",
    "pipe = make_pipeline(StandardScaler(), RandomForestClassifier(random_state=0))\n",
    "pipe"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is now an estimator which follows the model api. i.e, the standard scaler will be `.fit()` only when the pipe's `.fit()` method is called (more details [here](https://scikit-learn.org/stable/modules/compose.html)). This prevents our `StandardScaler` from ever being fit on the validation splits.\n",
    "\n",
    "All that is left to do is to provide our pipeline, dataset, and desired number of folds to the `cross_val_score()` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import cross_val_score\n",
    "\n",
    "scores = cross_val_score(pipe, X_train_val, y_train_val, cv=5)\n",
    "scores.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "🧠 List all the steps that happen inside the `cross_val_score()` function above.\n",
    "\n",
    "🤔 sklearn does take care of everything... but also defaults back to accuracy as metric. To measure f1-score, we can use the [`cross_validate()`](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.cross_validate.html) function, similar to `cross_val_score()` but with more control, including specifying the `metrics`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import cross_validate\n",
    "\n",
    "pipe = make_pipeline(StandardScaler(), RandomForestClassifier(random_state=0))\n",
    "metrics= ['precision', 'recall', 'f1']\n",
    "scores = cross_validate(pipe, X_train_val, y_train_val, scoring=metrics, cv=5)\n",
    "scores"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`cross_validate()` returns a lot of information, including the training & prediction times. We can average the fold results, to get the final cross-validation f1-score estimate:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scores['test_f1'].mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is significantly different than our hold-out set estimate! The cross-validation estimate is better because it is applied on more data, and benefits of the whole diversity of examples in the dataset.\n",
    "\n",
    "We suspect the variance of this estimate could be further reduced by using stratification inside the cross-validation splits... except it already is! The [documentation](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.cross_validate.html) mentions that [`StratifiedKFold`](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.StratifiedKFold.html#sklearn.model_selection.StratifiedKFold) is used by default.\n",
    "\n",
    "We want to repeat our SVM hyperparameter optimization with _cross-validation_ instead of a hold-out validation set. This should improve the quality of the metrics, and therefore of the hyperparameter search. We write an `evaluate_crossval()` function, and then once again loop through values of `C`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings \n",
    "\n",
    "def evaluate_crossval(clf, X_train_val, y_train_val):\n",
    "    metrics= ['precision', 'recall', 'f1']\n",
    "    # noisy warnings about 0 precision\n",
    "    with warnings.catch_warnings():\n",
    "        warnings.simplefilter(\"ignore\")\n",
    "        scores = cross_validate(clf, X_train_val, y_train_val, scoring=metrics, cv=5)    \n",
    "    return FPR(f1=scores['test_f1'].mean(), precision=scores['test_precision'].mean(), recall=scores['test_recall'].mean())\n",
    "\n",
    "\n",
    "c_values = [0.1, 1, 10, 100, 1000, 10000]\n",
    "svms = [SVC(kernel='rbf', C=c, random_state=0) for c in c_values]\n",
    "pipes = [make_pipeline(StandardScaler(), svm) for svm in svms]\n",
    "f1_scores = [evaluate_crossval(pipe, X_train_val, y_train_val).f1 for pipe in pipes]\n",
    "results = zip(c_values, f1_scores)\n",
    "\n",
    "for c, f1 in results:\n",
    "    print(f'SVM, C: {c}, f1-score: {f1:.3f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The optimal `C` value is the same, but the f1-scores are significantly lower. \n",
    "\n",
    "🧠 Which do you trust more: these f1-scores, or the ones from our first SVM hyperparameter search? why?\n",
    "\n",
    "🧠 Cite one disadvantage with using cross-validation instead of a hold-out validation set."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.4 Confidence Intervals\n",
    "\n",
    "After experimenting with models and hyperparameters, we have chosen our best wine tastiness detector candidate: a random forest 🌲. We'd like to evaluate it on the _test set_ , to measure its expected production performance.\n",
    "\n",
    "Since we're not using the validation set for evaluation here, we can combine it with the train set to optimize our model.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pipe = make_pipeline(StandardScaler(), RandomForestClassifier(random_state=0))\n",
    "pipe = pipe.fit(X_train_val, y_train_val)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that our model is trained on all non-test data, we can load the test set, and evaluate the random forest:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "test_df = pd.read_csv('wine_quality_test.csv')\n",
    "X_test, y_test = to_features(test_df)\n",
    "evaluate(pipe, X_test, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ℹ️ Note that we did not need to explicitly feature scale `X_test`. This is because we're using a `Pipeline`, so the `StandardScaler` is automatically applied by sklearn. \n",
    "\n",
    "The test score is higher than the cross-validation score. This is a sign that this f1-score estimate has high _variance_. 😕 This agrees with our previous results (which varied quite a lot), and our dataset size. This could be problematic for our tasty wine detector. What if the f1-score was really closer to 0.45? Maybe the International Sommelier Association has a minimum quality requirement of 0.5!\n",
    "\n",
    "Test set evaluations will always be _approximations_ of the population metric. However, we can get an idea of the _confidence_ in this estimation with the bootstrap. As described in the lecture slides, we will resample the test set with replacement, to measure n evaluations, and create a _distribution_ of f1-scores. \n",
    "\n",
    "Resampling with replacement is done with the [`resample()`](https://scikit-learn.org/stable/modules/generated/sklearn.utils.resample.html) function in sklearn:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.utils import resample\n",
    "\n",
    "seed = 1337\n",
    "pipe = make_pipeline(StandardScaler(), RandomForestClassifier(random_state=seed))\n",
    "pipe = pipe.fit(X_train_val, y_train_val)\n",
    "fprs = []\n",
    "for i in range(100):\n",
    "    X_test_sample, y_test_sample = resample(X_test, y_test, random_state=seed + i)\n",
    "    fpr = evaluate(pipe, X_test_sample, y_test_sample)\n",
    "    fprs.append(fpr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ℹ️ Note that we varied the seeds for each `resample()` call... or each sample would be exactly the same! \n",
    "\n",
    "We can now plot these metric bootstrap resampled distributions, which are an approximation of the sampling distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "sns.set()\n",
    "\n",
    "# split the fprs into f1s, precisions, and recalls\n",
    "f1s, precisions, recalls = list(zip(*fprs))\n",
    "\n",
    "fig = plt.figure(dpi=120)\n",
    "ax = fig.add_subplot()\n",
    "ax.hist(precisions, bins=50, range=[0,1], alpha=0.8, label='precision')\n",
    "ax.hist(recalls, bins=50, range=[0,1], alpha=0.8, label='recall')\n",
    "ax.hist(f1s, bins=50, range=[0,1], alpha=0.8, label='f1-score')\n",
    "ax.legend()\n",
    "ax.set_title('Bootstrap Distribution of F1-score, precision, and recall for Random Forest');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's a large variation in f1 scores! We can use [`np.percentile`](https://numpy.org/doc/stable/reference/generated/numpy.percentile.html) to approximate the confidence intervals:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def confidence_intervals(distribution, alpha=0.95):\n",
    "    top_p = ((1.0-alpha)/2.0) * 100\n",
    "    lower_limit = max(0.0, np.percentile(f1s, top_p))\n",
    "    bottom_p = (alpha+((1.0-alpha)/2.0)) * 100\n",
    "    upper_limit = min(1.0, np.percentile(f1s, bottom_p))\n",
    "    return lower_limit, upper_limit\n",
    "\n",
    "lower_limit, upper_limit = confidence_intervals(f1s, 0.95)\n",
    "\n",
    "print(f'F1-score estimate: {np.array(f1s).mean():.3f}, 95%CI: [{lower_limit:.3f}, {upper_limit:.3f}]')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our F1-score estimate isn't very precise. This is likely due to the fact that there are only 40 positive examples in our test set. 🙄 If we could acquire a larger labeled dataset, we would expect this confidence interval to be smaller. 💸\n",
    "\n",
    "🧠🧠 Would it be a good idea to use cross-validation to reduce this variance? why?\n",
    "\n",
    "There are other sources of error in the F1-score estimate. For example, the fact that bootstrap resampling the test set doesn't account for variations in the training set. Or that the model is trained on only 80% of the data, whereas we would ship a model retrained on 100% of the data to production. This creates a pessimistic bias in the metric measurement. More [advanced techniques](https://arxiv.org/pdf/1811.12808.pdf) exist to try solve these issues (e.g resampling training set too, and using bootstrap .632+ estimator). However, these biases and variances are tricky to eliminate, and an entire active statistical field of research is dedicated to these theoretical and empirical investigations. In the meantime, the methods outline in this lesson do a fair job of assessing model quality in practice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'Random Forest tasty wine detector metrics: f1-score: {np.array(f1s).mean():.3f}, precision: {np.array(precisions).mean():.3f}, recall: {np.array(recalls).mean():.3f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "🧠🧠 Based on these final results, could this tasty wine detector be useful in practice? For which usecases? Think of the type of errors it would make, based on the precision and recall values.\n",
    "\n",
    "💪💪 Estimate the confidence intervals of test f1-score for an SVM tasty wine detector.\n",
    "* Measure the f1-score, precision, and recall of an RBF kernel SVM with `C=100`\n",
    "* Use bootstrap resampling to estimate and plot their sampling distributions\n",
    "* store the f1-scores in a list called `f1s`\n",
    "* Then use percentiles to estimate their 95% confidence intervals, with the `confidence_intervals()` function\n",
    "* use a seed of `random_state=42` for all random operations. Increment this seed by one for the `random_state` of each successive bootstrap resample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# INSERT YOUR CODE HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "🧠🧠 Do you expect this SVM to perform differently to the random forest in a production setting? In what ways? Which one is most \"useful\"? How could we have incorporated this in our single-number metric?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Other Frameworks\n",
    "\n",
    "### 6.1 Metrics\n",
    "Keras supports all kinds of [metrics](https://www.tensorflow.org/api_docs/python/tf/keras/metrics), and they can be added to the model during the [compilation stage](https://keras.io/api/metrics/#usage-with-compile-amp-fit). \n",
    "\n",
    "Pytorch is a low-level library, and does not implement metrics. However, the sklearn metrics above can be used [directly](https://discuss.pytorch.org/t/calculating-precision-recall-and-f1-score-in-case-of-multi-label-classification/28265/3). Higher level apis like [ignite](https://pytorch.org/ignite/metrics.html) also implement various metrics.\n",
    "\n",
    "### 6.2 Dataset Splitting\n",
    "\n",
    "Keras and pytorch are deep learning libraries. Deep learning models are typically data-hungry, so there is less often a need for stratification or cross-validation. If needed however, the sklearn splitter classes (`KFold`, `StratifiedKFold`) can be used. \n",
    "\n",
    "Validation sets can directly be fed into keras models with the `validation_data` argument to the [`.fit()`](https://keras.io/api/models/model_training_apis/#fit-method) method.\n",
    "Splitters can be used directly in custom `Dataset` and `Dataloader` [classes](https://pytorch.org/tutorials/beginner/data_loading_tutorial.html) for pytorch.\n",
    "\n",
    "### 6.3 Confidence Intervals\n",
    "\n",
    "Likewise, we won't cover confidence intervals in keras or pytorch, since the exact same methods as section 5 can be used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. Summary\n",
    "\n",
    "Today, we learned about **advanced evaluation methods**. First, we defined the **imbalanced learning** regime, where some classes are rare in the dataset. We showed that accuracy is a bad metric for model quality when the classes are imbalanced, and how **balanced accuracy** can overcome these shortcomings. We then described a binary classification task mostly interested in finding the positive class. For these usecases, we showed that **precision**, **recall**, and **f1-score** metrics are good indicators of model quality. We alluded to different types of model error, i.e **false positives** and **false negatives**, and how single-number metrics must be designed with care to reflect the relative importance of these mistakes in the imbalanced learning regime. We saw how **stratification** can help maintain coherent class distributions across dataset splits, helping with the variance of metric estimates. **Cross-validation** was shown as a more efficient use of the dataset for training & validation. It rotates independent splits into **K-folds**, so that all the data is used for validation. Despite these techniques, metric estimates can still have high **variance**, and we mentioned several methods to estimate **confidence intervals**: binomial experiment MLE for accuracy, bootstrap resampling for held-out metrics, and Gaussian MLE for cross-validation. We highlighted that all these have some drawbacks, either in the expense of repeated computation, or in the accuracy of the estimations, and that model evaluation is an active research topic. Finally, we tried some of these techniques on a **tasty wine detection** task, where the rarity of the positive class and the dataset size presented challenges to assess model quality. \n",
    "\n",
    "\n",
    "# Resources\n",
    "\n",
    "## Core Resources\n",
    "\n",
    "- [**Slides**](https://docs.google.com/presentation/d/1HXO3G2Z9fAHOV7itzJ7jmUbggEDNZrg99rB7P6d_1OE/edit?usp=sharing)\n",
    "- [Intuition through visualization, ROC & AUC](https://www.pragmatic.ml/intuition-through-visualization-roc-auc/)  \n",
    "Visual explanation of ROC curves and AUC, metrics designed to measure the expertise of a model for a given problem.\n",
    "- [Machine learning yearning](https://www.deeplearning.ai/machine-learning-yearning/)  \n",
    "The reference for practical ML development\n",
    "- [sklearn on evaluation](https://scikit-learn.org/stable/modules/model_evaluation.html)  \n",
    "Official documentation on all sklearn evaluation methods\n",
    "- [raschka](https://arxiv.org/pdf/1811.12808.pdf)  \n",
    "Comprehensive and extremely detailed analysis of modern ML evaluation methods.\n",
    "\n",
    "\n",
    "### Additional Resources\n",
    "\n",
    "- [Google developer crash course - ROC AUC](https://developers.google.com/machine-learning/crash-course/classification/roc-and-auc)\n",
    "Short and sweet explanation of ROC curves and AUC \n",
    "- [Do machines actually beat doctors](https://lukeoakdenrayner.wordpress.com/2017/12/06/do-machines-actually-beat-doctors-roc-curves-and-performance-metrics/)  \n",
    "Practical demonstration of the value of ROC curves and AUC in the medical domain\n",
    "- [Confidence intervals for ML](https://machinelearningmastery.com/confidence-intervals-for-machine-learning/)  \n",
    "Solid explanation and coding examples of confidence estimation for ML models\n",
    "- [Estimation statistics should replace significance testing](https://www.nature.com/articles/nmeth.3729)  \n",
    "Classic paper about confidence intervals vs p-values and paired tests\n",
    "- [Wikipedia - binomial propotion confidence interval](https://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval)  \n",
    "Parametric method for confidence interval estimation of binomial experiments, like accuracy\n",
    "- [No unbiased estimator of the variance of k-fold crossvalidation](http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.9.3582)  \n",
    "Classic paper about the biases in cross-validation standard deviation as confidence estimation\n",
    "- [Why question machine learning evaluation methods?](https://www.aaai.org/Papers/Workshops/2006/WS-06-06/WS06-06-003.pdf)  \n",
    "Overview of limitations of current ML evaluation methods\n",
    "- [imbalanced-learn](https://github.com/scikit-learn-contrib/imbalanced-learn)  \n",
    "sklearn ecosystem framework for imbalanced learning methods\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}