1.1_Explainability.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualisations and Explainabilty of Machine Learning models for time series data\n",
    "\n",
    "In the 0_Intro.ipynb and 1_Baseline.ipynb notebooks, we create a machine learning pipeline for the classification of accelerometry data. For the pipeline developed, we have found the classification of activity labels in the test set is not perfect, hence we would like to make improvements to the model. This is a difficult task to do, and requires visualisations and model explainability to be done efficiently and effectively.\n",
    "\n",
    "## AWS info\n",
    "Run this notebook with the **Data Science** image and **ml.t3.2xlarge** instance type.\n",
    "\n",
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from imblearn.ensemble import BalancedRandomForestClassifier\n",
    "from sklearn import metrics\n",
    "from sklearn import preprocessing\n",
    "from sklearn import manifold\n",
    "import seaborn as sns\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "%matplotlib inline\n",
    "\n",
    "# For reproducibility\n",
    "np.random.seed(42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Path to your extracted windows\n",
    "DATASET_PATH = 'processed_data/'\n",
    "X_FEATS_PATH = 'X_feats.pkl'  # path to your extracted features, if have one\n",
    "print(f'Content of {DATASET_PATH}')\n",
    "print(os.listdir(DATASET_PATH))\n",
    "\n",
    "X = np.load(DATASET_PATH+'X.npy', mmap_mode='r')\n",
    "Y = np.load(DATASET_PATH+'Y.npy')\n",
    "T = np.load(DATASET_PATH+'T.npy')\n",
    "pid = np.load(DATASET_PATH+'pid.npy')\n",
    "X_feats = pd.read_pickle(DATASET_PATH+'X_feats.pkl')\n",
    "\n",
    "# As before, let's map the text annotations to simplified labels\n",
    "ANNO_LABEL_DICT_PATH = 'capture24/annotation-label-dictionary.csv'\n",
    "anno_label_dict = pd.read_csv(ANNO_LABEL_DICT_PATH, index_col='annotation', dtype='string')\n",
    "Y = anno_label_dict.loc[Y, 'label:Walmsley2020'].to_numpy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first task in understanding any dataset is to visualise the data. Here we plot 5 example plots for each label."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_acc_samples(X, Y, n_samples=5, plot_norm=False):\n",
    "    labels = np.unique(Y)\n",
    "    num_labels = len(labels)\n",
    "    \n",
    "    fig, axs = plt.subplots(num_labels, n_samples, figsize=(16,8))\n",
    "    fig.set_facecolor('white')\n",
    "    for i in range(num_labels):\n",
    "        axs[i,0].set_ylabel(labels[i])\n",
    "        idxs = np.where(Y==labels[i])[0]\n",
    "        for j in range(n_samples):\n",
    "            if plot_norm:\n",
    "                axs[i,j].plot(np.linalg.norm(X[idxs[j]], axis=1)-1)\n",
    "                axs[i,j].set_ylim([-2,2])\n",
    "                legend = ['norm']\n",
    "            else:\n",
    "                axs[i,j].plot(X[idxs[j]])\n",
    "                legend = ['x', 'y', 'z']\n",
    "            axs[i,j].set_xticks([])\n",
    "            axs[i,j].set_yticks([])\n",
    "    fig.legend(legend)\n",
    "    plot_desc = \"normalised\" if plot_norm else \"raw\"\n",
    "    plt.suptitle(\"Example {} plots for each label in Capture 24 dataset\".format(plot_desc))\n",
    "    fig.tight_layout()\n",
    "\n",
    "plot_acc_samples(X, Y, plot_norm=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our raw data contained 30 second segments of acceleration in x, y and z. The acceleration in these 3 axes can overlap with each other when plotting, limiting the visualisation of the signal. For this reason, it is sometimes more useful to plot the euclidean norm of acceleration, and subtracting 1g to remove gravity.\n",
    "\n",
    "Note: All plots of accleration norm minus one are limited between -2g and 2g."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_acc_samples(X, Y, plot_norm=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 1**: Inspecting the differences between the raw and normalised signal, what are the disadvantages to plotting only the normalisaed signal to observe patterns."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 2**: What are the noticeable differences between the examples for each label? Is this as you expected?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 3**: For simplicity of code, the plot_acc_samples function plots the first n accelerometry samples of a given label. Why might there be an issue with using these samples to represent the entire dataset? What one line of code could be added to address this issue?\n",
    "\n",
    "Hint: How are the X and Y arrays ordered? What is the meaning of the first n samples?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## t-SNE visualization\n",
    "\n",
    "Instead of plotting the raw or minimally processed accelerometry data, we can instead using unsupervised dimension reduction methods to view the data. While not perfect, if clusters of a certain label are clearly grouped in a t-SNE plot, it is expected that our model classifier should have good performance in predicting this label.\n",
    "\n",
    "Note: The data needs to be downsampled 100x, to run the t-SNE in a relevant time window."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def scatter_plot(X, Y):\n",
    "    unqY = np.unique(Y)\n",
    "    fig, ax = plt.subplots()\n",
    "    plt.title(\"t-SNE of extracted features of Capture24 data\")\n",
    "    for y in unqY:\n",
    "        X_y = X[Y==y]\n",
    "        ax.scatter(X_y[:,0], X_y[:,1], label=y, alpha=.5, s=10)\n",
    "    fig.legend()\n",
    "    \n",
    "print(\"Plotting t-SNE on extracted features...\")\n",
    "\n",
    "scaler = preprocessing.StandardScaler()  # PCA requires normalized data\n",
    "\n",
    "X_scaled = scaler.fit_transform(X[::100].reshape(X[::100].shape[0],-1))\n",
    "\n",
    "tsne = manifold.TSNE(n_components=2,  # project down to 2 components\n",
    "    init='random', random_state=42, perplexity=100, learning_rate='auto')\n",
    "X_tsne_pca = tsne.fit_transform(X_scaled)\n",
    "scatter_plot(X_tsne_pca, Y[::100])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Train/Test split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Hold out participants P101-P151 for testing (51 participants)\n",
    "test_ids = [f'P{i}' for i in range(101,152)]\n",
    "mask_test = np.isin(pid, test_ids)\n",
    "mask_train = ~mask_test\n",
    "X_train, Y_train, T_train, pid_train = \\\n",
    "    X_feats[mask_train], Y[mask_train], T[mask_train], pid[mask_train]\n",
    "X_test, Y_test, T_test, pid_test, X_raw_test = \\\n",
    "    X_feats[mask_test], Y[mask_test], T[mask_test], pid[mask_test], X[mask_test]\n",
    "print(\"Shape of X_train:\", X_train.shape)\n",
    "print(\"Shape of X_test:\", X_test.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train a random forest classifier\n",
    "\n",
    "*This may take a while*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Argument oob_score=True to be used for HMM smoothing (see later below)\n",
    "clf = BalancedRandomForestClassifier(\n",
    "    n_estimators=1000,\n",
    "    replacement=True,\n",
    "    sampling_strategy='not minority',\n",
    "    oob_score=True,\n",
    "    n_jobs=-1,\n",
    "    random_state=42,\n",
    "    verbose=1\n",
    ")\n",
    "clf.fit(X_train, Y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model performance\n",
    "\n",
    "The classification report gives standard metrics for measure the performance of a model, including accuracy and macro F1 score. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y_test_pred = clf.predict(X_test)\n",
    "print('\\nClassifier performance')\n",
    "print('Out of sample:\\n', metrics.classification_report(Y_test, Y_test_pred))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another way to visualise the model performance is the confusion matrix, comparing the true test labels to the predicted test labels. This gives us an understanding how the model incorrectly assigns labels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,8))\n",
    "conf = pd.DataFrame(metrics.confusion_matrix(Y_test, Y_test_pred, normalize='true'), index=clf.classes_, columns=clf.classes_)\n",
    "sns.set(font_scale = 2)\n",
    "sns.heatmap(conf, annot=True, fmt='.2f', ).set(title=\"Confusion matrix\", xlabel=\"Predicted\", ylabel=\"True\")\n",
    "plt.xticks(rotation=45)\n",
    "sns.set(font_scale = 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When observing the performance of a model, it is important to inspect the failure cases of the model. That is, the samples in which the model prediction does not match up to the true label. Inspecting this offers the possibility of identifying some potential causes of failure, such as:\n",
    "- Lack of model generalisability\n",
    "- Poor/inaccurate labelling of truth labels\n",
    "- Non trivial differences between labels "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_failure_cases(X, Y, Y_pred, label, n_samples=5):\n",
    "    labels = np.unique(Y_pred)\n",
    "    num_labels = len(labels)\n",
    "    \n",
    "    fig, axs = plt.subplots(num_labels, n_samples, figsize=(16,8))\n",
    "    fig.supylabel('Predicted label')\n",
    "    fig.set_facecolor('white')\n",
    "    for i in range(num_labels):\n",
    "        idxs = np.where((Y_pred==labels[i]) & (Y==label))[0]\n",
    "        axs[i,0].set_ylabel(\"{}: {:.1f}%\".format(labels[i], 100*len(idxs)/sum(Y==label)))\n",
    "        for j in range(n_samples):\n",
    "            axs[i,j].plot(np.linalg.norm(X[idxs[j]], axis=1)-1)\n",
    "            axs[i,j].set_ylim([-2,2])\n",
    "            axs[i,j].set_xticks([])\n",
    "            axs[i,j].set_yticks([])\n",
    "    fig.legend(['norm'])\n",
    "    plt.suptitle(\"Example normalised plots for labelled {} activity\".format(label))\n",
    "    fig.tight_layout()\n",
    "\n",
    "plot_failure_cases(X_raw_test, Y_test, Y_test_pred, label='sleep')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 4**: Plot the failure cases for the other labelled activities. What do you think is the main cause for these failure cases?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Feature selection\n",
    "\n",
    "The selection of features to be extracted for this classifaction task is a key factor in improving model performance. In previous notebooks, exercises were given to encourage you to expand the list of features extracted from the accelerometry windows. We note however that extracting more features requires more compute power, and takes a longer time to run. In some cases, these extract features may add minimal improvements to performance. We hence seek to find efficient ways to determine which of the extracted features can be removed/replaced while maintaining good model performance.\n",
    "\n",
    "### Feature correlation\n",
    "\n",
    "Feature correlation gives an indication of how a pair of features are associated with one another. A pair of features with a high correlation coefficient (close to 1) tends to make an inefficient model, as we only need one of these features to extract the information required for classification. A visualisation of the correlation of all extracted features can be seen in a correlation matrix, and any pairs of features with high values can be removed from the feature extraction pipeline."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,10))\n",
    "sns.heatmap(X_feats.corr().abs()).set(title=\"Correlation matrix of extracted features\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 5**: What is the correlation coefficient for the mean acceleration in the x and y axis (xmean and ymean)? Do you think both of these features need to be extracted for an efficient model? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Feature importance\n",
    "\n",
    "Feature importance is an explainable AI technique to reveal the relative significance of individual features on model outputs. There are many different methods that can be used to determine feature importance, however in this notebook, we will use GINI importance. When looking to optimise feature extraction (less compute power and time, better model performance), features with lowest importance should be removed first.\n",
    "\n",
    "#### GINI Importance\n",
    "GINI importance is a feature importance method that can be extracted directly from the BalancedRandomForestClassifier class. More information on how exactly it works can be found [here](https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-10-213).\n",
    "Note: GINI importance is known to perform poorly with highly correlated features. If relying only on GINI importance to determine which features to remove, highly correlated features should be removed first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(4,10))\n",
    "sns.barplot(x=clf.feature_importances_, y=X_feats.columns).set(title=\"GINI Feature importance\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise 5**: Which features have highest GINI importance? Does this match your expectation?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# References\n",
    "\n",
    "GINI importance\n",
    "\n",
    "- [A comparison of random forest and its Gini importance with standard chemometric methods for the feature selection and classification of spectral data](https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-10-213)\n",
    "- [Feature Importance Measures for Tree Models — Part I](https://medium.com/the-artificial-impostor/feature-importance-measures-for-tree-models-part-i-47f187c1a2c3)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.12 ('wearables_workshop')",
   "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.9.12"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "9c3f0ec06e35a01db1886840759b15ced8e4cb884230ad4ac6579f6dd8670105"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}