09_classification_metrics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Topics\n",
    "1. Review of model evaluation\n",
    "2. Model evaluation procedures\n",
    "3. Model evaluation metrics\n",
    "4. Classification accuracy\n",
    "5. Confusion matrix\n",
    "6. Metrics computed from a confusion matrix\n",
    "7. Adjusting the classification threshold\n",
    "8. Receiver Operating Characteristic (ROC) Curves\n",
    "9. Area Under the Curve (AUC)\n",
    "10. Confusion Matrix Resources\n",
    "11. ROC and AUC Resources\n",
    "12. Other Resources"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_This tutorial is derived from Data School's Machine Learning with scikit-learn tutorial. I added my own notes so anyone, including myself, can refer to this tutorial without watching the videos._"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Review of model evaluation\n",
    "\n",
    "- Need a way to choose between models: different model types, tuning parameters, and features\n",
    "- Use a **model evaluation procedure** to estimate how well a model will generalize to out-of-sample data\n",
    "- Requires a **model evaluation metric** to quantify the model performance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Model evaluation procedures\n",
    "\n",
    "1. **Training and testing on the same data**\n",
    "    - Rewards overly complex models that \"overfit\" the training data and won't necessarily generalize\n",
    "    \n",
    "2. **Train/test split**\n",
    "    - Split the dataset into two pieces, so that the model can be trained and tested on different data\n",
    "    - Better estimate of out-of-sample performance, but still a \"high variance\" estimate\n",
    "    - Useful due to its speed, simplicity, and flexibility\n",
    "3. **K-fold cross-validation**\n",
    "    - Systematically create \"K\" train/test splits and average the results together\n",
    "    - Even better estimate of out-of-sample performance\n",
    "    - Runs \"K\" times slower than train/test split"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Model evaluation metrics\n",
    "\n",
    "- **Regression problems:** Mean Absolute Error, Mean Squared Error, Root Mean Squared Error\n",
    "- **Classification problems:** Classification accuracy\n",
    "    - There are many more metrics, and we will discuss them today"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Classification accuracy\n",
    "\n",
    "[Pima Indian Diabetes dataset](https://archive.ics.uci.edu/ml/datasets/Pima+Indians+Diabetes) from the UCI Machine Learning Repository"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# read the data into a Pandas DataFrame\n",
    "import pandas as pd\n",
    "\n",
    "url = 'https://archive.ics.uci.edu/ml/machine-learning-databases/pima-indians-diabetes/pima-indians-diabetes.data'\n",
    "\n",
    "col_names = ['pregnant', 'glucose', 'bp', 'skin', 'insulin', 'bmi', 'pedigree', 'age','label']\n",
    "\n",
    "pima = pd.read_csv(url, header=None, names=col_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>pregnant</th>\n",
       "      <th>glucose</th>\n",
       "      <th>bp</th>\n",
       "      <th>skin</th>\n",
       "      <th>insulin</th>\n",
       "      <th>bmi</th>\n",
       "      <th>pedigree</th>\n",
       "      <th>age</th>\n",
       "      <th>label</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>6</td>\n",
       "      <td>148</td>\n",
       "      <td>72</td>\n",
       "      <td>35</td>\n",
       "      <td>0</td>\n",
       "      <td>33.6</td>\n",
       "      <td>0.627</td>\n",
       "      <td>50</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>85</td>\n",
       "      <td>66</td>\n",
       "      <td>29</td>\n",
       "      <td>0</td>\n",
       "      <td>26.6</td>\n",
       "      <td>0.351</td>\n",
       "      <td>31</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>8</td>\n",
       "      <td>183</td>\n",
       "      <td>64</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>23.3</td>\n",
       "      <td>0.672</td>\n",
       "      <td>32</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>89</td>\n",
       "      <td>66</td>\n",
       "      <td>23</td>\n",
       "      <td>94</td>\n",
       "      <td>28.1</td>\n",
       "      <td>0.167</td>\n",
       "      <td>21</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0</td>\n",
       "      <td>137</td>\n",
       "      <td>40</td>\n",
       "      <td>35</td>\n",
       "      <td>168</td>\n",
       "      <td>43.1</td>\n",
       "      <td>2.288</td>\n",
       "      <td>33</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   pregnant  glucose  bp  skin  insulin   bmi  pedigree  age  label\n",
       "0         6      148  72    35        0  33.6     0.627   50      1\n",
       "1         1       85  66    29        0  26.6     0.351   31      0\n",
       "2         8      183  64     0        0  23.3     0.672   32      1\n",
       "3         1       89  66    23       94  28.1     0.167   21      0\n",
       "4         0      137  40    35      168  43.1     2.288   33      1"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# print the first 5 rows of data from the dataframe\n",
    "pima.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- label\n",
    "    - 1: diabetes\n",
    "    - 0: no diabetes\n",
    "- pregnant\n",
    "    - number of times pregnant"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question:** Can we predict the diabetes status of a patient given their health measurements?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# define X and y\n",
    "feature_cols = ['pregnant', 'insulin', 'bmi', 'age']\n",
    "\n",
    "# X is a matrix, hence we use [] to access the features we want in feature_cols\n",
    "X = pima[feature_cols]\n",
    "\n",
    "# y is a vector, hence we use dot to access 'label'\n",
    "y = pima.label"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# split X and y into training and testing sets\n",
    "from sklearn.cross_validation import train_test_split\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,\n",
       "          intercept_scaling=1, max_iter=100, multi_class='ovr', n_jobs=1,\n",
       "          penalty='l2', random_state=None, solver='liblinear', tol=0.0001,\n",
       "          verbose=0, warm_start=False)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# train a logistic regression model on the training set\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "\n",
    "# instantiate model\n",
    "logreg = LogisticRegression()\n",
    "\n",
    "# fit model\n",
    "logreg.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# make class predictions for the testing set\n",
    "y_pred_class = logreg.predict(X_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Classification accuracy:** percentage of correct predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.692708333333\n"
     ]
    }
   ],
   "source": [
    "# calculate accuracy\n",
    "from sklearn import metrics\n",
    "print(metrics.accuracy_score(y_test, y_pred_class))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Classification accuracy is 69%"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Null accuracy:** accuracy that could be achieved by always predicting the most frequent class\n",
    "- We must always compare with this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    130\n",
       "1     62\n",
       "Name: label, dtype: int64"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# examine the class distribution of the testing set (using a Pandas Series method)\n",
    "y_test.value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.3229166666666667"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# calculate the percentage of ones\n",
    "# because y_test only contains ones and zeros, we can simply calculate the mean = percentage of ones\n",
    "y_test.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "32% of the "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6770833333333333"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# calculate the percentage of zeros\n",
    "1 - y_test.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6770833333333333"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# calculate null accuracy in a single line of code\n",
    "# only for binary classification problems coded as 0/1\n",
    "max(y_test.mean(), 1 - y_test.mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This means that a dumb model that always predicts 0 would be right 68% of the time\n",
    "- This shows how classification accuracy is not that good as it's close to a dumb model\n",
    "- It's a good way to know the minimum we should achieve with our models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    0.677083\n",
       "Name: label, dtype: float64"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# calculate null accuracy (for multi-class classification problems)\n",
    "y_test.value_counts().head(1) / len(y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comparing the **true** and **predicted** response values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True: [1 0 0 1 0 0 1 1 0 0 1 1 0 0 0 0 1 0 0 0 1 1 0 0 0]\n",
      "False: [0 0 0 0 0 0 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0]\n"
     ]
    }
   ],
   "source": [
    "# print the first 25 true and predicted responses\n",
    "print('True:', y_test.values[0:25])\n",
    "print('False:', y_pred_class[0:25])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Conclusion:**\n",
    "\n",
    "- Classification accuracy is the **easiest classification metric to understand**\n",
    "- But, it does not tell you the **underlying distribution** of response values\n",
    "    - We examine by calculating the null accuracy\n",
    "- And, it does not tell you what **\"types\" of errors** your classifier is making"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Confusion matrix\n",
    "\n",
    "Table that describes the performance of a classification model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[118  12]\n",
      " [ 47  15]]\n"
     ]
    }
   ],
   "source": [
    "# IMPORTANT: first argument is true values, second argument is predicted values\n",
    "# this produces a 2x2 numpy array (matrix)\n",
    "print(metrics.confusion_matrix(y_test, y_pred_class))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Small confusion matrix](images/09_confusion_matrix_1.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Every observation in the testing set is represented in **exactly one box**\n",
    "- It's a 2x2 matrix because there are **2 response classes**\n",
    "- The format shown here is **not** universal\n",
    "    - Take attention to the format when interpreting a confusion matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Basic terminology**\n",
    "\n",
    "- **True Positives (TP):** we *correctly* predicted that they *do* have diabetes\n",
    "    - 15\n",
    "- **True Negatives (TN):** we *correctly* predicted that they *don't* have diabetes\n",
    "    - 118\n",
    "- **False Positives (FP):** we *incorrectly* predicted that they *do* have diabetes (a \"Type I error\")\n",
    "    - 12\n",
    "    - Falsely predict positive\n",
    "    - Type I error\n",
    "- **False Negatives (FN):** we *incorrectly* predicted that they *don't* have diabetes (a \"Type II error\")\n",
    "    - 47\n",
    "    - Falsely predict negative\n",
    "    - Type II error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 0: negative class\n",
    "- 1: positive class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True [1 0 0 1 0 0 1 1 0 0 1 1 0 0 0 0 1 0 0 0 1 1 0 0 0]\n",
      "Pred [0 0 0 0 0 0 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0]\n"
     ]
    }
   ],
   "source": [
    "# print the first 25 true and predicted responses\n",
    "print('True', y_test.values[0:25])\n",
    "print('Pred', y_pred_class[0:25])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[118  12]\n",
      " [ 47  15]]\n"
     ]
    }
   ],
   "source": [
    "# save confusion matrix and slice into four pieces\n",
    "confusion = metrics.confusion_matrix(y_test, y_pred_class)\n",
    "print(confusion)\n",
    "#[row, column]\n",
    "TP = confusion[1, 1]\n",
    "TN = confusion[0, 0]\n",
    "FP = confusion[0, 1]\n",
    "FN = confusion[1, 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Large confusion matrix](images/09_confusion_matrix_2.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6. Metrics computed from a confusion matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Classification Accuracy:** Overall, how often is the classifier correct?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.692708333333\n",
      "0.692708333333\n"
     ]
    }
   ],
   "source": [
    "# use float to perform true division, not integer division\n",
    "print((TP + TN) / float(TP + TN + FP + FN))\n",
    "print(metrics.accuracy_score(y_test, y_pred_class))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Classification Error:** Overall, how often is the classifier incorrect?\n",
    "\n",
    "- Also known as \"Misclassification Rate\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.307291666667\n",
      "0.307291666667\n"
     ]
    }
   ],
   "source": [
    "classification_error = (FP + FN) / float(TP + TN + FP + FN)\n",
    "\n",
    "print(classification_error)\n",
    "print(1 - metrics.accuracy_score(y_test, y_pred_class))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Sensitivity:** When the actual value is positive, how often is the prediction correct?\n",
    "- Something we want to maximize\n",
    "- How \"sensitive\" is the classifier to detecting positive instances?\n",
    "- Also known as \"True Positive Rate\" or \"Recall\"\n",
    "- TP / all positive\n",
    "    - all positive = TP + FN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.241935483871\n",
      "0.241935483871\n"
     ]
    }
   ],
   "source": [
    "sensitivity = TP / float(FN + TP)\n",
    "\n",
    "print(sensitivity)\n",
    "print(metrics.recall_score(y_test, y_pred_class))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Specificity:** When the actual value is negative, how often is the prediction correct?\n",
    "- Something we want to maximize\n",
    "- How \"specific\" (or \"selective\") is the classifier in predicting positive instances?\n",
    "- TN / all negative \n",
    "    - all negative = TN + FP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.907692307692\n"
     ]
    }
   ],
   "source": [
    "specificity = TN / (TN + FP)\n",
    "\n",
    "print(specificity)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our classifier\n",
    "- Highly specific\n",
    "- Not sensitive"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**False Positive Rate:** When the actual value is negative, how often is the prediction incorrect?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0923076923077\n",
      "0.0923076923077\n"
     ]
    }
   ],
   "source": [
    "false_positive_rate = FP / float(TN + FP)\n",
    "\n",
    "print(false_positive_rate)\n",
    "print(1 - specificity)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Precision:** When a positive value is predicted, how often is the prediction correct?\n",
    "\n",
    "- How \"precise\" is the classifier when predicting positive instances?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.555555555556\n",
      "0.555555555556\n"
     ]
    }
   ],
   "source": [
    "precision = TP / float(TP + FP)\n",
    "\n",
    "print(precision)\n",
    "print(metrics.precision_score(y_test, y_pred_class))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Many other metrics can be computed: F1 score, Matthews correlation coefficient, etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Conclusion:**\n",
    "\n",
    "- Confusion matrix gives you a **more complete picture** of how your classifier is performing\n",
    "- Also allows you to compute various **classification metrics**, and these metrics can guide your model selection\n",
    "\n",
    "**Which metrics should you focus on?**\n",
    "\n",
    "- Choice of metric depends on your **business objective**\n",
    "    - Identify if FP or FN is more important to reduce\n",
    "    - Choose metric with relevant variable (FP or FN in the equation)\n",
    "- **Spam filter** (positive class is \"spam\"): \n",
    "    - Optimize for **precision or specificity** \n",
    "        - precision\n",
    "            - false positive as variable\n",
    "        - specificity\n",
    "            - false positive as variable\n",
    "    - Because false negatives (spam goes to the inbox) are more acceptable than false positives (non-spam is caught by the spam filter)\n",
    "- **Fraudulent transaction detector** (positive class is \"fraud\"): \n",
    "    - Optimize for **sensitivity**\n",
    "        - FN as a variable\n",
    "    - Because false positives (normal transactions that are flagged as possible fraud) are more acceptable than false negatives (fraudulent transactions that are not detected)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  7. Adjusting the classification threshold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 0, 0, 0, 0, 0, 0, 1, 0, 1])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# print the first 10 predicted responses\n",
    "# 1D array (vector) of binary values (0, 1)\n",
    "logreg.predict(X_test)[0:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.63247571,  0.36752429],\n",
       "       [ 0.71643656,  0.28356344],\n",
       "       [ 0.71104114,  0.28895886],\n",
       "       [ 0.5858938 ,  0.4141062 ],\n",
       "       [ 0.84103973,  0.15896027],\n",
       "       [ 0.82934844,  0.17065156],\n",
       "       [ 0.50110974,  0.49889026],\n",
       "       [ 0.48658459,  0.51341541],\n",
       "       [ 0.72321388,  0.27678612],\n",
       "       [ 0.32810562,  0.67189438]])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# print the first 10 predicted probabilities of class membership\n",
    "logreg.predict_proba(X_test)[0:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Row: observation\n",
    "    - Each row, numbers sum to 1\n",
    "- Column: class\n",
    "    - 2 reponse classes there 2 columns\n",
    "        - column 0: predicted probability that each observation is a member of class 0\n",
    "        - column 1: predicted probability that each observation is a member of class 1 \n",
    "- Importance of predicted probabilities\n",
    "    - We can rank observations by probability of diabetes\n",
    "        - Prioritize contacting those with a higher probability\n",
    "- predict_proba process\n",
    "    1. Predicts the probabilities\n",
    "    2. Choose the class with the highest probability \n",
    "- There is a 0.5 classification threshold\n",
    "    - Class 1 is predicted if probability > 0.5\n",
    "    - Class 1 is predicted if probability < 0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.36752429,  0.28356344,  0.28895886,  0.4141062 ,  0.15896027,\n",
       "        0.17065156,  0.49889026,  0.51341541,  0.27678612,  0.67189438])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# print the first 10 predicted probabilities for class 1\n",
    "logreg.predict_proba(X_test)[0:10, 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# store the predicted probabilities for class 1\n",
    "y_pred_prob = logreg.predict_proba(X_test)[:, 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# allow plots to appear in the notebook\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# adjust the font size \n",
    "plt.rcParams['font.size'] = 12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x11c2b1128>"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEgCAYAAABFO1+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcHVWd9/HPNwmGNRDWyJJEEBTiQERkBAI0oIIoAjIC\nboAwjgsPqIg6DyCJbG6MOo+i+Iggi4yggKjogAItqAMDLjhsIhDWhCWQlRCC5Dd/nHND5aa7695O\n963q7u/79bqv7ltVt86vTtWt361zalFEYGZm1pdRVQdgZmb152RhZmalnCzMzKyUk4WZmZVysjAz\ns1JOFmZmVsrJoiKSZko6qeo4hhJJoyWdL2mOpJck7VF1TACSJklaJmnXnt5XEM90SX+rouxc/lGS\nXhyA+Vwg6bqSaaZLuq+3siXtmbeVTduZj63MyWIA9bVx553HewuDdgK+1uJ8d8ufnzgQcQ5hhwCH\nA28HXgn8vtpwVlC8YOkRYAJwaysfHKT1W+UFVNHB8r8CvKmPsn8HvDIiZkGfdd08H2sypuoARqqI\neKaNyUUHv/ySVouIVf5lOAi2AR6PiJZ2wu0YgGVW459IV7o+1eZna391bB23i4hYDCzuY/zfWXFd\n9FjXZfMxH1lUprkZStKBkv4o6TlJcyXdImkHSZOAm/JkD+VfRTcUPneipAckvSDpfkkfbypnfUk/\nkrRI0ixJp+YjoF8VprlR0nmSTpM0C3g4D39PjmOepKcl/VzS1oXPNZpb3iPpP3Ps90jaQ9Kmkq7J\n5d4laVoLddLrski6ETgN2DKX+WAv82jE9D5Jv5a0OM/zsB6meW8jxjxvJL1a0o/zOnhW0rWSXtdU\nxqGS/ibpeUm/BbbvJYZdC8M2yvX+RP7cPbnJpGz9Hi7pT/kzMyX9m6Q1C+PHSvp2XkfPSPoWMLaF\nul4m6fi8rIskPSbp+B6mOU7SDyTNAy7Kw1+T621hfv1U0lY9lLGPpDtz7LdI2qEwbj1JF0t6OK+j\neyWd0Eusn8jxPSfpcknjC+P6bHJTaoZalrfHXuta0ozm+Uh6i6Tf5vgeU2oCXb8wfru83c8tbOfv\n66Pah7aI8GuAXsAFwHW9jFsGvLfwfiZwUv5/E+AF4FPAJOA1pOaWKaRfQgcALwE7AhsD6+XPHQs8\nBxwDbAX8C/A88MFCOT8F7gX2ALYFzgfmFuMEbgTmA98CXgtMycOPJDX5TAZ2AH4C3AeMyeMn5eX6\nW47x1cCVwCzgOuDAPOxHpAQ0uo+663NZgPVITQUPABsBG/Qyn0ZMj+U63Bo4Hfg7sEPTNI8A78nv\nJ+W6nQ18E9guf/bfgacb5QGvz/M6I48/CHgwr59dC/Mvvl8duAe4Hdgrj98LeHfJ+j0KeAZ4b/7M\nNODPwIWF5f0a8ATwDtKR11fyuryvZFtdBswBPpbX0XHAi8ABTdM8nad5VV4vq+d1+Stgaq6PG5q2\niyPz8tyeY34d8DPgUWBsYZv/DGm7mpSXcQFwZNP3aT5pu9uOtA3fB1xRmGZ6cVlz2UsL7/fMsWxK\n+nHcW103z2dv0vb4MWBL4A3A9UB3YZo7gEtI39fJwL7A/lXvhwZt/1Z1AMPplTfuF4GFPbxeovdk\nMTWPn9jLfHfraTxpZ/eFpmFfBe7P/2+dv/BdhfFj8ueak8W9LSzf+nl+u+T3jZ3ucYVpdsrDPlEY\n1li+7fqYd5/Lkt+v8IXuZT6NmGY0Df8deSdbmOakpmmmA79vGibgfuD4/P4S4OamaY5l5WSxrPD+\nGFITxyvbXL8zgX9pGrZ7nve6wJqkhHp00zS3tVBPy4DvNw37AfCbpmn+f9M0xwCLgPGFYRvn5Xt/\nft9IFsXtbj3S9+CDfcT0deDapu/TAmDtwrC35Li27GmboI9kUVLXzfO5ETiraZqJuezt8/t5wBFl\n35vh8nIz1MC7hdQssUPTS3185i+kX+J3SboyNw9s3lchktYBNgdubhr1G2CypNVJRxJBoaM1Uhvu\n7T3M8g89lDE1x/OgpAWkX5RB2hk2x9/wRP77P03DRNqp9HdZ2nVL0/vfkY7Uim5rev9GYKdC88pC\n0s5qEunXN6Q6be5Y/y19r98dgbsjYnarwUvaMJf71aZ4fknaYb2a9Ev/FcB/9RBPK/pTR9uRlmVu\nY0BEPAX8tYfP3lKYZh7p6GoKgJJ/zU1sT+dl+wgrb1t3R8SiphgbcQymNwKfaKr7u0jbf6Mp9mzg\ne0rNuNMlvX6QY6qUO7gH3vMRMbN5oNT7viQilgFvk7QT8GbSWT9flPRPEfGLAYiplc7T54pvJK0B\nXEvagR8FPJlH3U3aQRUVOz2jj2F1+3HyXNP7UcCvSUcKzStsfkciWjEWgOOB7h7GP0Zq/ugrSQ2E\n5joaKCcCnwU+QWpaWwicAOw/SOW1axTwJeDiHsY9ARARZ0i6BNiP1Gx1kqQvRcSpnQuzc+r25R3R\nIuL2iPhiROxJ+lX9wTxqaf47ujDtQtIOo/lagy5gZkQsIe3YAXZpjJQ0mtT+WmZbYEPg5Ii4KSL+\nCmzAIOycWlyWdjWfBrkrL9dHb24n/fJ9PCIebHo1zl67O8+raBp9J+Q/ANup93P9e1q/T5Ha+F/b\nQywPRsRSUv/N0h7i2a1kORua62g3yuvoLtKyFDt6NyElrv9pmvZNhWnWI21Td+VBuwP/GREXRsQd\nEfEgqc+l2baS1m6KMVqIszcr1XUvbif13fVU98vPmoqIhyLi3Ig4FDgV+Gg/46o9J4sakLSLpFMk\n7SxpC0n7kJqyGl+sh0lND/vns2rG5eFfAI6T9M9KZ/F8GPgwcCZARNwP/Bw4R+kMpW2B7wDjKD/a\neJjU6X68pC1zTF/PcQyGPpelH45ROktra0mnkXZc/1bymW+SdiI/lTRN6aymaZLOkNTY8X0N2CUP\n21rSwaRfxH35D1J9/jSfITRZ0t6SDs3je1u/J5Pq/yRJUyRtI+kgSefC8tM9zwXOkHRAHv8l0o67\nFe+QdGyu7+NIHe5nl3zmUlLH+GWSXi/pDcAPSYnt8qZpvyxpd0n/QDqTakGuC0jNVl2SunI9ng7s\n3EuZF+Xl34O0jq7OyaVVxR84vdV1s1OBA5XOPtshfwf2UzprcKyktSR9U9JeeX2+nnSEcVcv8xv6\nqu40GU4v+j4bqrmD+0Fe7uDeDriGdBbR86SOzS+Szy7J05xI+kK+CNxQGP4p0i/MF0gdscc1lTue\n9CVeRDp8/nx+f3Vhmhto6sjMw99F+lIvJv063p30y+yIPH6Fs37ysM3ysD0KwzbJw/Yuqb+yZWmn\ng/t9pE7KxXmehzVNs0LchXFbkJoeniysi4uASYVpDiWdAfY8qb/gAPo4GyoP2xj4Pumc/8WkX8ZH\nFMb3tn7fSWqnX0TqUP0jcEph/OrAt0lnuM0lJY8zW6inZaQmrqtITU2PAx/va5stDN+a9CNkQX5d\nTe5wzuOPzNvJm/NyNupph8I040hJZh7pjKtv5G3zwebvEykZz8p1cDkrdq631cHdW133tG2RjmKu\nIzVBNvosvkr6kT2WdELAA3l9PkFKhJtVvR8arJdypXSEpG7gH/NKEvBYRGybx+1D+tWwBalD9oMR\n8UjHghshJI0inUp7dUR8uup4BprSufQzgWkRUacrvGtF0jLS2UuXVh2LDQ2dboYK4GMRMS4i1ikk\nig2AK0iH3euTfsVe1uHYhqXcDHBIPoyeSvq1Non0K9fMrCVVnA3VUwfpu4A7I+JKSFdTAnMkbRMR\nvrnXqhkNnEI6zfJF4E7S+e/Dt211CNw6owZcR9aWTjdD3UhqnxepLfyUiPiNpK8Dq0XEsYVp/wJM\nj4irOhagmZn1qNNHFp8hdXgtJd1m4ae5aWRtVr7x2gJgnc6GZ2ZmPelosoiI4tWgF0lq3G56Eens\niKJ1SWcgrECSD5/NzPohIvp9nVRdrrO4i3T/IAAkrUVqY++xXb3qU8jq8po+fXrlMdTl5bpwXbgu\n+n6tqo4lC0nrSnprvqBltNKtfHcn3evmKmCKpIMljSWd8/zncOe2mVktdLIZajXSbZ1fQ7pI5l7g\nwIh4AEDSIcA5pLt63kq6vbSZmdVAx5JFRMyh98v5iYgbSPeOsRZ1dXVVHUJtuC5e5rp4meti4HT0\n1NmBICmGWsxmZlWTRAyDDm4zM6sxJwszMyvlZGFmZqWcLMzMrJSThZmZlXKyMDOzUk4WZmZWysnC\nzMxKOVmYmVkpJwszMyvlZGFmZqWcLIahCRMmI6nS14QJk6uuBjMbQL6R4DAkCai6jjQgD1wxs4Hh\nGwmamdmgc7IwM7NSThZmZlbKycLMzEo5WZiZWSknCzMzK+VkYWZmpZwszMyslJOFmZmVcrIwM7NS\nThZmZlbKycLMzEqNqTqA4WjChMk8+eTDVYdhZjZgfNfZQVD9XV+rLj/FUPf1ZDaS+K6zZmY26Jws\nzMyslJOFmZmVcrIwM7NSThZmZlbKycLMzEpVkiwkbS3peUkXFYbtI+keSYskXS9pYhWxmZnZyqo6\nsvgm8N+NN5I2BK4ATgbWB/4AXFZNaGZm1qzjyULS4cBc4PrC4IOBOyPiyohYCswAdpC0TafjMzOz\nlXU0WUgaB3weOIF0mXHDFOCOxpuIWAzcn4ebmVnFOn1kcRrw3YiY1TR8bWB+07AFwDodicrMzPrU\nsRsJSpoKvBmY2sPoRcC4pmHrAgt7mteMGTOW/9/V1UVXV9eAxGhmNlx0d3fT3d09YPPr2I0EJX0c\nOIOUAEQ6mhgF3AOcCxwVEdPytGsBTwNTI+K+pvn4RoLlEVRcfoqh7uvJbCRZ1RsJdjJZrM6KRw+f\nBiYBHyEljb8BRwO/AE4HpkXErj3Mx8miPIKKy08x1H09mY0kq5osOtYMFRFLgCWN95IWAUsi4tn8\n/hDgHOAS4Fbg8E7FZmZmffPzLAaBjyxSDHVfT2YjiZ9nYWZmg87JwszMSjlZmJlZKScLMzMr5WRh\nZmalnCzMzKyUk4WZmZVysjAzs1JOFmZmVsrJwszMSjlZmJlZKScLMzMr5WRhZmalnCxsWJowYTKS\nKn1NmDC56mowGzC+Rfkg8C3KUwxVrqfq1wFUXQdmRb5FuZmZDTonCzMzK+VkYWZmpZwszMyslJOF\nmZmVcrIwM7NSThZmZlbKycLMzEo5WZiZWSknCzMzK+VkYWZmpZwszMyslJOFmZmVcrIwM7NSThZm\nZlbKycLMzEqNqToAG67G5gcQmdlw4GRhg+QFqn9aoJkNFDdDmZlZqZaThaSPS9pwMIMxM7N6aufI\nYm/gIUk/l3SYpLHtFibpYkmzJc2TdK+kYwrj9pF0j6RFkq6XNLHd+ZuZ2eBQROvtypI2AA4H3g+8\nFrgCuCgibmrx89sBD0bEEknbAL8B9gceAR4AjgZ+DpwB7B4Ru/Qwj2gn5iqkjt2q2+urrqOqY6i6\n/BRD3bdVGzkkERH97sxrK1k0Fbw9cDHwOuBR4LvAv0fEohY//xrgRuB4YDxwZERMy+PWBOYAUyPi\nvqbPOVmUR1Bx+XWIoeryUwx131Zt5FjVZNF2B3duLroA6AaeBI4APgC8HvhlC58/R9JzwD3ALOAX\nwBTgjsY0EbEYuD8PNzOzirV86qyks0lNUPOBi4BTIuLxwvhbgLll84mIYyX9H2AXoAtYCqwNPNU0\n6QJgnVbjMzOzwdPOdRarAwdHxG09jYyIFyXt1MqMcjvS7yV9APgosAgY1zTZusDCnj4/Y8aM5f93\ndXXR1dXVSrFmZiNGd3c33d3dAza/lvssJG0GLI6IuYVh44E1ImJWvwqXvktKFHcBRxX6LNYCnsZ9\nFv2NoOLy6xBD1eWnGOq+rdrI0ck+i58AmzcN2xy4qpUPS9oon3K7lqRRkvYlNWv9Os97iqSD8ym5\n04E/NycKMzOrRjtHFgsiormpCEnzI2LdFj6/IfBjYHtSknqYdPbU+Xn83sA5wETgVtKRxiM9zMdH\nFuURVFx+HWKouvwUQ923VRs5OnbqrKT7gf0i4v7CsFcD10XElv0NoF1OFi1FUHH5dYih6vJTDHXf\nVm3k6GQz1PnAFZLeIWk7SQeQjhTO62/hZmY2NLRzNtQXgReBs4EtSBfinQd8dRDiMjOzGun3FdxV\ncTNUSxFUXH4dYqi6/BRD3bdVGzlWtRmqredZ5Ft07EC6iG65Rie1mZkNT+1cwX0ScCrpthyLC6OC\n1J9hZmbDVDtnQz0FvDki/jK4IZXG4Wao8ggqLr8OMVRdfoqh7tuqjRydPBvqeeDe/hZkZmZDVzvJ\n4nPANyS9Ml+Bvfw1WMGZmVk9tNMMtSz/W/yASPcFHD3QgfURh5uhyiOouPw6xFB1+SmGum+rNnJ0\n8myoV/W3EDMzG9ravs4iNzttEhGzByek0vJ9ZFEeQcXl1yGGqstPMdR9W7WRo2Md3JLWk3QpsIT0\nFDskvVPSGf0t3MzMhoZ2OqfPJT0lbxLp6XYA/wUcNtBBmZlZvbTTZ7EPsGl+Il4ARMTTkjYenNDM\nzKwu2jmymA9sWBwgaSJQSd+FmZl1TjvJ4jzSLcr3AkZJ2gW4kNQ8ZWZmw1g711kIOB74MKnf4hHg\nO6Sn3XXslA+fDdVSBBWXX4cYqi4/xVD3bdVGjo49Ka8unCxaiqDi8usQQ9Xlpxjqvq3ayNGxi/Ly\nM7J7FBE39DcAMzOrv3aaoWY2DdoIeAXwmJ/BvSIfWdQhhqrLTzHUfVu1kaNjRxYRscLtPiSNBk4B\nFva3cDMzGxpWqc9C0hjSkcWEgQuptEwfWZRHUHH5dYih6vJTDHXfVm3k6OTzLHryFmBZ6VRmZjak\ntdPB/Sgr/lRbE1gd+NhAB2VmZvXSzu0+3t/0/jngvohYMIDxmJlZDfk6i0HgPos6xFB1+SmGum+r\nNnJ08jqLi2nh2xcRR/Q3GDMzq6d2OrjnAQcBo4HH8mcPzMMfKLzMzGyYaafPYhvg7RFxc2OApGnA\n5yJi3wGPzMzMaqOdK7jnAxtGxIuFYasBz0TEuEGKr6c43GdRHkHF5dchhqrLTzHUfVu1kaOT11n8\nCThL0hq54DWAM4E/97dwMzMbGtpJFkcBuwHzJT1JehjSNODIQYjLzMxqpO1TZyVtAWwKzI6IRwYl\nqr7LdzNUeQQVl1+HGKouP8VQ923VRo6O3u5D0gZAF7BnRDwiaVNJm/e3cDMzGxpaThaS9gT+CrwP\n+FwevDXw7UGIy8zMaqSdI4uvA4dFxH7A3/OwW4GdW/mwpFdIOk/SQ5LmS/qjpP0K4/eRdI+kRZKu\nlzSxjdjMzGwQtZMsJkfE9fn/RkPsUlq/VmMM6bndu0fEuqSjk8slTczNW1cAJwPrA38ALmsjNjMz\nG0TtXJR3t6R9I+LawrA3A//TyocjYjFwWuH9Nfnpe28ANgTujIgrASTNAOZI2iYi7msjRjMzGwTt\nHFl8CviBpAuBNSR9B/g+8On+FCxpE1Kfx13AFOCOxricWO7Pw83MrGItJ4uIuAXYnrRzPx+YCewc\nEbe1W2h+wt4lwPfzkcPapOs2ihYA67Q7bzMzG3gtNUPl521fD+wbEV9elQKVLkK4BHgBOC4PXgQ0\n3zJkXXp5vveMGTOW/9/V1UVXV9eqhGRmNux0d3fT3d09YPNr595QDwOvjYjnV6lA6XxgIrB/RCzN\nwz4EHBkR0/L7tYCnganNfRa+KK+lCCouvw4xVF1+iqHu26qNHJ28KO/zwLclTZI0WtKoxqvVGUg6\nF3gt8M5GosiuAqZIOljSWGA68Gd3bpuZ1UM7RxbLCm8bHxIQETG6hc9PBB4ClgAvFebz4Yj4D0l7\nA+eQjjpuBY7q6XYiPrJoKYKKy69DDFWXn2Ko+7ZqI8eqHlmUJgtJEyLiCUmTepsmIh7ubwDtcrJo\nKYKKy69DDFWXn2Ko+7ZqI0cnksWC4vMqJF0ZEe/qb4GrysmipQgqLr8OMVRdfoqh7tuqjRyd6LNo\nnnlXfwszM7OhqZVk4Z9GZmYjXCvXWYyRtBcvH2E0vycibhiM4MzMrB5a6bN4iL6PLiIithzIoPri\nPouWIqi4/DrEUHX5KYa6b6s2cgx6B3fdOFm0FEHF5dchhqrLTzHUfVu1kaOjT8ozM7ORycnCzMxK\nOVmYmVkpJwszMyvlZGFmZqWcLMzMrJSThZmZlXKyMDOzUk4WZmZWysnCzMxKOVmYmVkpJwszMyvl\nZGFmZqWcLMzMrJSThZmZlXKyMDOzUk4WZmZWysnCzMxKOVmYmVkpJwszMyvlZGFmZqWcLMzMrJST\nhZmZlXKyMDOzUk4WZsPUhAmTkVTpa8KEyVVXgw0QRUTVMbRFUtQ9ZklAlTFWXX4dYqi6/BRDldtq\n9dshVF0H9jJJRIT6+3kfWZiZWSknCzMzK9XRZCHpWEm3SVoi6fymcftIukfSIknXS5rYydjMzKx3\nnT6yeBw4HfhecaCkDYArgJOB9YE/AJd1ODYzM+vFmE4WFhE/AZD0RmCzwqh3AXdGxJV5/AxgjqRt\nIuK+TsZoZmYrq0ufxRTgjsabiFgM3J+Hm5lZxeqSLNYG5jcNWwCsU0EsZmbWpKPNUH1YBIxrGrYu\nsLCniWfMmLH8/66uLrq6ugYrLjOzIam7u5vu7u4Bm18lF+VJOh3YLCKOzu8/BBwZEdPy+7WAp4Gp\nzX0WviivpQgqLr8OMVRdforBF+X5ory6GFIX5UkaLWl1YDQwRtJYSaOBq4Apkg6WNBaYDvzZndtm\nZvXQ6T6LU4DFwGeB9+X/T46IOcAhwFnAs8BOwOEdjs3MzHoxrO4NNW/ePHbeeS+efXZuh6N62ahR\n4umnH8JNMFXHUHX5KQY3Q7kZqi5WtRmqLh3cA+LZZ5/l8cfnsHjxzZXFsMYanwUeqqx8q5OxeYdt\nNvQNq2QBMGrUasDkysqXfLavNbxA9UdXZgOjLtdZmJlZjTlZmJlZKScLMzMr5WRhZmalnCzMzKyU\nk4WZmZVysjAzs1JOFmZmVsrJwszMSjlZmJlZKScLMzMr5WRhZmalnCzMzKyUk4WZmZVysjAzs1LD\n7nkWZlYn1T8AapNNJvHEEw9VGsNw4GRhZoOo6gdAwZNP+iFQA8HNUGZmVsrJwszMSjlZmJlZKScL\nMzMr5WRhZmalnCzMzKyUk4WZmZVysjAzs1JOFmZmVsrJwszMSjlZmJlZKScLMzMr5WRhZmalnCzM\nzKyUk4WZmZVysjAzs1K1ShaSxku6StIiSTMlvafqmMzMrGbJAvgWsATYCHg/8G1J21YbUp11Vx1A\njXRXHUCNdFcdQI10Vx3AsFGbZCFpTeBdwCkR8XxE/A64GvhAtZHVWXfVAdRId9UB1Eh31QHUSHfV\nAQwbtUkWwDbAixHxQGHYHcCUiuIxM7NsTNUBFKwNLGgatgBYp9UZjBkzhiVLnmDcuAMGNLB2LF16\nR2Vlm5kNFkVE1TEAIGkq8NuIWLsw7FPAHhFxYGFYPQI2MxtiIkL9/WydjizuA8ZI2qrQFLUDcFdx\nolVZWDMz65/aHFkASLoUCOBDwI7Az4BdI+KeSgMzMxvh6tTBDXAssCbwFHAJ8BEnCjOz6tXqyMLM\nzOqpbkcWQHtXckv6pKTZkuZJOk/Sap2MdbC1WheSjpB0u6T5kh6R9CVJtVy//dWfK/wlXS9p2Uiu\nC0mvkvQzSQskPSXpi52MdbC1WRdnSHpM0lxJN0jarpOxDiZJx0q6TdISSeeXTNv2frOuX6CWruSW\ntC/wGWAvYBKwFfD5DsbZCa1e1b4G8HFgA+AfgX2AEzsVZIe0dYW/pPeSTuIYjofPrX5HVgN+Bfwa\n2BjYnNTEO5y0WheHAkcBuwHrA7cAF3cuzEH3OHA68L2+Jur3fjMiavUi9Vm8AGxVGHYhcFYP0/4A\nOKPwfi9gdtXLUEVd9PDZTwJXV70MVdUFMA64F9gZeAkYVfUyVFEXpJNFflN1zDWpi88APyy83w5Y\nXPUyDEKdnA6c38f4fu0363hk0c6V3FPyuOJ0G0saP4jxddKqXNW+B02nHQ9x7dbFWaRfnE8OdmAV\naKcu3gQ8LOkXkp7OTS+v60iUndFOXfwQ2ErS1vmI6yjgl4MfYu30a79Zx2TRzpXcawPzm6ZTL9MO\nRf26ql3S0cAbgLMHKa4qtFwXknYCdgW+0YG4qtDOdrE5cBjwdeCVwC+AqyXV6RqrVdFOXcwGfgf8\nFXgOOAQ4YVCjq6d+7TfrmCwWkZoQitYFFrYw7bqk9umeph2K2qkLACQdBJwJ7BcRzw5ibJ3WUl1I\nEnAO8PFIx9jD8SLOdraL50l3RrguIv4eEWeT+rWGy92c26mL6cAbgc2A1YHTgBslrT6oEdZPv/ab\ndUwWy6/kLgxb6Uru7K48rmEq8GREzB3E+DqpnbpA0n7Ad4B3RMTdHYivk1qti3Gko6rLJM0G/puU\nMB6TtFtHIh187WwXf2F4dvA3tFMXO5D6LGZHxLKIuBAYT+q7GEn6t9+sujOmlw6YS0mdMGsC04C5\nwLY9TLcvMIv0K2k8cCNwZtXxV1QXewNzgGlVx1yDuti48NoJWAZMAMZUvQwV1MU2pF+Se5N+HH4S\n+NsIrYtTgZvydiHS4w8WAuOqXoYBqofRpCOms4CLgLHA6B6m69d+s/IF7GWhxwNX5Y38IeCwPHwL\nUvva5oVpPwE8AcwDzgNWqzr+KuoCuAFYmoctzH+vqTr+qraLwmcmMczOhmq3LoCDcoKYl7eTlXak\nQ/nVxndkLKkfa1aui9uBt1Qd/wDWw3TSD6OXCq9Tcz0sXNX9pq/gNjOzUnXsszAzs5pxsjAzs1JO\nFmZmVsrJwszMSjlZmJlZKScLMzMr5WRhZmalnCysJZImFR8ilO9i+oEOlDtdUkeeOSDpSEk39/Oz\ne0p6tI/x35Z0ck/TSrpT0h79KbcfcV4g6VlJt7Qwbb/XeX4I0d6rGq/Vx3C586QBkh4i3crg76S7\nav4ncGxELB6gIpZfwRkR+7cY00zgmIi4YSDK7YBVKavXz0bER3ubNiKW3zJc0nTSsxmOWIU4eiRp\nGumhWJtGxJIWP9b2Ol9VA7TN2ADzkcXwEsDbI2IcsCPpvkin9DRhvjvriDLcHq3aD5OBh9pIFGbL\njfQvz3C4gemsAAAG3UlEQVQkgIiYTXqwy+sAJN2Ynz/8W0nPAa+SNE7S9yTNkvSopNMbSUTSKEln\n5wfm3A+8fYVC0vyOLrz/kKS783Oe75Q0VdJFwESg8fznE/O0b5L0u/wc5D9J2rMwn8mSupWeJX4t\nsGGvC5qbcyT93xzng/lRqo3xF0j6lqRrJC0EuvIyX6T0LOqZjaahglGSvpGfTXx3sSlF0lGFZbxf\n0r+sHFKfsZzWy3LMlLS30uMuTwIOy2X8SdI/Sbq9afoTJF3Vy7xeKelqSc9Iuk/SP+fhRwPfBXbJ\n857ew2dbXueStlR6vvmcXJeXSGq+VfjOku7KsXxP0isK83pHXr65eZv8hzy8P9vMUZIeyNM/oBae\nzW79UPXNr/wa0BuJzQT2zv9vAdwJzMjvbyTdZO21pB8JY0g3X/sW6U6VG5KeSfyhPP1HgLuBTYH1\nSDegW35Dvjy/o/P/7wYeBXbM77cEtijEtFchxk1Jd8fdN7/fJ7/fIL//PfAVYDVgd9KN4C7qZXn3\nBF4sTL8H6WZyW+fxF5DuQPqm/H4s6W6cV5HuUDqJ9CCcD+bxR+b5HU+6g+ehpButrZfHvw2YnP/f\nndTUN7WNWE4rTPtIL+ttenF5gVfk+nlNYdgfgYN6qZObSDfLW410G+qngK7C8t3Ux/bTzjrfKq+7\nMaTnY3QDX21apr8U5vXbwvK/nvQEw514+e6vM8k3s2tnm8nrcT7w6jxuE4bZjRLr8vKRxfDzE0nP\nknYaNwJfKIz7fkTcGxHLSA+sfxvwyYhYEhFzSE9TOzxP+27g6xExKyLmNc2n2THAlyPijwAR8WBE\nFDt7i01e7yfdDffaPO31pLt/7i9pC9IO5NSIeDEibgZ+VrK8AXwuT38TcA1pJ99wdUQ0OnNfJD01\n7l8jYnFEPAz8G2ln1fBkRPy/iHgpIi4nJZO351h/GREP5f9vBq4jJY1WY2lbRCwFLiPVG5KmkJLc\nNc3TStoc2AX4bI7hDtIdRVvt/2h5nUfEAxFxfaQHKj0DfI2UBIu+UZjXmUDjF/+HgHMj4vZILiY9\nR/tNxcUp/N/rNpPHvwT8g6TVI+LJiLinxeW1NjhZDD8HRsT6EfGqiDguIl4ojCvuwCeRfn3OVjo7\nZi5wLrBRHr9p0/QP91HmFsADfYwvmgQcmstslLsb6ZGfmwJzI+L5FsslT19sg384z6ehuAwbkn4J\nP9I0/WaF9483zX/5/CS9TdJ/5WaVuaRkW2wmK4ulvy4CGk1a7wcuj4gXe5huU+DZWPGEhubl60vL\n61zSxpL+Q9JjkuYBl7Byk+FjTfNq1MUk4FNN28Dm9F5XvW4zeVkPAz5K2pZ/Juk1rS2utcPJYvjp\nq+O6eLbOo8ASUvPP+hExPiLWi4jt8/jZpCTQMKmP+T5KapYoK7Mx7UW5zEa560TEl3OZ4yWtUZh+\nYh/l0sv0s3opfw7p6KK4LJNYMUE071gnArNye/uPgS8DG0XEeFKfULG+y2JpxUpnVEXErcBSSbuT\nkkZvpxLPAtaXtFZTDM0JsDftrPOzSM9OmBIR65GSWPO21zyvRl08SnrYTnEbWDsiLsvj29lmiIhf\nRcRbSQ+4+iupb8YGmJPFCBURT5CaUb4maR0lW+rl8/0vB46XtJmk8cBn+5jdecCJknYEkLRVblKC\n1Da9ZWHaS4ADJL01d6iurtRRvWlEPEJqXvi8pNWUTvU8oGRRVJh+d1KT0eW9LPOyPO5MSWtLmkR6\nclxx57uJpOMkjZH0blIfzzWkvoNXAHMiYpmktwFv7W8sfXgSmCytdLbaxcA3gaUR8ftelu8xUp/P\nFySNlbQ9qYmw1etU2lnn65D6ZBZK2gz4dA/THJvntT6p4/6Hefh3gY9I2hlA0lqS9i8kuZa3mXyE\n805Ja5J+CCwiNUvZAHOyGF76ukagp3FHkHaAdwPPAj8i/TqD9IW+FriDtAO/orf5RcSPSW3Sl0pa\nQOpAXj+P/gLwudx8cELeoR1I2nk8TWqeOJGXt8X3kdqunwE+B1zY9yIzm9SJPYu0U/xwRPytj2U+\nHlgMPEjq17kkIi4ojL8F2Jp0FHI6cEhEzIuIRfmzP8p9QocDV7cRS1+Kcf6IlHSe0YpnQV1MOrOt\nbMf/HuBVOYYrSH0oN7YQA7SxzoHPk551Po/Ur9TTtJeSfpDcT3pS35kAEfEHUr/FN3Nd3kfqfG9o\nZ5sZBZxAOnqaQzqxoPmaFhsAflKeDVn59MmLI6KsqWrIk7Q66Rf3jhHRav+Q2YDxkYXZ0PAx4DYn\nCquKb/dhVnNKt78AOKjSQGxEczOUmZmVcjOUmZmVcrIwM7NSThZmZlbKycLMzEo5WZiZWSknCzMz\nK/W/6rx1g5UWJ7QAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11c3d4b70>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# histogram of predicted probabilities\n",
    "\n",
    "# 8 bins\n",
    "plt.hist(y_pred_prob, bins=8)\n",
    "\n",
    "# x-axis limit from 0 to 1\n",
    "plt.xlim(0,1)\n",
    "plt.title('Histogram of predicted probabilities')\n",
    "plt.xlabel('Predicted probability of diabetes')\n",
    "plt.ylabel('Frequency')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- We can see from the third bar\n",
    "    - About 45% of observations have probability from 0.2 to 0.3\n",
    "    - Small number of observations with probability > 0.5\n",
    "    - This is below the threshold of 0.5\n",
    "    - Most would be predicted \"no diabetes\" in this case\n",
    "- Solution\n",
    "    - **Decrease the threshold** for predicting diabetes\n",
    "        - **Increase the sensitivity** of the classifier\n",
    "            - This would increase the number of TP\n",
    "                - More sensitive to positive instances\n",
    "                - Example of metal detector\n",
    "                    - Threshold set to set off alarm for large object but not tiny objects\n",
    "                    - YES: metal, NO: no metal\n",
    "                    - We lower the threshold amount of metal to set it off\n",
    "                    - It is now more sensitive to metal\n",
    "                    - It will then predict YES more often"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/ritchieng/anaconda3/envs/py3k/lib/python3.5/site-packages/sklearn/utils/validation.py:386: DeprecationWarning: Passing 1d arrays as data is deprecated in 0.17 and willraise ValueError in 0.19. Reshape your data either using X.reshape(-1, 1) if your data has a single feature or X.reshape(1, -1) if it contains a single sample.\n",
      "  DeprecationWarning)\n"
     ]
    }
   ],
   "source": [
    "# predict diabetes if the predicted probability is greater than 0.3\n",
    "from sklearn.preprocessing import binarize\n",
    "# it will return 1 for all values above 0.3 and 0 otherwise\n",
    "# results are 2D so we slice out the first column\n",
    "y_pred_class = binarize(y_pred_prob, 0.3)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.36752429,  0.28356344,  0.28895886,  0.4141062 ,  0.15896027,\n",
       "        0.17065156,  0.49889026,  0.51341541,  0.27678612,  0.67189438])"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# print the first 10 predicted probabilities\n",
    "y_pred_prob[0:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.,  0.,  0.,  1.,  0.,  0.,  1.,  1.,  0.,  1.])"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# print the first 10 predicted classes with the lower threshold\n",
    "y_pred_class[0:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[118  12]\n",
      " [ 47  15]]\n"
     ]
    }
   ],
   "source": [
    "# previous confusion matrix (default threshold of 0.5)\n",
    "print(confusion)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[80 50]\n",
      " [16 46]]\n"
     ]
    }
   ],
   "source": [
    "# new confusion matrix (threshold of 0.3)\n",
    "print(metrics.confusion_matrix(y_test, y_pred_class))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- The row totals are the same\n",
    "- The rows represent actual response values\n",
    "    - 130 values top row\n",
    "    - 62 values bottom row\n",
    "- Observations from the left column moving to the right column because we will have more TP and FP "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7419354838709677\n"
     ]
    }
   ],
   "source": [
    "# sensitivity has increased (used to be 0.24)\n",
    "print (46 / float(46 + 16))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.6153846153846154\n"
     ]
    }
   ],
   "source": [
    "# specificity has decreased (used to be 0.91)\n",
    "print(80 / float(80 + 50))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Conclusion:**\n",
    "\n",
    "- **Threshold of 0.5** is used by default (for binary problems) to convert predicted probabilities into class predictions\n",
    "- Threshold can be **adjusted** to increase sensitivity or specificity\n",
    "- Sensitivity and specificity have an **inverse relationship**\n",
    "    - Increasing one would always decrease the other\n",
    "- Adjusting the threshold should be one of the last step you do in the model-building process\n",
    "    - The most important steps are \n",
    "        - Building the models\n",
    "        - Selecting the best model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8. Receiver Operating Characteristic (ROC) Curves \n",
    "\n",
    "**Question:** Wouldn't it be nice if we could see how sensitivity and specificity are affected by various thresholds, without actually changing the threshold?\n",
    "\n",
    "**Answer:** Plot the ROC curve.\n",
    "- Receiver Operating Characteristic (ROC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEgCAYAAABM0P/cAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcXFWZ//HPlx0kkARERJYIgixRA4ILi0ZQUZAR8ScK\nChNhcGNmRAaXQQRFVFRUHBdEBdlkcQFF0cEFGgFBQSagAQEhQEAgQBoSIGzJ8/vjnM69Kaqrb1XX\n1t3f9+tVr666deueU09V31P3POfeo4jAzMysGSv0ugJmZjb2uPEwM7OmufEwM7OmufEwM7OmufEw\nM7OmufEwM7OmufGwcUXShpJ+L+lRSUu6VOYxkm6teXxLk9u4VNJ321+73pC0iaSlknbsYplLJe1f\nejxJ0gWSHpa0RNLG4y3OveTGo49I+kH+B1gq6RlJ8ySdLmmDOutuKuk0SXdLelLSPfnxpnXWXV3S\nUZKul/SYpAclXSXp3yWt3p131zVHAusCLwWe38VyyydMfRl4VRfLXkbSrZKO7kXZdXT7JLL1gZ+U\nHn8QeCWwI+m7MA94G3B4l+s1Lrnx6D9/AJ4HbATsB2wL/Ki8gqRtgWuBDYB3AZsB78yPr5X00tK6\nk4A/AocC3wBeDbwcOAF4B/CGzr6d5UlaucNFbA78OSJuj4j5rW5E0kqtvjYiHo+IBa2+fhxRNwuL\niPkR8VRp0ebAnIi4MT8XEfFwRDw6mnK68B0eGyLCtz65AT8AflOz7N+BJcCapWXXA/8HrFCz7orA\nDcB1pWXfAB4DNh6mzLUa1Oc5wInAXcATwO3AJ/JzmwBLgR1rXnMrcHTp8VLgP4AfAg8D5wBXACfX\nKe8m4NjS43fl97kYmAt8BVijQX2X5lgN/T01L18fOBcYBB4HLgVeXnrda/Nr9gAuz+u8f5gyVgVO\nyu/lIeDbwOeBW0rrHAPcWno8DfgpcE/+LG4A3lOz3UuB7wNfAB4AHgFOBlapWe8/cpwWAzeTjrRW\nKG2jNgYb5+deRPpVPggsAC4Gppe2Oyl//+7Nn/WdwAkjfF+fm19zX67PTcCs4b4fwHHAjTkGd+U4\nrlW1DsDO+buzMN/+D3hDzee/f74/txSHpcAlefkA8N0KMV2x9Pxc4LPAt4AHgat6va/oh1vPK+Bb\n6cOoaTxIRxKXAU8Bq+dlL83/DPsNs4335H+a6aRffg9RZ0ddsT4DwD+AvfIOcEfgoPzcJrmcKo3H\nA8CHgBeSjpIOyfVaubTeK/L2NsuPZ+V19s9l7QzMBk5vUN/1gCuBM/OObVJe/ifgOtJR1zakhmQB\nMDU/P9R43AjsmcvbYJgyvkbaWb4F2ILURfUIz248yo+n5/c/Pcfg0PyZvra0zqUUDcaLcz3uB75S\nWufTeUf2L7mObwLuAD6Tn59CauC/lGOxXv4OrEfaIX8T2Jr0i/zr+XNZJ7/2f0g74+2BDUndbgc3\niPVqpB3utcDrcn1eB7yj9P2obTyOzN+hjfO6NwI/KD0/bB1IP4weyvHeNH+P3grsVPNdG2o81smf\n80D+Lkwuxfm7VWOa15lL+rFwNKkR3rLX+4p+uPW8Ar6VPozUeDwNLCL9Ohv69fjF0jrvyMteNsw2\nts2ve3v+p1kKHNZCXXbL5Ww7zPPNHHnU/tJbm/Tr/u2lZd8Eriw9ngu8r+Z1u+Ttrd2g3rU7h6H3\n8eLSslWAfwJH5cdDjcf+I8RkDdKv04Nqll9Dg8ZjmG39jFKjnut9O6DSskNynFbPt8eAN9Zs5wBg\ncLj4l+rzx5plIv0w+M9SfU5t4vtxcK7b85v5ftSsszewuCYmdesATM6f42sabG+5z5D6R/LLvh9N\nxHQu8Ntm/4fG+63lfl3rmKuBA0lf7H2B1wOfanFbo+lz3o70D/R/o9jGkGvKDyLiEUkXkv5Jf5rz\nC+8EPgkgaV3Szuerkr5SeqlISdgXAX+pWPbWwEMRcXOp/Kck/Yl0FLJscW0969iM1PBcVbP8CtKR\nQl15UMIxpKOV5+dtrELakZX9OfLeKruS1E22GemX/uqkeJVfsyKwiqR1IuKhYaqwA7C9pEU1y1cj\nHYVA6n77qaTtgUuA/wUurqlP2XbAjRFx7zDPP4ukfYAPkz6/tUg511UkrR8R9zWqQ0Q8LOkU4DeS\nLiEdkV8QEU2NaquxDdVj+udRlDMuufHoP4sjYm6+/2lJLyL9Kn9fXnYLaSc6nZT7qDWdtCP8O6lb\nYpC0A223pflvbQNVL5n4WJ1lZwDnS1qHdETxHOC8/NzQQI7/JHU71Lq7qZpWV6+etVppkE8gdf19\nhPT5PQZ8lbQDrVreUEz+H+noolajBP0KwO9I3WW19X8EICJ+I2kjYHdgJnAWcIOk3Ro0IJVJeiVp\n4MfngCNI38tXA6eRGtIR6xAR75N0IvDGfPuspEMj4nstVquZmFb5bkwoHm3V/z4NvFfSdgARcT3w\nN+Cjkpb7/CStCHwUuD4i5uR/+rOBd0vapN7GJQ23A/sLMGWo3DoeyH+XDSOWtB7wgkrvKiVsF5BG\nlB0A/DIihnZk80nDKreMNGqq9vbU8Jt9ljnAOpK2LNVzVdIQzr82sR2A20i5itpzF3Ya4XW7AD+M\niJ9GxF9J3SBb1FlvBy3/E3gnUuL4NtL7eIKUE6oXk6Ed/FOkX85l15J+Zd9T53XLjlYijUQ6LyI+\nSDqSmsnwPzz+Amxdbxj5MHYCHoiIYyLimoj4B2lE4XJGqkOkkVMnRsQewCkUP6paUTWmVocbjz6X\n/8l+QRrRM2QWqVvn15J2ySfG7QL8mpRonFVa95OkX1VXSzpE0kslTZP0NkmXkf4565V7Cak75jxJ\n/5Jfs6Okg/PzT5C6VT6Wt/ly4HTSP2OV97WENPLqg6RRTqfXrPJJ4D8lHSlpG0lbSNpb0neqbL/m\nfVwDnJ3rP5101LMqUN7WiEcUEfF4fs1xkvbKdfoiKcHdyM3AWyXtIGlr4LuUGt2SdYBvSdpS0p7A\nscB3ImJxRDxG+g58XtKHctlbS3qnpONL25gL7CRpo3xUB+nIdUXgQkk75xP4dpZ0nKRXAeT7b8vb\n3Zw08GIRaVRUPeeQRkNdKGm3/P3YVdK+DWLwXEkHSXqhpANJn/0yjeogaTNJx0vaSelkv1eTGuU5\nw5Q3oiZiavX0OuniW3GjToIvL381NclCUj/4D0hdOE+ShoH+AHhhndevDhxF6uZ6jDRq5WrSUcpq\nDerzHNKonHsofgF/rPT8i0j99otIO4e9Sd0y5YT5EoZJRJNGji0hjQRaoc7z/0JqoB4ljXa5jpzk\nblDnS3h2gv55pCOwBfn9X0ppIAApYb6EYUZY1WxrNdIQ08F8+w6pK6bRaKsNSQ37ohzLY4DvkYeP\nlur9feCLpOGgQyOvVq0p/6Ach8fz53gVpWHFpHN4rs3Pl4fqbkQahXY/xdDnM4BN8vNHkYYQL8zv\n61Lg1SPEYj1St9P8XN6NwIH5uU2oGY0HfCZ/1ouAX5LyXOU6DlsH0nDrn5Ias8Wk7/13yCPq6n3X\nqJ8wr/f9GCmmtwNH9nr/0G835eB0haRDSb+KXwKcHREHNVj3I8DHSDu+nwAfjIinu1FPMzNrrNvd\nVveQTrY5pdFKknYnNRxD48c3I/1qMTOzPtDVxiMifhYRF9J4ZAikoaqnRMTfIyVRjwXe2/EKmplZ\nJf2aMN+G5YehXg+sJ2lKj+pjZmYl/dp4rEkef54tJI2GmdSb6piZWVm/niT4KMufQLU26cS32jNk\nkeSx2GZmLYiIlq9C0a9HHnOAl5UezwDuj4jBeiv3eshav9yOOeaYntehX26OhWPhWBQ3qLdsdLra\neEhaUdJqpBOWVpK0aj4rutYZwMGStsp5jqNIY7atgTvuuKPXVegbjkXBsSg4Fu3T7SOPo0gn4nwc\neHe+/8l8NuwiSRsCRMTFpMtKX0o6mek20mU6zMw6bupUkMbPbUoHhhp19STBTpAUY/09tMvAwAAz\nZ87sdTX6gmNRcCwKVWMhwXjfrUgiRpHzcONhZlbDjcfI+jVhbi0YGBjodRX6hmNRcCwKjkX79OtQ\nXTOboKZOhcG64yq7pxM5gvHG3VZm1lcmQpdRP3C3lZmZdZ0bj3HE/bkFx6Iw2lh0e9hqJ7uM/L1o\nH+c8zKyhwUF3I9mzOedhZg05BzE+OedhZmZd58ZjHHF/bsGxKOcqBvo2B9Ft/l60jxsPs3FqKFdx\n6aXpb6u3BSPN+2kTknMeZuOUcxXWiHMeZmbWdW48xhH35xYci4JjUXAs2seNh5mZNc05D7NxyjkP\na2S0OY8RzzCXtC7wBtKc4pOBh4Hrgd9FxAOtFmxm7Ve+Iu14GmJr/WfYbitJW0g6F7gFOASYRGo4\nJuXHN0s6V9KLu1JTG5H7cwsTNRZDw3PLQ2wnaizqcSzap9GRxw+BrwDvjYjFtU9KWg14G3Am8IrO\nVM/MzPqRcx5m44jzHFZVV87zkHSopKmtFmJmZuNL1aG6ewB3SfqZpLdLWqWTlbLWuD+34FgUHIuC\nY9E+lRqPiNgT2BS4FPgEcK+k70jasZOVMzOz/tRSzkPSDOB0YDpwB/Bd4BsR8Xhba1etLs552LhX\nHoLbyJQpvpChVTPanEdTjYek1wLvIY2yuoHUgNwFfBiYHBEzW61Iq9x42ETgRLi1W7cS5sdLugM4\nmXSksV1E7BoRp0fEpcC+wPatVsLaw/25Bcei4FgUHIv2qTqH+WRgv4i4qt6TEfGUpFe1r1pmZtbP\nKnVbSfpqRBxeZ/kJEXFER2pWkbutbLyqvdSIcxnWTl3JeUhaGBFr1Vm+ICJ6ev6HGw8br5znsE7q\naM5D0oGSDgRWknTA0ON8+zTwYKsFW/u5P7fgWBQci4Jj0T4j5TwOyX9XAd5XWh7A/cB7O1EpMzPr\nb1W7rY6PiE90oT5Nc7eVjQVVz9Moc57DOqmr53n0IzceNhY4f2H9pmM5D0kLS/eXSlpSc1sqaUmr\nBVv7uT+34FgUHIuCY9E+jRLmLyvd3xzYouY2tKwySVMkXSDpUUlzJe3XYN3jJN0taVDSJZK2bqYs\ns06ZOjUdSTRz86x+Nt5UzXlMj4i/jbow6Zx89yBgO+Ai4NURcVPNevsCXwV2Il3+5HPA7hHx8jrb\ndLeVdZW7oGw86MrlSYDLJN0g6ROSNm6lIElrAPsAR0XE4oi4Evg5cECd1acBV0TEnbllOAvYqpVy\nzcys/ao2HusDnwJmADdKulzSByWt00RZWwBPR8RtpWXXA9vUWfdcYDNJm0taGZgF/LqJsiYk9+cW\nHIuCY1FwLNqn0rWtIuJp0lHCzyU9h3RV3fcBXwNWq1jWmsDCmmULgUl11r0XuBK4GXgGmAfsWrEc\ns1FrNLTW+Quz6hdGBCDPIPhG4K2knMXVTbz8UaD2EidrA4vqrHsMsAPwAtLJiAcAl0raOiKeqF15\n1qxZTJs2DYDJkyczY8YMZs6cCRS/NCbC45kzZ/ZVfcby48HBmUQM/zz0V31HejykX+rTq8dDy/ql\nPt18PDAwwGmnnQawbH85GlUT5m8E9gf2Bv5B6lY6JyLuqVxQynksALYZ6rqSdAZwd0QcWbPuL4Df\nRMQ3SssGgd0i4rqadZ0wt7ZzUtzGu24lzL8J3Am8MiK2j4gTmmk4APIsg+cDx0paQ9LOwF7AmXVW\nvwZ4h6T1lBxAOkr6RzNlTjS1vzInsqqxGG7Y7XjqmvL3ouBYtE/VnEdT53M0cChwKjCfdFHFD0TE\nTZI2AuYAW0fE3cAXgecCs4E1SI3GPhFRmzMxG5XBQR9hmLVi2G4rSZ+IiOPz/aOH20BEHNuhulXi\nbisbDXdP2UQ12m6rRkcem5Xub95qAWZmNv4Mm/OIiENK9w8Y7tadaloV7s8tOBYFx6LgWLRPpYS5\npPnDLP9ne6tjZmZjQdWhuosiYlLNspWA+yOimbPM2845DxsN5zxsoupkzgNJl5JmDVxN0iU1T29I\ncycJmpnZODFSt9VZwNmkS4T8sHQ7C/gI6TIl1icman9u/XM1Bnyp9Gyifi/qcSzap+GRR0ScAiDp\n6nZckt2sE+qdqzEwAKUrUphZmzU6z2O/iDgn3z9wuA1ExBkdqlslznmY8xZmzevYHOaSLo6I3fP9\ny4d5fUTEa1otvB3ceJgbD7PmdezaVkMNR76/yzC3njYctjz35xYci4JjUXAs2qfqeR5T81VxkbSC\npAMk7S+p5VbLzMzGrqrneVwNfCgirpP0BdIoq6dJl03/rw7XcaS6udtqgnO3lVnzOpbzqClkEJga\nESFpHrAzaXKnv0bEBq0W3g5uPMavRrP5lU2ZAgsWdL4+ZuNJt+bzWAKsLGk6sCgi7iRN7LRmqwVb\n+423/tyhIbgj3eo1HOMtFqPhWBQci/apOg3txaTZA9fNfwG2Js01bmZmE0zVbqvVgfeS8hynRcTT\nkl4HbBARP+xwHUeqm7utxinnMsw6pys5j37mxmP8cuNh1jldyXlImizpWEkXSrqkfGu1YGs/9+cW\nHIuCY1FwLNqnas7jbFJy/MfA452rjpmZjQVVcx4LgfUi4onOV6k57rYav9xtZdY53Rqq+1egp+dz\nmJlZ/6jaePwW+LWkj0k6sHzrZOWsOe7PLTgWBcei4Fi0T9Wcx27AfGCvmuUB9PSS7GZm1n0eqmt9\nyzkPs87pVs4DSVMk7Sfp8Px4fUnOg5iZTUBVz/PYBbgFOBj4TF68JfCdDtXLWuD+3IJjUXAsCo5F\n+1Q98vg68O6IeD3wTF52NfCKjtTKJoSpU1PX1HC3KVN6XUMzG07lS7JHxJR8f0FETJW0AvBARKzT\n6UqOUDfnPMYo5zTMeqdbOY+/S3p9zbJdgb+1WrCZmY1dVRuPI4BzJZ0CrC7pW6Qhuh/rWM2sae7P\nLTgWBcei4Fi0T6XGIyKuBGYAt5EajXuBV0fEnzpYNzMz61MtnechaVJELOpAfZrmnMfY5ZyHWe90\nNOch6d2S3lB6vK2kO4CHJc2RtHmrBZuZ2dg1UrfVx4AHSo+/D/wB2A64AjihQ/WyFrg/t+BYFByL\ngmPRPiNd22pj4AYASRsCLwPeGBEPSfoocGuH62dmZn2oYc5D0gPAhhHxpKR9gaMjYnp+bgXg4YhY\nq3Jh0hTgVOANpCOaIyPinGHWfSHwP8BrgSeAUyPiE3XWc85jjHLOw6x3On2ex+XAZyVtDfw78MvS\nc1sC9zdZ3rdJDcFzgfcAJ0naqnYlSSuTLgP/O2A9YEPgrCbLMjOzDhmp8fgw8CrgL6TLkhxfeu5f\ngd9ULUjSGsA+wFERsTgP//05cECd1WcB90TE1yPiiYh4KiJ8QuII3J9bcCwKjkXBsWifhjmPiJgH\nvGaY5z7eZFlbAE9HxG2lZdeTuqVqvQq4U9KvgB1IMxn+pxsQM7P+MGzOQ9K6EfHgiBuovt7OwI8i\nYoPSsn8D9o+IXWvWvRiYSZp86hLgMOCDwIsj4pmadZ3zGKOc8zDrndHmPBodeVwu6bfAmcC15T20\nJAEvBw4kzTK4TYWyHgVqk+trA/VONlwMXBERQ91iJ0g6CtiKdBSynFmzZjFt2jQAJk+ezIwZM5g5\ncyZQHKb6sR/7sR9P5McDAwOcdtppAMv2l6PR6MhjVdKv/fcDLyBdmmQRMAnYDLgLOBn4bkQ8OWJB\nKeexANhmqOtK0hnA3RFxZM26xwI75kvADy17GNglIv5as66PPLKBgYFlX5p+M3UqDA4uv2zKFFiw\noDPl9XMsus2xKDgWhY4deeQG4UTgxDxs9iXAZGAQuCEi7mymoIh4XNL5wLGSDiGdaLgXsGOd1c8C\nDpe0KzBAStw/ANzUTJnWPwYH3UVlNp50dQ7zmvM8HgQ+HhHnSdoImANsHRF353X3Br5MGtZ7HXBo\nRDyr8fCRx9jg/IZZfxntkUdXG49OcOMxNrjxMOsv3ZoMysaAoeSYORZljkXBsWgfNx5mZta0prut\nJK0XEfM7VJ+mudtqbHC3lVl/6Uq3laS1JZ0haTEwNy/bS9JnWi3YzMzGrqrdVieRLmi4BfBUXvYn\nYL9OVMpa4/7cgmNRcCwKjkX7jDSfx5DXky7N/pSkAIiI+ZKe17mqmZlZv6qU85D0D2DniLhP0oKI\nmJrPzfh9RGzR8Vo2rptzHn2m22eTm1nzujVU91Tgx5J2AVaQtAPwA9LlScyWM3Q2efnmhsNsfKna\neHwB+BlwCrAacDbwv8DXOlQva4H7cwuORcGxKDgW7VM157FORHwF+Ep5oaR1SZcZMTOzCaRqzmNh\nvbnKh/IfHalZRc55dEe9PMZwnN8w63+dnM9juXLqFLwmsLTVgm1s8VVxzaysYc5D0lxJtwOrS7q9\nfAP+CfyiK7W0StyfW3AsCo5FwbFon5GOPP6NdNRxIXBIaXkA90fEnE5VzMzM+lfVnMdaEbGwC/Vp\nmnMe3eFrU5mNL13JeUTEQknTgV2AdSnlQCLi2FYLNzOzsanqhREPBv4M7AF8EtgB+ASwTeeqZs1y\nf27BsSg4FgXHon2qniT4CWCPiNgLWJz/7gs81rGaWVdMnZq6pEa6TZnS65qaWT9p+jwPSQ+R5hUP\n4CGf5zG2OZdhNjF16zyPuyVtEhF3ArcCe5LOLH+61YLNzGzsqtpt9RVger5/HPAj4HLgc52olLXG\n/bkFx6LgWBQci/apOtrqlNL9X0qaAqwaEY90rGZWWTOXDqnlXIaZtaLpOcyXvVBaBzgiIv67vVVq\nuh4TPufhvIWZNavj83lI+ldJX5P0IUkr5fnMvwzcAWzXasFmZjZ2jXRtqy8BxwMbA0eTJoW6FtgI\n2CUidu94Da2u8hDboa4n9+cWHIuCY1FwLNpnpJzHu4DXRMStkrYC5gD7RcR5na+aNeKr3JpZLzXM\nedTO4yHp0YhYsys1q2ii5jyc5zCz0ej0eR6StBHFtayeqXlMRNzVauFmZjY2jZQwfw4pMT50Wwu4\ns/R4bqcqNtGNdNmQekNs3Z9bcCwKjkXBsWifkY48Vu5KLexZnNMws37W8nke/WK85jyc0zCzTur4\neR5mZma13Hj0maFcRyuXDXF/bsGxKDgWBceifapeVde6xLkOMxsLKuc8JK1EmkHwBRHxE0mrA0TE\n4g7Wr0q9xlXOw7kOM+uGruQ8JG0D/B04EzgtL96NdLmSyiRNkXSBpEclzZW0X4XX/F7SUknuYjMz\n6xNVd8gnAcdFxIsoJoAaAHZpsrxvA0+QZiJ8D3BSvuxJXZL2J3Wt+bd4Be7PLTgWBcei4Fi0T9XG\n4yXA6fl+AETEo8AaVQuStAawD3BURCyOiCuBnwMHDLP+WqSLMX60ahlmZtYdVecwnw0cFBHXSVoQ\nEVMlbQ+cFBE7VCpImgFcUb42lqTDgddGxFvrrP9N4BZSA3M7sHJELK2znnMeZmZN6tZ5HkcDF0n6\nFLCKpI8CP8nLq1oTWFizbCEwqXbF3DDtCHyjie2bmVmXVJ2G9kJJ9wKHAFcCLwbeGRF/aqKsR0nX\nxipbG1hUXiBJwLeAD0dE5McNzZo1i2nTpgEwefJkZsyYwcyZM4Gij3OsPIYBBgZae325P7df3k+v\nHg8t65f69PLx7NmzOeyww/qmPr18fOKJJ47p/cNoHg8MDHDaaacBLNtfjkbVbqspEdHiLNnLtrEG\nsADYJiJuy8vOAO6OiCNL660NPATMJ129d0VgXeA+4B05V1LerrutsoGBgWVfmonOsSg4FgXHojDa\nbquqjcfjwO+AHwIXtnpuh6SzSQn3Q0hT2P4C2DEibqpZb73Sw42BPwMbAA9GxDM16475xmPq1HRy\nIKQzyxcs6G19zGz861bO44WkxuMjwH2SzpT0ZkkrNlneoaQRWvOBs4APRMRNkjaStFDShgARMX/o\nBjxAanDm1zYc48XQWeURbjjMbGyo1HhExP0R8T8R8SpgW+Bm4ATgn80UFhGDEfG2iFgzIqYNTWcb\nEfMiYq2IuLvOa+6MiBXrjbSy5ZX7+yc6x6LgWBQci/Zp5azttfNtEvBYe6tjZmZjQdWcxxbAfsD+\npIbjx8A5EfHHzlZvZOMh5+FzO8ys27qVMH8EuAA4B/hdRCxptcB2c+NhZta8biXMnxcRsyLi4n5q\nOGx57s8tOBYFx6LgWLTPsCcJStovIs7JD/cd7ly9iDijExUzM7P+NWy3laSLI2L3fP/yYV4fEfGa\nTlWuCndbmZk1rys5j37mxsPMrHndmgzqmmGWX91qwdZ+7s8tOBYFx6LgWLRP1YT5lsMs36JdFZmI\npk5NRx1TpvS6JmZmzWnYbSVpaJrZd5Oua1U2DVg1InbqTNWqGcvdVu6uMrNeGW231UiXZL9nmPsB\n/AU4r9WCzcxs7Kp6kuCeEXFRF+rTNB95FHy56YJjUXAsCo5FoWNHHpJ2Ks2dsUhS3SG5EfGHVguf\nKMqXXC9zrsPMxqpG53n8PSK2zPfnDfP6iIiNO1W5KsbCkYdzG2bWb3yehxsPM7OmdevaVrWF7iJp\nx1YLtc7wGPaCY1FwLAqORftUPUlwQNLO+f4RwPnATyV9vJOVMzOz/lR1tNVDpCvrPiPpVmBvYBFw\neURs0uE6jlQ3d1uZmTWp0+d5DFkBWCppU2CliJiTC5/aasFmZjZ2Vc15/BE4EfgSaVIockPyUIfq\nNS50+/Ij7s8tOBYFx6LgWLRP1cZjFvAEcDNwTF62NfCNDtRp3BgcTN1VCxb0uiZmZu3lobod5FyH\nmfWrbl2SfSVJn5J0i6TH8t9PSVq51YLNzGzsqtpt9UVgT+Aw4BX575uB4ztUrzFhKKcx3K3blx9x\nf27BsSg4FgXHon2qjrbaF9g2Ih7Mj+fkCaJmA//VkZqNAUM5DTOziabqeR7/BKZHxILSsnWAv0XE\n8ztYvxH1MufhnIaZjVXdujzJT4ALJe0maXNJrycN2f1JqwWbmdnYVbXx+CjwB+AUYA7wPeDKvNz6\nhPtzC45FwbEoOBbtUynnERFPAkfmm5mZTXAjzWG+OeloYzpwHXBQRNzVpbpV4pyHmVnzOp3z+CZp\n7vJZwIOkS5RMeN2+7IiZWb8ZqfF4Oelo40LgEOCVna9S/+vXy464P7fgWBQci4Jj0T4jNR6rRMRi\ngIhYBKzH/BMlAAASvklEQVTe+SqZmVm/Gynn8STwudKi/wa+UF4nIo7tTNWq6VbOY+rUdMQBqbuq\n3446zMya0emcx4+AzUu3n9Q8flEzhUmaIukCSY9Kmitpv2HWO1DStZIekXSXpC9KamnK3HYZ6qrq\nx+4qM7Nua7hDjogDRrgd2GR53yZd2v25wHuAkyRtVWe91YEPA+uQ8iy7AUc0WdaE4/7cgmNRcCwK\njkX7VL221ahJWgPYB9g651GulPRz4ABqzh+JiJNLD++V9ENgZrfqamZmjXVtPg9JM4ArImLN0rLD\ngddGxFtHeO0FwE0R8ayTFDuZ83Cew8zGq27NYd4OawILa5YtBCY1epGkg0hDhg/uUL2G5avmmpnV\n183G41FgrZplawOLhnuBpL1Jo712K1/Rt9asWbOYNm0aAJMnT2bGjBnMnDkTKPo4J8Ljcn9uP9Sn\nl4+HlvVLfXr5ePbs2Rx22GF9U59ePj7xxBMn9P7htNNOA1i2vxyNyt1Wkl4HvAt4XkTsLWk7YFJE\nXFbx9WsAC4BtIuK2vOwM4O5huqPeBJwO7BERf2mw3Y51W421y48MDAws+9JMdI5FwbEoOBaF0XZb\nVZ3P40Ok0U6nAh+NiLUlTQe+ExE7N1HZs4Egna2+HfALYMeIuKlmvV1Jw4T3jogrRthm041HOZfR\niPMcZjZedavxuA14Q0TcLmkwIqZIWhGYHxHrNFHZKaQG6A2ka2V9PCLOk7QR6VLvW0fE3ZIuAXYm\nDesVqcG5PCL2rLPNphuPsXZEYWbWbt2aDGoScGe+P7TbXQl4qpnCImIwIt4WEWtGxLSIOC8vnxcR\na0XE3fnxrhGxSl42Kf99VsNhyyv39090jkXBsSg4Fu1TtfG4gmefpHcoUCnf0S98NVwzs/ao2m31\nAuCXpCOQTYBbSEcde0TEvR2t4ch1q9xt5e4qM7OkKzmPXJCAHYGNgXnAVRGxpNWC28WNh5lZ87qV\n8yCSKyPinIi4oh8aDlue+3MLjkXBsSg4Fu1T6SRBSXMpEuXLiYhN21ojMzPre1VzHrvVLHo+8B/A\nORHR06lp3W1lZta8ruU86hT8fOBXEbFtq4W3gxsPM7PmdS3nUcdiwF1WfcT9uQXHouBYFByL9qma\n8zi6ZtEawJ7Ab9peIzMz63tVcx5n1ix6DJgNnBYRT3SiYlW528rMrHkdn88jX8Pqt8CPet1QmJlZ\nfxgx55HP5/iGG47+5/7cgmNRcCwKjkX7VE2YXyRpj47WxMzMxoyqOY9zgb1JF0icR+mEwYg4qGO1\nq8A5DzOz5nVrDvNbgS+3WoiZmY0vDY88JO0XEed0sT5N85FHwVNsFhyLgmNRcCwKnT5J8ORWN2xm\nZuPXSEceiyJiUhfr0zQfeZiZNa/TOY8VJb2ONI94XRFxSauFd8PUqTA4mO57BkEzs/YY6chjCWnu\n8uEaj+j1JdlHOvKYSEcb7s8tOBYFx6LgWBQ6feTxWK8bBzMz6z8jHXksjIi1ulifpvnIw8yseZ0e\nbdXyhntp6tTUaEjOc5iZdULDxqPfR1oNZ3AwHW1EwIIFva5N9/i6PQXHouBYFByL9hnNZFBmZjZB\ntTwNbb8o5zyGhuVOmTKxjjjMzJrVsznM+0W58XBy3Mysml7OYW59xv25Bcei4FgUHIv2ceNhZmZN\nc7eVmdkE5G4rMzPrOjce44j7cwuORcGxKDgW7ePGw8zMmuach5nZBOSch5mZdV1XGw9JUyRdIOlR\nSXMl7ddg3Y9IulfSw5K+L2nlbtZ1LHJ/bsGxKDgWBceifbp95PFt4AngucB7gJMkbVW7kqTdgY8B\nrwM2ATYDPtPFeo5Js2fP7nUV+oZjUXAsCo5F+3St8ZC0BrAPcFRELI6IK4GfAwfUWf1A4JSI+HtE\nPAIcC7y3W3Udqx5++OFeV6FvOBYFx6LgWLRPN488tgCejojbSsuuB7aps+42+bnyeutJ8uwcZmZ9\noJuNx5rAwpplC4F6c4asCTxSs56GWdeyO+64o9dV6BuORcGxKDgW7dO1obqSZgBXRMSapWX/Bbwm\nIt5as+5s4LiI+El+vA4wH1g3IgZr1vXgXDOzFoxmqO5K7azICG4BVpK0Wanr6mXAnDrrzsnP/SQ/\nngHcX9twwOjevJmZtaZr3VYR8ThwPnCspDUk7QzsBZxZZ/UzgIMlbZXzHEcBP+hWXc3MrLFuD9U9\nFFiD1AV1FvCBiLhJ0kaSFkraECAiLga+BFwKzAVuAz7d5bqamdkwxvzlSczMrPv6/vIkPiu9UDUW\nkg6UdK2kRyTdJemLkvr+s25GM9+L0mt+L2npRI6FpBdK+kU+0p8v6fhu1rXTmozFcZLuljQo6RJJ\nW3ezrp0m6VBJ10h6QtKpI6zb9L5zLPwT+az0QqVYAKsDHwbWAV4J7AYc0a1KdknVWAAgaX/SAJHx\neKhd9X9kZeC3wO+A9YANSd3H40nVWOwLzAJ2AqYCV1M//zqW3QN8Fjil0Uot7zsjom9vpPzIk8Bm\npWWnA5+vs+4PScN7hx6/Dri31++hF7Go89qPAD/v9XvoVSyAtYC/A68AlgAr9Po99CIWwCHAZb2u\nc5/E4mPAuaXHWwOP9/o9dCgunwVObfB8S/vOfj/y8FnphWZiUes11B8SPVY1G4vPk36R3t/pivVA\nM7F4FXCnpF9JeiB31UzvSi27o5lYnAtsJmnzfEQ2C/h156vYl1rad/Z74+Gz0gvNxGIZSQcBLwdO\n6FC9eqFyLCRtD+wIfKML9eqFZr4XGwLvBE4Eng/8Cvi5pG6e79VJzcTiXuBK4GbgMeDtwOEdrV3/\namnf2e+Nx6OkLoeytYFFFdZdm9S/XW/dsaiZWAAgaW/gc8CbImJBB+vWbZViIUnAt4APRzoeH48n\nlDbzvVhMusrDbyLimYg4gZQXGzZXNMY0E4tjgB2AFwCrkS6+eqmk1Tpaw/7U0r6z3xuPZWell5aN\ndFb6kGHPSh+jmokFkt4EnAy8JSJu7EL9uqlqLNYiHXWdJ+le4M+kBuRuSTt1paad18z34gbG54CB\nIc3E4mWknMe9EbE0Ik4HppByHxNNa/vOXidzKiR7ziYldNYAdgYGga3qrLc78E/Sr6gppBMMP9fr\n+vcoFrsCDwI797rOfRCL9Uq37YGlwPrASr1+Dz2IxRakX5m7kn44fgS4dYLG4mjgD/l7IdLUEIuA\ntXr9HtoYixVJR1WfJ121Y1VgxTrrtbTv7PkbrBCAKcAF+Ut/B/DOvHwjUt/chqV1DwPuAx4Gvg+s\n3Ov69yIWwCXAU3nZovz3ol7Xv1ffi9JrNmGcjbZqNhbA3rnBeDh/T561Yx3Ltyb+R1Yl5cH+mWNx\nLfCGXte/zbE4hvRjaUnpdnSOxaLR7jt9hrmZmTWt33MeZmbWh9x4mJlZ09x4mJlZ09x4mJlZ09x4\nmJlZ09x4mJlZ09x4mJlZ09x4TFCSzpR0dK/rMRJJf290KRFJF1eZCGqskbSqpJskrdvrurRL7Wcp\n6Yw8EdMVkl4r6a8VtnGgpItaLH99SXPG0YUge8qNxxgn6Q5Jj+eZ4Rblv+v3qC5nSnoy1+HBvGPf\nfDTbjIgtI+LKvP3P1s6IFhG7R8Q5oymjlqQV84yDQ/G8S9KXmnj9bpLmjrIaHwR+GxEP5m3uKunS\nPDvkLaPcNpJeI+mPeea4ByX9QdKM0W63kZrPciawC7B+ROwcEZdFxEsqbOOMiNgzb2Poc9q4Yvn3\nAZcD/9bym7Bl3HiMfQHsGRFrRcSk/Pe+HtbncxGxFukSCAsYYRazPhbA1vm97AYcIOlfK75WjP4C\nhO9n+ZntHgO+R5rEaFQkTQZ+TrpM/xTSpdqPI13SplumAXMj4slRbKOVOJ9Niq2NkhuP8eFZlxpX\n8uM8L/GCPPHPlnVfLD1X0kW5C+EhSQOl514g6Xyl+a5vk/ShKhWKiMXAOcD0vJ1VJf2PpH9Kmifp\nK0PdByOUPy//St6TtON8dz4iuCY/f3nuylgt/yrfovTa5+Wjsin58b9Imp3L+YOkRhNpaSiuEXEr\n8EfS1UaHtn2wpBvzkcmtkg7Oy9cCLgQ2Lh25rJs/jyMl/SPH8mxJaw/zebyQdKnwa0vx/FNEnE26\nXtNovZg0adL5kTwR6TLtN5be22WSvpWPTObkI4Wh+q0t6dT8Wd4labkpSyW9X6nLbaGkGyS9JC8f\n+iwPAU4CdsnrfLL2aE3Sxkpzkc/Pt6+V6nZpXu2y/Hfoc9gnl7t7aTur5O//0Gd9FbClpOe3IY4T\nmhuP8e0XpPmI1wf+xvBzNH8UuI00t8PzgKNg2XwYvwT+RJo86A3AEZJeN1LBkiYB+wPX5UXHANuR\nGpNtSXNH/3ej8ssi4iLgS8AP8xHWDjXPP0G6IF45//FO4HcRMShpB9Il6g8izVl9KhUnQlKaA3sn\n0gUFh9wHvDkfmRwCfEPS9IhYCOwF3FU6EnyQNNHQm0lXet2QdOG+4SaoeglwW3TuwnM3AyvmBmD3\nYRqxHYEbSZ/JccD5uWGENO/5Y8ALSZe830PSewGU8k9HAvvl2OxDOgJdJiK+B/w7cHmOz+eGnsrb\nWBG4iHSJ9U1IR7E/Km8i/30NqYHfKm/nfNK0sweU1t2LdIQzJ5f9NHA7y1+C3FrgxmN8+Fn+dbVA\n0vkA+RflGRHxeEQ8RZrs5uWSVq/z+qeBDYBpkSYJuiIvfzUwKSK+GBFLIuJ20k73XQ3q8t+SFpDm\nDF+FtLOG1JAcExEL8s70WIp/8uHKb9Y5uZwh+5Muzw1pB//tiLgux+a0vHy5RqjGDZIeJc138Bvg\nu0NPRMRFEXFnvj8A/J7Uhz+c9wNHRsR9+fP4LPCOYdadTAcnMYuIh0mNmEhXUJ2ff+WvU1rtnxHx\nrfy5nwPMBd4saQNSN97hEfFkRDwAfJ3iO3EwcHxEzM5l/SMi7mmyijuSGq1PRMTiXM5VDdYvH3mf\nBexV+p6/h2f/aFpEirGNghuP8eGtETE13/YBkLSCpC/lrqaHSb+aA6g3eucLwF3A73MXzBF5+SbA\nJqWGaZB0lPC8BnX5Qq7HCyJin4i4Ky/fIJcx5E5S1wzA8cOU36zfAWtL2lbSpqT5CS4svZeP17yX\n9Ut1qOclEbEm6Wjm1cBzhp6Q9BZJV+dutkHSUVmjkVEbA78YKp80MdNSSevVWXeQUUyfLOl7pS6z\nurGMiJsi4r0RsRHw0ly/r5ZWubvmJXeSPsNNSJczv78Ux2+S5sWAdJRwG6OzIXBHK0deEXE36Uh5\nn9xd+UZSnqNsEunS4zYKHrI2PtSbXvVA4E3AzIiYl39VPlBv3Yh4lNStcnjuGx6Q9CdgHnBLRDTK\nDVR1D2nHM9T1s0leRkQsqld+RFxeW9VGBUTEEkk/Jh1xPAJcmHMv5PfymYj4chN1Hsp5nCfpbaTu\ntI8qTVX6Y2Bf0jwpSyX9giK29eo5D9g/Iq6pUO4NwGaS1OIO9BDSkVbV9W+WdAbpOzNkw5rVNibN\nfTEPeCwipg6zuXmkrtLRmEf60TLS+x/uuTNIR7WTgMsiYv7QE5JWBjYFrh9lHSc8H3mMX5OAJ4FB\nSc8hzSZW958t/4reND9cBDxDmkTmKuApSYcrJbxXlDRd0nYt1Odc4GhJ60h6LmlHfGaD8pfU2cb9\npFE6jZxDynXsx/K/OL8HHCpp+1zmmrncet149RwPfCA3wqsCK5NmawxJbyF15ZTrua6kNUvLTga+\nIGmjXP56kvaqV1DuDruLlE8gry9Jq5K6AlfIn0dLP/4kbSXpI7kLCqWhru8ifd5Dni/pQ/kzfxdp\nh/u/+Zf9ZUoDHiblem0maajL7vvAx5SH/Up6kaRGR3f1XAU8BHxe0upKgyF2rF0pIpaSPoNNa546\nH3glcCipISl7FXBzRNzbZJ2shhuPsW+4X18/AO4l/Vr8K9Aoj/Bi4BJJi0jj4E+MiCsjYgmwB/AK\n0iif+cB3GL5LpdGvxM+Qfu39DZhN2kEc36D8P9bZ5nnAqrm75Op6ZebXPUPqQvpNafmfSOdOnFTK\nyby7QX1rtzubNOLqiIh4hHSk9DPSTm4f0uCEoXXnAD8F7sh1XZfUJfRrUtfcI6TPY/sG5Z/M8kcC\nuwKLc5kvBB4HftXg9Y0sInXDXZNjfgXwF+DjpXX+CGxDSnYfDeyT3zekPMJzSAn1BaRk9vMAIuJc\n4IukeeMfIcVhSn5dpaOo/L17C2k+8XmkLrO3D7P6McA5Oc5759c/TorTxvlv2btJ32EbJc8kaNaH\n8lHGdcBr8wCDbpZ9MPDuiNi1m+W2k9Lw4Y0i4qDSsvVJebEZEfFMzyo3TjjnYdaHIp08145c04ST\nuxbfS8pJLRPp5NnpPanUOORuKzMbNyR9gNTNdUFEXD3S+tY6d1uZmVnTfORhZmZNc+NhZmZNc+Nh\nZmZNc+NhZmZNc+NhZmZNc+NhZmZN+/9kU7IazikDdgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11c3eb128>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# IMPORTANT: first argument is true values, second argument is predicted probabilities\n",
    "\n",
    "# we pass y_test and y_pred_prob\n",
    "# we do not use y_pred_class, because it will give incorrect results without generating an error\n",
    "# roc_curve returns 3 objects fpr, tpr, thresholds\n",
    "# fpr: false positive rate\n",
    "# tpr: true positive rate\n",
    "fpr, tpr, thresholds = metrics.roc_curve(y_test, y_pred_prob)\n",
    "\n",
    "plt.plot(fpr, tpr)\n",
    "plt.xlim([0.0, 1.0])\n",
    "plt.ylim([0.0, 1.0])\n",
    "plt.rcParams['font.size'] = 12\n",
    "plt.title('ROC curve for diabetes classifier')\n",
    "plt.xlabel('False Positive Rate (1 - Specificity)')\n",
    "plt.ylabel('True Positive Rate (Sensitivity)')\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- ROC curve can help you to **choose a threshold** that balances sensitivity and specificity in a way that makes sense for your particular context\n",
    "- You can't actually **see the thresholds** used to generate the curve on the ROC curve itself"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# define a function that accepts a threshold and prints sensitivity and specificity\n",
    "def evaluate_threshold(threshold):\n",
    "    print('Sensitivity:', tpr[thresholds > threshold][-1])\n",
    "    print('Specificity:', 1 - fpr[thresholds > threshold][-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sensitivity: 0.241935483871\n",
      "Specificity: 0.907692307692\n"
     ]
    }
   ],
   "source": [
    "evaluate_threshold(0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sensitivity: 0.725806451613\n",
      "Specificity: 0.615384615385\n"
     ]
    }
   ],
   "source": [
    "evaluate_threshold(0.3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 9. AUC\n",
    "AUC is the **percentage** of the ROC plot that is **underneath the curve**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.724565756824\n"
     ]
    }
   ],
   "source": [
    "# IMPORTANT: first argument is true values, second argument is predicted probabilities\n",
    "print(metrics.roc_auc_score(y_test, y_pred_prob))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- AUC is useful as a **single number summary** of classifier performance\n",
    "- Higher value = better classifier\n",
    "- If you randomly chose one positive and one negative observation, AUC represents the likelihood that your classifier will assign a **higher predicted probability** to the positive observation\n",
    "- AUC is useful even when there is **high class imbalance** (unlike classification accuracy)\n",
    "    - Fraud case\n",
    "        - Null accuracy almost 99%\n",
    "        - AUC is useful here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.73782336182336183"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# calculate cross-validated AUC\n",
    "from sklearn.cross_validation import cross_val_score\n",
    "cross_val_score(logreg, X, y, cv=10, scoring='roc_auc').mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use both of these whenever possible\n",
    "\n",
    "1. **Confusion matrix advantages:**\n",
    "\n",
    "    - Allows you to calculate a **variety of metrics**\n",
    "    - Useful for **multi-class problems** (more than two response classes)\n",
    "\n",
    "2. **ROC/AUC advantages:**\n",
    "\n",
    "    - Does not require you to **set a classification threshold**\n",
    "    - Still useful when there is **high class imbalance**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 10. Confusion Matrix Resources\n",
    "\n",
    "- Blog post: [Simple guide to confusion matrix terminology](http://www.dataschool.io/simple-guide-to-confusion-matrix-terminology/) by me\n",
    "- Videos: [Intuitive sensitivity and specificity](https://www.youtube.com/watch?v=U4_3fditnWg) (9 minutes) and [The tradeoff between sensitivity and specificity](https://www.youtube.com/watch?v=vtYDyGGeQyo) (13 minutes) by Rahul Patwari\n",
    "- Notebook: [How to calculate \"expected value\"](https://github.com/podopie/DAT18NYC/blob/master/classes/13-expected_value_cost_benefit_analysis.ipynb) from a confusion matrix by treating it as a cost-benefit matrix (by Ed Podojil)\n",
    "- Graphic: How [classification threshold](https://media.amazonwebservices.com/blog/2015/ml_adjust_model_1.png) affects different evaluation metrics (from a [blog post](https://aws.amazon.com/blogs/aws/amazon-machine-learning-make-data-driven-decisions-at-scale/) about Amazon Machine Learning)\n",
    "\n",
    "\n",
    "### 11. ROC and AUC Resources\n",
    "\n",
    "- Lesson notes: [ROC Curves](http://ebp.uga.edu/courses/Chapter%204%20-%20Diagnosis%20I/8%20-%20ROC%20curves.html) (from the University of Georgia)\n",
    "- Video: [ROC Curves and Area Under the Curve](https://www.youtube.com/watch?v=OAl6eAyP-yo) (14 minutes) by me, including [transcript and screenshots](http://www.dataschool.io/roc-curves-and-auc-explained/) and a [visualization](http://www.navan.name/roc/)\n",
    "- Video: [ROC Curves](https://www.youtube.com/watch?v=21Igj5Pr6u4) (12 minutes) by Rahul Patwari\n",
    "- Paper: [An introduction to ROC analysis](http://people.inf.elte.hu/kiss/13dwhdm/roc.pdf) by Tom Fawcett\n",
    "- Usage examples: [Comparing different feature sets](http://research.microsoft.com/pubs/205472/aisec10-leontjeva.pdf) for detecting fraudulent Skype users, and [comparing different classifiers](http://www.cse.ust.hk/nevinZhangGroup/readings/yi/Bradley_PR97.pdf) on a number of popular datasets\n",
    "\n",
    "### 12. Other Resources\n",
    "\n",
    "- scikit-learn documentation: [Model evaluation](http://scikit-learn.org/stable/modules/model_evaluation.html)\n",
    "- Guide: [Comparing model evaluation procedures and metrics](https://github.com/justmarkham/DAT8/blob/master/other/model_evaluation_comparison.md) by me\n",
    "- Video: [Counterfactual evaluation of machine learning models](https://www.youtube.com/watch?v=QWCSxAKR-h0) (45 minutes) about how Stripe evaluates its fraud detection model, including [slides](http://www.slideshare.net/MichaelManapat/counterfactual-evaluation-of-machine-learning-models)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [py3k]",
   "language": "python",
   "name": "Python [py3k]"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}