chierici_practical_part2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "TJfNUgelLKKI"
   },
   "source": [
    "# <center>Machine learning from scratch - Part II</center>\n",
    "## <center>WebValley 2019 - AI for Health @ Casez, Italy</center>\n",
    "### <center>Marco Chierici & Margherita Francescatto</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 pandas as pd ## for reading text files and manipulating data frames\n",
    "from sklearn import neighbors ## kNN classifier\n",
    "from sklearn import svm ## SVM classifier\n",
    "from sklearn.ensemble import RandomForestClassifier ## RF classifier\n",
    "from sklearn.model_selection import cross_val_score ## needed to train in CV\n",
    "from bokeh.io import output_notebook, show\n",
    "from bokeh.plotting import figure\n",
    "np.random.seed(42) ## set random seed just in case"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "qOlPPJmmLKKO"
   },
   "source": [
    "Define files to read:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##  for convenience, define the data directory as a variable\n",
    "DATA_DIR = \"NB_data/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "cnidnnC4LKKP"
   },
   "outputs": [],
   "source": [
    "DATA_TR = DATA_DIR + \"MAV-G_272_tr.txt\"\n",
    "DATA_TS = DATA_DIR + \"MAV-G_272_ts.txt\"\n",
    "LABS_TR = DATA_DIR + \"labels_tr.txt\"\n",
    "LABS_TS = DATA_DIR + \"labels_ts.txt\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "yd4YsVYnLKKS"
   },
   "source": [
    "Read in the files as _pandas dataframes_ (they are conceptually like R data frames):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "SBBKnUO-LKKS"
   },
   "outputs": [],
   "source": [
    "data_tr = pd.read_csv(DATA_TR, sep = \"\\t\")\n",
    "data_ts = pd.read_csv(DATA_TS, sep = \"\\t\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "2QiDbifgLKKV"
   },
   "source": [
    "Since we already looked at the data in the first part of the dataset, we move directly to the juicy stuff."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "nY4Jcd9XLKKW"
   },
   "source": [
    "We drop the first column from the train and test expression sets, since it's just the sample IDs..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "-Ejm4N6xLKKW"
   },
   "outputs": [],
   "source": [
    "data_tr = data_tr.drop('sampleID',axis =1)\n",
    "data_ts = data_ts.drop('sampleID',axis =1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Q6DrKsxlLKKZ"
   },
   "source": [
    "...and store the data into Numpy arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "Mua2Ajr-LKKa"
   },
   "outputs": [],
   "source": [
    "x_tr = data_tr.values\n",
    "x_ts = data_ts.values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "6yvTIil4LKKc"
   },
   "source": [
    "Now we read in the files containing labels and select the column with the CLASS target to do our first round of analyses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "69bKNl0QLKKe"
   },
   "outputs": [],
   "source": [
    "labs_tr = pd.read_csv(LABS_TR, sep = \"\\t\")\n",
    "labs_ts = pd.read_csv(LABS_TS, sep = \"\\t\")\n",
    "class_lab_tr = labs_tr[['CLASS']]\n",
    "class_lab_ts = labs_ts[['CLASS']]\n",
    "y_tr = class_lab_tr.values.ravel()\n",
    "y_ts = class_lab_ts.values.ravel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "arErpTdvLKKh"
   },
   "source": [
    "In the previous part of the tutorial, we focused on the k-NN classifiers. In the previous lecture, however, we explored theoretical aspects related to two other broadly used classifiers: Support Vector Machines (SVMs) and Random Forests (RFs). In this second 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": {
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import MinMaxScaler\n",
    "## first you need to create a \"scaler\" object\n",
    "scaler = 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": "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(random_state=42, verbose=1, 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": {
    "colab_type": "text",
    "id": "waXcgaLpLKKm"
   },
   "source": [
    "As previously done with the k-NN classifier, we fit the SVM and get the predictions for the test data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "Qqc3TmFBLKKn",
    "outputId": "d9ef6c64-9f18-4bea-9167-decaa0ca1820"
   },
   "outputs": [],
   "source": [
    "## fit the model and get the predictions\n",
    "svc.fit(x_tr, y_tr)\n",
    "class_pred_ts = svc.predict(x_ts)"
   ]
  },
  {
   "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. to access the functions, we need to load the metrics module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "Ku0JSF_ALKKs",
    "outputId": "94585c0e-534a-445d-d0ba-92a9bf3a9388"
   },
   "outputs": [],
   "source": [
    "from sklearn import metrics\n",
    "print('MCC = ', metrics.matthews_corrcoef(class_lab_ts, class_pred_ts))\n",
    "print('ACC = ', metrics.accuracy_score(class_lab_ts, class_pred_ts))\n",
    "print('SENS = ', metrics.recall_score(class_lab_ts, class_pred_ts))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "1JnMF6UoLKKw"
   },
   "source": [
    "We can also give a look at the classification report."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "whSZnHGALKKx",
    "outputId": "2c471734-3504-4af7-8ebb-74e5a02be301"
   },
   "outputs": [],
   "source": [
    "print(metrics.classification_report(class_lab_ts, class_pred_ts))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "K6ZZSMPkLKKz"
   },
   "source": [
    "Exercise: **one-shot Random Forest classification**. _Hint:_ the RF classifier is implemented in the Scikit learn class RandomForestClassifier, from _sklearn.ensemble_ module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "ZT6XjB20LKK0"
   },
   "outputs": [],
   "source": [
    "## space for exercise\n"
   ]
  },
  {
   "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": {
    "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('C = ', C)\n",
    "    svc = svm.SVC(kernel = 'linear', C=C)\n",
    "    svc.fit(x_tr, class_lab_tr.values.ravel())\n",
    "    class_pred_ts = svc.predict(x_ts)\n",
    "    print('MCC = ', metrics.matthews_corrcoef(class_lab_ts, class_pred_ts))\n",
    "    print('ACC = ', metrics.accuracy_score(class_lab_ts, class_pred_ts))\n",
    "    print('SENS = ', metrics.recall_score(class_lab_ts, class_pred_ts), \"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Hhn8hFrJLKK-"
   },
   "source": [
    "From C = 1e-03 the classification performance reaches a plateau. C = 1e-04 yields the highest MCC on the test set: when tuning the parameters we would consider this as the best choice for the problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "oh8VvvpcLKK_"
   },
   "source": [
    "**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. In order to do this we define a small function svc_param_selection that does the work for us."
   ]
  },
  {
   "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\n",
    "\n",
    "def svc_param_selection(X, y, nfolds):\n",
    "    Cs = [0.001, 0.01, 0.1, 1, 10]\n",
    "    gammas = [0.001, 0.01, 0.1, 1]\n",
    "    param_grid = {'C': Cs, 'gamma' : gammas}\n",
    "    grid_search = GridSearchCV(svm.SVC(kernel='rbf'), param_grid, cv=nfolds, n_jobs=4)\n",
    "    grid_search.fit(X, y)\n",
    "    grid_search.best_params_\n",
    "    return grid_search.best_params_\n",
    "\n",
    "svc_param_selection(x_tr, y_tr, 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "tI_6v1jNLKLF"
   },
   "source": [
    "## Feature ranking"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "HcQ4gJjgLKLG"
   },
   "source": [
    "As mentioned in the lecture, random forests have a built-in tool for feature ranking"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "2lZAaTXJLKLH",
    "outputId": "2155231c-e50c-4c06-82c4-6b6a5f7c4ee2"
   },
   "outputs": [],
   "source": [
    "# Build a forest and compute the feature importances\n",
    "rf = RandomForestClassifier(n_estimators=250)\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": [
    "class_pred_ts = rf.predict(x_ts)\n",
    "print('MCC = ', metrics.matthews_corrcoef(class_lab_ts, class_pred_ts))\n",
    "print('ACC = ', metrics.accuracy_score(class_lab_ts, class_pred_ts))\n",
    "print('SENS = ', metrics.recall_score(class_lab_ts, class_pred_ts))"
   ]
  },
  {
   "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": [
    "importances = rf.feature_importances_\n",
    "indices = np.argsort(importances)[::-1]\n",
    "\n",
    "# Print the feature ranking\n",
    "print(\"Feature ranking (top 10 features):\")\n",
    "for f in range(10):\n",
    "    print(\"%d. feature %d (%f)\" % (f + 1, indices[f], importances[indices[f]]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "VX9otFjOLKLX"
   },
   "source": [
    "Would be nice to know to which genes they actually correspond. If you remember the gene names are the column names of the pandas dataframe containing the training/test data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "2fSkitN7LKLY",
    "outputId": "73191a71-9657-4582-ede6-7fb14cd3fc05"
   },
   "outputs": [],
   "source": [
    "columnsNamesArr = data_tr.columns.values\n",
    "for i in range(10):\n",
    "    print(columnsNamesArr[indices[i]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "kkgGzniAP6y4"
   },
   "source": [
    "## Extra exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "uSW_fBLwP-fO"
   },
   "source": [
    "The classification task considered so far (CLASS) is quite easy, since the classes reflect extreme disease outcomes (favorable vs unfavorable).\n",
    "\n",
    "A more interesting task could be the prediction of Event-Free Survival (EFS). To do this, an extended version of the dataset is provided in the `/data/marco` directory:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "QfYqaV0ZP9Y-"
   },
   "outputs": [],
   "source": [
    "DATA_TR = DATA_DIR + \"full_MAV-G_498_tr.csv\"\n",
    "DATA_TS = DATA_DIR + \"full_MAV-G_498_ts.csv\"\n",
    "LABS_TR = DATA_DIR + \"full_labels_tr.txt\"\n",
    "LABS_TS = DATA_DIR + \"full_labels_ts.txt\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "L7i0MA0FSa7Z"
   },
   "source": [
    "Read the data in (take note of the input data format) and prepare the `x_tr`, `x_ts`,  `y_tr`, `y_ts` Numpy arrays, as before, using \"EFS\" as target variable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "UNpBL3UDS8dI"
   },
   "source": [
    "Recalling concepts from the first practical, perform an explorative PCA analyisis, plotting the first two components."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "zNfPACYqTOT4"
   },
   "source": [
    "Train a kNN classifier in one-shot mode: fit the model on the training set and predict the labels on the test set. Compute performance metrics using the provided true labels of the test set."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "6OkabT2FTvPT"
   },
   "source": [
    "Experiment with different classifier(s) and/or different parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "tHhMWGSLUXS5"
   },
   "source": [
    "Try tuning the parameter(s) (e.g. using GridSearchCV) and find the optimal parameter set."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "_vk_ChapUzbP"
   },
   "source": [
    "Using the optimal parameters, run a 10x iterated 5-fold cross-validation on the training set; compute the average cross-validation metrics."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "NjyflIl5VF2Y"
   },
   "source": [
    "Using the optimal parameters, train a model on the whole training set and predict the labels of the test set. Compute the metrics and compare them with the average cross-validation metrics. What do you expect? Use the trained model to rank the features and inspect the top ones."
   ]
  }
 ],
 "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}