chierici_practical_part3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "TJfNUgelLKKI"
   },
   "source": [
    "# <center>Machine learning from scratch - Part III</center>\n",
    "## <center>WebValley ReImagined 2020</center>\n",
    "### <center>Marco Chierici</center>\n",
    "#### <center>FBK/MPBA</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "t1GLtJ7YLKKJ"
   },
   "source": [
    "Recap. We are using a subset of the SEQC neuroblastoma data set [Zhang et al, Genome Biology, 2015] consisting of 272 samples (136 training, 136 test). The data was preprocessed a bit to facilitate the progress of the tutorial."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "6Dhu_6_ULKKK"
   },
   "source": [
    "We start by loading the modules we need to process the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "jT_SUG13LKKL"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pylab as plt ## for plotting\n",
    "import pandas as pd ## for reading text files and manipulating data frames\n",
    "from sklearn import metrics\n",
    "from sklearn import neighbors ## kNN classifier\n",
    "from pathlib import Path ## for creating paths in a neat way\n",
    "%matplotlib inline\n",
    "np.random.seed(42) ## set random seed just in case"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Data loading"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "qOlPPJmmLKKO"
   },
   "source": [
    "Let's start from scratch by reloading the Neuroblastoma data, preparing them for a classification task on the \"CLASS\" label (see part2 notebook):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##  for convenience, define the data directory as a variable\n",
    "DATA_DIR = Path(\"data\")\n",
    "DATA = DATA_DIR / \"MAV-G_272.txt\"\n",
    "LABS = DATA_DIR / \"labels.txt\"\n",
    "data = pd.read_csv(DATA, sep = \"\\t\")\n",
    "\n",
    "# We drop the first column from the train and test expression sets, since it's just the sample IDs...\n",
    "data = data.drop('sampleID', axis=1)\n",
    "# ...and store the data into Numpy arrays.\n",
    "X = data.values\n",
    "# Now we read in the files containing labels and select the column with the CLASS target\n",
    "labs = pd.read_csv(LABS, sep = \"\\t\")\n",
    "class_lab = labs[['CLASS']]\n",
    "y = class_lab.values.ravel()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# following our consolidated practice, print the shape of the data as a sanity check:\n",
    "X.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "13LEke6Js720"
   },
   "source": [
    "# 2. Data partitioning"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "-IIhs_0Exo1B"
   },
   "source": [
    "### Hold-out strategy\n",
    "\n",
    "The idea behind data partitioning is to split your original data set into a **train** portion (for developing your machine learning model) and a **test** portion (for evaluating the performance of the trained model).\n",
    "\n",
    "The simplest and most straightforward way to partition your data set is to randomly split it in two groups (*hold-out strategy*).\n",
    "\n",
    "You achieve this using scikit-learn's function `train_test_split`, in the `model_selection` submodule.\n",
    "\n",
    "For example, let's split the data (X) into 80% train and 20% test (note the argument `test_size=0.2`). Use a random_state of your choice:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "tqbLHGVP2cAH"
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "x_tr, x_ts, y_tr, y_ts = train_test_split(X, y, test_size=0.2, random_state=78)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Data preprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import preprocessing\n",
    "## first you need to create a \"scaler\" object\n",
    "scaler = preprocessing.MinMaxScaler(feature_range=(-1,1))\n",
    "## then you actually scale data by fitting the scaler object on the data\n",
    "scaler.fit(x_tr)\n",
    "x_tr = scaler.transform(x_tr)\n",
    "x_ts = scaler.transform(x_ts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "FErPMN8F6sb9"
   },
   "source": [
    "# 4. Supervised Learning"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scikit-learn provides you access to several models via a very convenient _fit_ and _predict_ interface."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "T7xRxdjO6vTY"
   },
   "source": [
    "## 4.1 k-NN classifier"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "LnHfIehl62Js"
   },
   "source": [
    "Let's fit again a k-NN model on the whole training data and then use it to predict the labels of the test data. This time we'll use a different number of neighbors, k:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "6P24VUwQk3XG"
   },
   "outputs": [],
   "source": [
    "knn = neighbors.KNeighborsClassifier(n_neighbors=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "Cg8TpDATk3XI",
    "outputId": "e9658389-474c-4bf5-f196-11d4518311b7"
   },
   "outputs": [],
   "source": [
    "knn.fit(x_tr, y_tr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "B5pHIVW2k3XM"
   },
   "outputs": [],
   "source": [
    "y_pred_knn = knn.predict(x_ts) # predict labels on test data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Crr_2DwX7f1Z"
   },
   "source": [
    "_In general, a classifier has **parameters** that need to be tuned. Default choices are not good in all situations._\n",
    "\n",
    "_For example, in k-NN the main parameter is the **number of neighbors** used in the nearest neighbors algorithm._"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "L88ilJJG-Cnf"
   },
   "source": [
    "## Performance assessment: Confusion matrix\n",
    "\n",
    "Confusion matrix recap:\n",
    "\n",
    "|      |  |  Predicted  |    |\n",
    "|------|-----------|----|----|\n",
    "|      |           | 0 | 1  |\n",
    "| True | 0        | TN | FP |\n",
    "|      | 1         | FN | TP |\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "JISD2EVQ9Q9Z"
   },
   "outputs": [],
   "source": [
    "conf = metrics.confusion_matrix(y_ts, y_pred_knn)\n",
    "conf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Kh9MqHB5cC43"
   },
   "source": [
    "The total number of class 0 test samples (AN = All Negatives) should be equal to the sum of the first row of the confusion matrix, i.e., TN + FP:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "pZVN8GKKdOhy"
   },
   "outputs": [],
   "source": [
    "np.sum(y_ts == 0) # total number of \"class 0\" samples in the test set"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "SoBqoaTGcDVy"
   },
   "source": [
    "Similarly for class 1, i.e., AP = All Positives = TP + FN:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "1PVj7JbxdVk0"
   },
   "outputs": [],
   "source": [
    "np.sum(y_ts==1) # total number of \"class 1\" samples in the test set"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "kva3wkz5ddap"
   },
   "source": [
    "To compute the metrics, we'll take the scikit-learn shortcut instead of recomputing them by hand. I'm also throwing in a couple more metrics: precision and F1-score."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"MCC = {metrics.matthews_corrcoef(y_ts, y_pred_knn):.3f}\")\n",
    "print(f\"Accuracy = {metrics.accuracy_score(y_ts, y_pred_knn):.3f}\")\n",
    "print(f\"Sensitivity = {metrics.recall_score(y_ts, y_pred_knn):.3f}\")\n",
    "print(f\"Precision = {metrics.precision_score(y_ts, y_pred_knn):.3f}\")\n",
    "print(f\"F1-score = {metrics.f1_score(y_ts, y_pred_knn):.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "arErpTdvLKKh"
   },
   "source": [
    "So far we focused on the k-NN classifiers. However, we already explored theoretical aspects related to two other broadly used classifiers: Support Vector Machines (SVMs) and Random Forests (RFs). In this part of tutorial, the first thing we want to do is assessing how these two alternative classification methods perform on our neuroblastoma dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.2 Support Vector Machines"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "zMP7DMy3LKKh"
   },
   "source": [
    "We start with SVM. We first rescale the data, import the relevant model and create an instance of the SVM classifier. Use the same seed used previously."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# reminder: x_tr, x_ts were previously rescaled using the MinMaxScaler.\n",
    "# Now we want to standardize the features, so we need to recreate the original arrays:\n",
    "x_tr, x_ts, y_tr, y_ts = train_test_split(X, y, test_size=0.2, random_state=78)\n",
    "\n",
    "scaler = preprocessing.StandardScaler()\n",
    "## then you actually scale data by fitting the scaler object on the data\n",
    "scaler.fit(x_tr)\n",
    "x_tr = scaler.transform(x_tr)\n",
    "x_ts = scaler.transform(x_ts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "oY93_AikLKKj"
   },
   "outputs": [],
   "source": [
    "## import support vector classifier (SVC) and create an instance\n",
    "from sklearn.svm import SVC\n",
    "svc = SVC(kernel='linear')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Pseq_nh7LKKl"
   },
   "source": [
    "Note that the specification _kernel = 'linear'_ implies that a linear kernel will be used. If you remember from the lecture, this means that a linear function is used to define the decision boundaries of the classifier. Alternatives include _‘poly’_ and _‘rbf’_ for polynomial or gaussian kernels respectively. Scikit-learn offers an alternative implementation of linear SVMs. You can find more details in Scikit User Guide. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Tip:** always start experimenting with the linear kernel, which is the simplest one; try more complex kernels later on (Occam's razor...)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "waXcgaLpLKKm"
   },
   "source": [
    "As previously done with the k-NN classifier, we fit an SVM model on the training data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "Qqc3TmFBLKKn",
    "outputId": "d9ef6c64-9f18-4bea-9167-decaa0ca1820"
   },
   "outputs": [],
   "source": [
    "svc.fit(x_tr, y_tr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simple as that. Now we have a fitted SVM model that we can use to make predictions on new data.\n",
    "\n",
    "Of course there are a few parameters that may require tuning: more on this later."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Making predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "Qqc3TmFBLKKn",
    "outputId": "d9ef6c64-9f18-4bea-9167-decaa0ca1820"
   },
   "outputs": [],
   "source": [
    "y_pred_svm = svc.predict(x_ts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now have a set of predictions for each entry in the test set. Since we also have the actual labels for each record in the test set, we can use them to assess the performance of the SVM model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "lI7vGjulLKKr"
   },
   "source": [
    "Now we give a look at the classification metrics introduced in the first part of the tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "Ku0JSF_ALKKs",
    "outputId": "94585c0e-534a-445d-d0ba-92a9bf3a9388"
   },
   "outputs": [],
   "source": [
    "print(f\"MCC = {metrics.matthews_corrcoef(y_ts, y_pred_svm):.3f}\")\n",
    "print(f\"ACC = {metrics.accuracy_score(y_ts, y_pred_svm):.3f}\")\n",
    "print(f\"SENS = {metrics.recall_score(y_ts, y_pred_svm):.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.3 Random Forest"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thanks to scikit-learn simple API, we can easily create a RandomForest instance and fit it to our training data. Note that the Random Forest is robust to feature rescaling (and also supports categorical features), so it is not necessary to apply `MinMaxScaler` or `StandardScaler`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.ensemble import RandomForestClassifier as RFC\n",
    "# x_tr, x_ts, y_tr, y_ts = train_test_split(X, y, test_size=0.2, random_state=78)\n",
    "\n",
    "clf = RFC(n_estimators=500, random_state=42)\n",
    "clf.fit(x_tr, y_tr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Making predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred_rfc = clf.predict(x_ts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute a few metrics using the actual test set labels:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "ZT6XjB20LKK0",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print(f\"MCC = {metrics.matthews_corrcoef(y_ts, y_pred_rfc):.3f}\")\n",
    "print(f\"ACC = {metrics.accuracy_score(y_ts, y_pred_rfc):.3f}\")\n",
    "print(f\"SENS = {metrics.recall_score(y_ts, y_pred_rfc):.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How did you do? Note that (not surprisingly) there is an element of randomness in a \"random forest\". \n",
    "\n",
    "I'm getting an accuracy of about 94-98%, but \"your mileage may vary\": the precise number will be different each time you run the algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Probabilities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Depending on the algorithm, some of the supervised models in scikit-learn can also provide the **probability** that a record belongs in each category. \n",
    "\n",
    "For random forest, we can obtain these probabilities by calling the `predict_proba` method on the fitted model. Because we have two categories, 0 (negative) and 1 (positive), scikit-learn will return two category probabilities (the sum across all categories will add up to 1).\n",
    "\n",
    "Looking at the prediction probabilities may be useful to understand on which samples the classifier is \"unsure\" (i.e., the probabilities are around 0.5 for both classes)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prob_rfc = clf.predict_proba(x_ts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(prob_rfc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Trying different classifiers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "scikit-learn makes it easy to try different classifiers. \n",
    "\n",
    "Neural networks, naive bayes, random forest, logistic regression, support vector machines, and other algorithms all have a very similar interface in sklearn (of course the underlying mathematics can be dramatically different!) \n",
    "\n",
    "1. you create an instance of the model (in the variable, say, `clf`), optionally tuning the parameters; \n",
    "2. you fit the model on some training data (`clf.fit(x_tr, y_tr)`); \n",
    "3. you predict the labels of unseen test data (`clf.predict(x_ts)`);\n",
    "4. you evaluate the performance.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: neural net!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To use a neural net, you only need to make a few changes to this workbook.\n",
    "\n",
    "First, load the appropriate library\n",
    "\n",
    "from sklearn.neural_network import MLPClassifier\n",
    "\n",
    "Then, call a neural network classifier rather than a random forest model\n",
    "\n",
    "clf = MLPClassifier()\n",
    "\n",
    "Proceed as usual to fit and make predictions.\n",
    "\n",
    "_Et voilà!_ In seconds you are making classification predictions using a neural net model rather than a random forest!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Train a neural net and assess its performance on the test set;\n",
    "* Compare the performance of the classifiers you just used;\n",
    "* Given the performance metrics you assessed, can you say there is a \"best\" classification algorithm? Do you think this algorithm will perform better on a different task / data?\n",
    "\n",
    "Hint: the effectiveness depends on a number of factors, mainly the **classifier parameters** and the **data**. Some data sets just match very well to a particular prediction algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "HtPbufFSpSxL"
   },
   "source": [
    "\n",
    "* Compare the metrics of the different classifiers. What can you say about this classification task? Do the classifiers learn something?\n",
    "* Which classifier has the best accuracy and MCC?\n",
    "* Are you confident enough that such classifier is able to *generalize* beyond its training set?\n",
    "* Do you know if there is a more robust way to assess the performance of the models?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "rqkAcvnpxpCx"
   },
   "source": [
    "# 5. Cross-validation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "NqrxSWmp8Zg5"
   },
   "source": [
    "So far we fitted several models on our training data, then we predicted the labels on the test data, computing a set of metrics.\n",
    "\n",
    "How can we know if our model is going to generalize well on new unseen data?\n",
    "\n",
    "We need a more robust way to _estimate_ the model performance and its generalization capabilities.\n",
    "\n",
    "We already used the **hold-out strategy**, with scikit's `train_test_split` function, to split our `X` data into one training/test partition.\n",
    "\n",
    "Partitioning the dataset once is not enough. The partitions depend on the random seed used in the splitting function.\n",
    "\n",
    "More robust strategies involve splitting the data in **multiple (complementary) subsets**.\n",
    "\n",
    "One of such strategies is the **k-fold cross-validation**, introduced during the lecture:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "jnBd81Sk35QM"
   },
   "source": [
    "![k-fold cv](https://www.researchgate.net/profile/B_Aksasse/publication/326866871/figure/fig2/AS:669601385947145@1536656819574/K-fold-cross-validation-In-addition-we-outline-an-overview-of-the-different-metrics-used_W640.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "rmolsnknM8cS"
   },
   "source": [
    "Example of a 5-fold cross-validation (CV) with scikit-learn:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "551-qJzL8WKh"
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import StratifiedKFold\n",
    "skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "rMmmyxp5NX4P"
   },
   "source": [
    "A \"stratified\" 5-fold CV means that the folds are made by preserving the percentage of samples for each class.\n",
    "\n",
    "_Recap:_ the same random_state will generate the same splits. This is useful for reproducibility."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "H_VtvZEINnE6"
   },
   "source": [
    "To actually get the splitting indices and create the folds, we need to iterate over the `skf` object. Note that here I am using a Random Forest: feel free to experiment with other classifiers!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "-uoahY6yNcIv"
   },
   "outputs": [],
   "source": [
    "## get the number of splitting operations\n",
    "N = skf.get_n_splits(X, y)\n",
    "\n",
    "## reinitialize a classifier\n",
    "clf = RFC(n_estimators=500)\n",
    "\n",
    "## create empty lists to store the CV metrics\n",
    "acc_list = []\n",
    "mcc_list = []\n",
    "\n",
    "## split data and iterate over the splits,\n",
    "## computing classifier accuracy & MCC on each test partition\n",
    "n_fold = 1\n",
    "for (idx_tr, idx_ts) in skf.split(X, y):\n",
    "    print(\"### Fold \", n_fold)\n",
    "    X_train, Y_train = X[idx_tr], y[idx_tr]\n",
    "    X_test, Y_test = X[idx_ts], y[idx_ts]\n",
    "    print()\n",
    "    clf.fit(X_train, Y_train)\n",
    "    Y_test_pred = clf.predict(X_test)\n",
    "    acc = metrics.accuracy_score(Y_test, Y_test_pred)\n",
    "    mcc = metrics.matthews_corrcoef(Y_test, Y_test_pred)\n",
    "    print(f\"Accuracy on TEST set: {acc:.3f}\")\n",
    "    print(f\"MCC on TEST set: {mcc:.3f}\")\n",
    "    print()\n",
    "    ## append values to lists\n",
    "    acc_list.append(acc)\n",
    "    mcc_list.append(mcc)\n",
    "    \n",
    "    n_fold = n_fold + 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Q: how are the computed metrics on the different folds?*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get an estimate of the predictive performance of our model, we can average over the cross-validation metrics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## note: we need to convert the lists to numpy arrays before computing the means\n",
    "acc_avg = np.mean(np.array(acc_list))\n",
    "mcc_avg = np.mean(np.array(mcc_list))\n",
    "\n",
    "print(f\"Average cross-validation accuracy: {acc_avg:.3f}\")\n",
    "print(f\"Average cross-validation MCC: {mcc_avg:.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Remarks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want to be cool and have an even better estimate of how the model can generalize on new data, you can **repeat the cross-validation several times** (\"iterated cross-validation\"). Each time you need using a different random_state for generating the splits.\n",
    "\n",
    "Moreover, the average alone is not sufficient as you should also capture the **dispersion** of the values around the average. This can be done, for example, by assessing the standard deviation or, even better, by computing **confidence intervals**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "tI_6v1jNLKLF"
   },
   "source": [
    "# 6. Feature ranking"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One type of insight you can gain from a machine learning model is the feature importance (or weight). Which genes are most likely to influence the classification of our neuroblastoma patients? Are some features pivotal, and others largely ignored?\n",
    "\n",
    "Because a Random Forest model branches repeatedly on different features, the model becomes \"aware\" of which features are particularly influential in classifiying a patient. \n",
    "\n",
    "Scikit-learn allows us to read this information off of a trained Random Forest model through the `feature_importances_` attribute (mind the trailing underscore!)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "2lZAaTXJLKLH",
    "outputId": "2155231c-e50c-4c06-82c4-6b6a5f7c4ee2"
   },
   "outputs": [],
   "source": [
    "# Retrain a random forest\n",
    "rf = RFC(n_estimators=500)\n",
    "rf.fit(x_tr, y_tr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "7gH9yIYOLKLJ"
   },
   "source": [
    "For the sake of completeness make the predictions and check the classification performance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "rspvHmO0LKLK",
    "outputId": "7b131d8f-ebc8-4d03-9f38-ad90de735367"
   },
   "outputs": [],
   "source": [
    "y_pred_rfc = rf.predict(x_ts)\n",
    "print(f\"MCC = {metrics.matthews_corrcoef(y_ts, y_pred_rfc):.3f}\")\n",
    "print(f\"ACC = {metrics.accuracy_score(y_ts, y_pred_rfc):.3f}\")\n",
    "print(f\"SENS = {metrics.recall_score(y_ts, y_pred_rfc):.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "m3k78HePLKLT"
   },
   "source": [
    "Now extract the feature importances and display the first 10:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "7g9k5EHsLKLU",
    "outputId": "aa26094b-0e4a-48f0-be91-ecd2874ab204"
   },
   "outputs": [],
   "source": [
    "# get the importances\n",
    "importances = rf.feature_importances_\n",
    "# sort by decreasing importance\n",
    "indices = np.argsort(importances)[::-1]\n",
    "# get the gene names\n",
    "genes = data.columns.values\n",
    "# print the feature ranking\n",
    "print(\"Feature ranking (top 10 features):\")\n",
    "for f in range(10):\n",
    "    print(genes[indices[f]], \" - \", importances[indices[f]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A **stem plot** is a common way to visually represent this kind of information:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_feat = 10\n",
    "plt.figure()\n",
    "plt.title(\"Feature importances\")\n",
    "plt.stem(range(n_feat), importances[indices[:n_feat]])\n",
    "plt.xticks(range(n_feat), genes[indices[:n_feat]], rotation=\"vertical\")\n",
    "plt.xlim([-1, n_feat])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 7. Further topics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "6khMtHGSLKK3"
   },
   "source": [
    "## Parameter tuning"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "XlHdQ4xSLKK4"
   },
   "source": [
    "As mentioned in the lecture, Scikit learn offers a very useful and flexible tool for parameter tuning called _GridSearchCV_. While the tool is very sophisticated and efficient, it is useful to at least try an example _by hand_ to understand what is happening in the background.\n",
    "\n",
    "For this example we use a linear SVM and try to tune the C parameter. You might remember from the lectures that the paramenter C essentially controls how much we want to avoid misclassifying each training example. Large values of C result in smaller margins, i.e. closer fitting to the training data. As mentioned in the classes, the drawback is over-fitting, resulting in poor generalization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# first of all, let's get a clean train/test split of the original data\n",
    "x_tr, x_ts, y_tr, y_ts = train_test_split(X, y, test_size=0.2, random_state=78)\n",
    "# rescale\n",
    "scaler = preprocessing.StandardScaler()\n",
    "scaler.fit(x_tr)\n",
    "x_tr = scaler.transform(x_tr)\n",
    "x_ts = scaler.transform(x_ts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "XG01P5fdLKK6",
    "outputId": "099e6404-c7fd-414a-b49a-4092af095c57",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## define the sequence of C values we want to use in the search of the best one\n",
    "C_list = [0.000001, 0.00001, 0.0001, 0.001, 0.01, 0.1]\n",
    "for C in C_list:\n",
    "    print(f'C = {C}')\n",
    "    svc = SVC(kernel='linear', C=C)\n",
    "    svc.fit(x_tr, y_tr)\n",
    "    class_pred_ts = svc.predict(x_ts)\n",
    "    print(f'MCC = {metrics.matthews_corrcoef(y_ts, class_pred_ts):.3f}')\n",
    "    print(f'ACC = {metrics.accuracy_score(y_ts, class_pred_ts):.3f}')\n",
    "    print(f'SENS = {metrics.recall_score(y_ts, class_pred_ts):.3f}', '\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Hhn8hFrJLKK-"
   },
   "source": [
    "Depending on the splits I generated with train_test_split, I get the highest MCC for C=1e-5, which I take as the optimal parameter in this setting. You may get different values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "oh8VvvpcLKK_"
   },
   "source": [
    "**Optional exercise:** as you already saw in the lectures, there are many parameters that can be tuned, also when considering only one simple classifier. For example, if you consider SVM with 'rbf' kernel, you could check performance changes with different values of C **and** gamma, for example using two nested loops."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "BPtC-EBSLKK_"
   },
   "outputs": [],
   "source": [
    "## space for exercise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "hJC3MFEALKLB"
   },
   "source": [
    "As we mentioned, Scikit offers fully automated parameter tuning engine. We illustrate its power with a simple example on our data. We use GridSearchCV to search through a grid of C and gamma parameter options for SVM with 'rbf' kernel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "utM1ALBfLKLC",
    "outputId": "d96dc041-2f6f-4f1a-bca5-70310d1f79ee"
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import GridSearchCV, StratifiedShuffleSplit\n",
    "\n",
    "C_range = [0.001, 0.01, 0.1, 1, 10, 100]\n",
    "gamma_range = ['auto', 0.01, 1]\n",
    "# the parameter grid is defined as a dictionary {<parameter>: <values>} for each parameter\n",
    "param_grid = dict(gamma=gamma_range, C=C_range)\n",
    "# define the type of cross-validation for the grid search\n",
    "cv = StratifiedShuffleSplit(n_splits=10, test_size=0.5, random_state=42)\n",
    "# create a GridSearchCV object\n",
    "grid = GridSearchCV(SVC(kernel=\"rbf\"), param_grid=param_grid, cv=cv, n_jobs=4)\n",
    "# go!\n",
    "grid.fit(x_tr, y_tr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The best parameters and the corresponding average cross-validated score are available in the `best_params_` and `best_score_` attributes, respectively:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(grid.best_params_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(grid.best_score_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Note: GridSearchCV by default does not maximize the MCC, but the accuracy. It is however possible to specify a custom metric to be maximized (or minimized) using sklearn's make_scorer() function._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from sklearn.metrics import make_scorer\n",
    "# mcc_scorer = make_scorer(metrics.matthews_corrcoef)\n",
    "# grid = GridSearchCV(SVC(kernel=\"rbf\"), param_grid=param_grid, cv=cv, n_jobs=4, scorer=mcc_scorer)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Train a model with the best parameters and predict on the test set, computing a few metrics:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf = SVC(kernel=\"rbf\", C=grid.best_params_['C'], gamma=grid.best_params_['gamma'])\n",
    "clf.fit(x_tr, y_tr)\n",
    "y_pred = clf.predict(x_ts)\n",
    "print(f'MCC = {metrics.matthews_corrcoef(y_ts, y_pred):.3f}')\n",
    "print(f'ACC = {metrics.accuracy_score(y_ts, y_pred):.3f}')\n",
    "print(f'SENS = {metrics.recall_score(y_ts, y_pred):.3f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Implementing a basic Data Analysis Protocol"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a final example, we implement a 10x 5-fold Cross-validation schema with a simple feature ranking. \n",
    "\n",
    "For each CV iteration, a Random Forest model is trained on the training portion of the data, then features are ranked according to the Random Forest importances; a series of Random Forest models are built upon an increasing number of the ranked features (i.e., 1, 5, 10, etc.) and evaluated on the test data in terms of MCC.\n",
    "\n",
    "The average MCC over the 10x5 CV iterations is computed for the different feature set sizes. We choose the feature set size that maximizes the average MCC.\n",
    "\n",
    "This basic example is meant as a starting point for building more complex pipelines, i.e., with more feature steps, confidence intervals for MCC, computation of a unified ranked feature list (as in Jurman et al., _Bioinformatics_ , 2008)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "CV_N = 10 # number of CV iterations\n",
    "CV_K = 5 # number of CV folds\n",
    "FEATURE_STEPS = [1, 5, 10, 25, 50, 100]\n",
    "# prepare output MCC array\n",
    "MCC = np.empty((CV_K*CV_N, len(FEATURE_STEPS)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for n in range(CV_N):\n",
    "    print(\"~~~ Iteration %d ~~~\" % (n+1))\n",
    "    skf = StratifiedKFold(n_splits=CV_K, shuffle=True, random_state=n)\n",
    "    for i, (idx_tr, idx_ts) in enumerate(skf.split(x_tr, y_tr)):\n",
    "        print(\"Fold %d\" % (i+1))\n",
    "        X_train, Y_train = x_tr[idx_tr], y_tr[idx_tr]\n",
    "        X_test, Y_test = x_tr[idx_ts], y_tr[idx_ts]\n",
    "        \n",
    "        clf = RFC(n_estimators=500, random_state=n)\n",
    "        clf.fit(X_train, Y_train)\n",
    "        ranking = np.argsort( clf.feature_importances_ )[::-1]\n",
    "        \n",
    "        for j, s in enumerate(FEATURE_STEPS):\n",
    "            v = ranking[:s] # consider the top s ranked features\n",
    "            X_tr_fs, X_ts_fs = X_train[:, v], X_test[:, v] # extract them from internal train and test data\n",
    "            clf.fit(X_tr_fs, Y_train) # train a classifier on the reduced train dataset\n",
    "            yp = clf.predict(X_ts_fs) # predict on the reduced test dataset\n",
    "            MCC[(n*CV_K)+i, j] = metrics.matthews_corrcoef(Y_test, yp) # evaluate the model performance\n",
    "        \n",
    "    print()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# np.save(\"MCC_CV\", MCC)\n",
    "# MCC = np.load(\"MCC_CV.npy\")\n",
    "MCC.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "MCC_avg = np.mean(MCC, axis=0)\n",
    "MCC_max = np.max(MCC_avg)\n",
    "n_feats = FEATURE_STEPS[np.argmax(MCC_avg)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# average MCC for each feature step\n",
    "for nf, mcc in zip(FEATURE_STEPS, MCC_avg):\n",
    "    print(\"nf = %d, MCC = %.2f\" % (nf, mcc))\n",
    "\n",
    "print()\n",
    "\n",
    "print(\"Best MCC = %.2f with %d features\" % (MCC_max, n_feats))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.title(\"Average MCC\")\n",
    "plt.plot(FEATURE_STEPS, MCC_avg, 'o-')\n",
    "plt.ylim((0.0, 1.0))\n",
    "plt.xlabel(\"Feature steps\")\n",
    "plt.ylabel(\"MCC\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "name": "practical_neuroblastoma_partII_v0.3.ipynb",
   "provenance": [],
   "version": "0.3.2"
  },
  "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}