BIDS_06_Clustering_2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "591c8891",
   "metadata": {},
   "source": [
    "Aim to complete parts of this tutorial on your own *before* the practical session (after the lecture).\n",
    "\n",
    "Use the practical session to get help for any aspect you do not understand or were unable to complete."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "059df924",
   "metadata": {},
   "source": [
    "# Clustering 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79e03dc1",
   "metadata": {},
   "source": [
    "Learning objectives\n",
    "1. Apply GMM and DBSCAN (density-based clustering) to different data sets and interpret the outputs using the popular python library [sklearn](https://scikit-learn.org/stable/)\n",
    "2. Interpret the results to learn about the data structure\n",
    "3. Explore differences based on changing clustering parameters\n",
    "4. Visualise the clustering results\n",
    "\n",
    "Optional learning objective (view the video online about OPTICS before you start)\n",
    "\n",
    "5. Apply OPTICS (extension of DBSCAN) to different data sets, interpret the results, explore differences and visualise the results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e126d434",
   "metadata": {},
   "source": [
    "## Import specific packages and functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "959c3008",
   "metadata": {},
   "outputs": [],
   "source": [
    "# IMPORTS\n",
    "import pandas as pd \n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.decomposition import PCA, SparsePCA, KernelPCA\n",
    "from sklearn.preprocessing import StandardScaler, PowerTransformer, RobustScaler\n",
    "from sklearn.cluster import DBSCAN, OPTICS\n",
    "from sklearn.neighbors import NearestNeighbors\n",
    "from sklearn.mixture import GaussianMixture\n",
    "from sklearn.manifold import TSNE, MDS\n",
    "from sklearn import metrics\n",
    "import umap\n",
    "import warnings\n",
    "warnings.simplefilter('ignore', category=FutureWarning)\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ccf56a3",
   "metadata": {},
   "source": [
    "## Load in dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "140495b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in the dataset\n",
    "covid19_proteomics = pd.read_excel('../Data-main/COVID19_proteomics.xlsx')\n",
    "\n",
    "# Scale the data with standard scaling (0 mean and unit variance)\n",
    "covid19_proteomics_scaled = StandardScaler().fit_transform(covid19_proteomics.iloc[:, 3:])\n",
    "\n",
    "# Create a dataframe with adenoma/carcinoma status and the scaled data\n",
    "scaled_df = pd.DataFrame(covid19_proteomics_scaled)\n",
    "scaled_df['COVID19'] = covid19_proteomics.COVID19"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c0f3c28",
   "metadata": {},
   "source": [
    "## Gaussian Mixture Models (GMM)\n",
    "\n",
    "A Gaussian mixture model is a probabilistic model that assumes all the data points are generated from a mixture of a finite number of Gaussian distributions with unknown parameters. One can think of mixture models as generalizing k-means clustering to incorporate information about the covariance structure of the data as well as the centers of the latent Gaussians. Read more [here](https://scikit-learn.org/stable/modules/mixture.html).\n",
    "\n",
    "We will use the [`sklearn.mixture.GaussianMixture()`](https://scikit-learn.org/stable/modules/generated/sklearn.mixture.GaussianMixture.html#sklearn.mixture.GaussianMixture) object to fit the model. Similar to dimensionality reduction and other clustering methods, one of the parameters to set for GMM models is the `n_components` (number of components). Another parameter is the `covariance_type` which can be one of the following options: \n",
    "- ‘full’: each component has its own general covariance matrix\n",
    "- ‘tied’: all components share the same general covariance matrix\n",
    "- ‘diag’: each component has its own diagonal covariance matrix\n",
    "- ‘spherical’: each component has its own single variance\n",
    "\n",
    "Here we will fit a GMM model and explore through visualisations how "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bdd69da9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# enter your CID here, or date of birth, or another number of your choosing to use as random state\n",
    "CID = 0\n",
    "\n",
    "# remember to check the documentation of each algorithm if setting the random_state is needed\n",
    "# for this tutorial and all upcoming tutorials..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67ee5d82-2ace-40db-a011-7aed66607938",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Split the data into train and test sets\n",
    "X_train, X_test, y_train, y_test = train_test_split(scaled_df.iloc[:,:-1], scaled_df.COVID19, test_size=0.2, random_state=CID)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eba6e837-3cd2-40fb-8c32-929f1756fec6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# fit a Gaussian Mixture Model with two components\n",
    "gmm = GaussianMixture(n_components=2, covariance_type=\"full\")\n",
    "y_pred = gmm.fit_predict(X_train)\n",
    "\n",
    "# here we take the first 2 proteomic features for the x and y axis \n",
    "sns.scatterplot(X_train.iloc[:, 0], X_train.iloc[:, 1], c=y_pred, s=30, cmap='vlag')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b100e090-ced2-4daf-8bab-4146f83c20b7",
   "metadata": {},
   "source": [
    "We can use PCA or another dimensionality reduction method (see previous tutorials on dimension reduction) to create components that encopass information from all our variables, and not just 2 random variables we select for visualisation.\n",
    "\n",
    "After we visualise the PCA results with the current GMM model, we will explore the different covariance options and how (or if) they change our clusters.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8871d996-3d62-492a-bc35-46fb5ad262a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "pca = PCA(n_components=2)\n",
    "pca_df = pca.fit_transform(X_train)\n",
    "sns.scatterplot(pca_df[:, 0], pca_df[:, 1], c=y_pred, s=30, cmap='vlag')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "620b238b-a40c-4a78-8b75-ceea2cad6662",
   "metadata": {},
   "outputs": [],
   "source": [
    "# function for fitting GMM models with the 4 different covariance options\n",
    "def test_gmm_covariance (cov):\n",
    "    gmm = GaussianMixture(n_components=2, covariance_type=cov)\n",
    "    y_pred = gmm.fit_predict(X_train)\n",
    "    #pca = PCA(n_components=2)\n",
    "    #pca_df = pca.fit_transform(X_train)\n",
    "    return y_pred\n",
    "\n",
    "cov_list = ['full', 'tied', 'diag', 'spherical']\n",
    "\n",
    "# creating subplots \n",
    "plt.figure(figsize=(10,7))\n",
    "for n, cov_type in enumerate(cov_list):\n",
    "    ax = plt.subplot(2, 2, n+1)\n",
    "    sns.scatterplot(pca_df[:, 0], pca_df[:, 1], \n",
    "                    c=test_gmm_covariance(cov_type), s=30, cmap='vlag', ax=ax)\n",
    "    ax.set_title(cov_type.upper())\n",
    "    if n>1:\n",
    "        ax.set_xlabel(\"PC1\")\n",
    "    if (n % 2) == 0:\n",
    "        ax.set_ylabel(\"PC2\")\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4db4b3fc",
   "metadata": {},
   "source": [
    "### Compare the result from GMM with k-Means with 2 clusters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c40cc2c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code goes here, reuse materials from BIDS 5\n",
    "\n",
    "# Initialise a k-Means model\n",
    "\n",
    "# Fit a k-Means model to X_train\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c81ba4a",
   "metadata": {},
   "source": [
    "Once you have your data ready, prompt an code-LLM to write the visualisation code for you."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "514f561f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the results side to side with GMM as comparison (use subplot similar to the above example)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23b6edd5",
   "metadata": {},
   "source": [
    "How are these two methods (GMM and k-Means) similar? When do you observe differences?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61e92025",
   "metadata": {},
   "source": [
    "## Density-Based Spatial Clustering of Applications with Noise (DBSCAN)\n",
    "\n",
    "DBSCAN is a density-based clustering non-parametric algorithm. Given a set of points in some space, it groups together points that are closely packed together (points with many nearby neighbours - this should ring a bell about a method from BIDS 5), marking as outliers points that lie alone in low-density regions (whose nearest neighbours are too far away). Original paper found [here](https://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.121.9220).\n",
    "\n",
    "We will use the [`sklearn.cluster.DBSCAN()`](https://scikit-learn.org/stable/modules/generated/sklearn.cluster.DBSCAN.html#sklearn.cluster.DBSCAN) object to fit our model. Below are some of the parameters that will need to be set, you will notice they are similar to other clustering techniques (i.e. k-means): \n",
    "\n",
    "- `eps`: The maximum distance between two samples for one to be considered as in the neighbourhood of the other. This is not a maximum bound on the distances of points within a cluster. This is the most important DBSCAN parameter to choose appropriately for your data set and distance function.\n",
    "- `min_samples`: The number of samples (or total weight) in a neighbourhood for a point to be considered as a core point. This includes the point itself.\n",
    "- `metric`: The metric to use when calculating distance between instances in a feature array.\n",
    "\n",
    "Before fitting our DBSCAN model we will examine the best values for epsilon by using the `NearestNeighbors` function. After plotting the distances we get, we will find the optimal value for epsilon at the point of maximum curvature:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5dd60d17-d8b3-458e-be0e-d1ddac84dac8",
   "metadata": {},
   "outputs": [],
   "source": [
    "NN = NearestNeighbors(n_neighbors=2)\n",
    "nbrs = NN.fit(scaled_df.iloc[:, :-1])\n",
    "distances, indices = nbrs.kneighbors(scaled_df.iloc[:, :-1])\n",
    "\n",
    "distances = np.sort(distances, axis=0)\n",
    "distances = distances[:,1]\n",
    "# plt.plot(distances)\n",
    "sns.lineplot(data=distances)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae9d5322",
   "metadata": {},
   "outputs": [],
   "source": [
    "dbscan = DBSCAN(\n",
    "    eps=19, \n",
    "    min_samples=3, \n",
    "    metric='euclidean'\n",
    ").fit(scaled_df.iloc[:, :-1])\n",
    "\n",
    "labels = dbscan.labels_\n",
    "n_clusters_ = len(set(labels))\n",
    "n_noise_ = list(labels).count(-1)\n",
    "\n",
    "print(\"Estimated number of clusters: %d\" % n_clusters_)\n",
    "print(\"Estimated number of noise points: %d\" % n_noise_)\n",
    "print(\"Homogeneity: %0.3f\" % metrics.homogeneity_score(scaled_df.iloc[:, -1], labels))\n",
    "print(\"Completeness: %0.3f\" % metrics.completeness_score(scaled_df.iloc[:, -1], labels))\n",
    "print(\"Silhouette Coefficient: %0.3f\" % metrics.silhouette_score(scaled_df.iloc[:, :-1], labels))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1428329a-62ff-469f-9985-bec1e4dae24d",
   "metadata": {},
   "outputs": [],
   "source": [
    "tsne_df = TSNE(\n",
    "    n_components=2, \n",
    "    perplexity=10,\n",
    ").fit_transform(scaled_df.iloc[:, :-1]) # don't forget you should optimise the t-SNE parameters as well (see BIDS 4)\n",
    "\n",
    "ax = sns.scatterplot(tsne_df[:, 0], tsne_df[:, 1], c=labels, s=30, cmap='rainbow')\n",
    "ax.set_xlabel(\"tSNE 1\")\n",
    "ax.set_ylabel(\"tSNE 2\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb09d987-3226-488d-a594-b84c3dcf6f2d",
   "metadata": {},
   "source": [
    "As you can see the clusters found by DBSCAN don't really correspond to the t-SNE visualisations. Why do you think this is? This about the data that is used as input for both methods.\n",
    "\n",
    "Now, try reducing the dimensions of the dataset with PCA and see if your results are more informative if you use the PCA scores (or another dimension reduction method) as input to DBSCAN. Alternatively, use the coordinates output from t-SNE as input to DBSCAN as well. Make sure to also tweak the parameters of the DBSCAN model to check how sensitive it is to even small changes! You can also rerun the GMM models with t-SNE or kernel t-SNE decompositions and look at how the visualisations change.\n",
    "\n",
    "Feel free to utilise and modify the code in the cell below to do a side-to-side comparison of the parameter you are changing. The example below still uses t-SNE as the dimensionality reduction method, make sure whichever method you pick that you set the same `random_state` parameter (to your CID), so that the embeddings are the same across all iterations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80d5fd49-9011-4aed-9bc6-62ce4eea9d7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_dbscan(eps):\n",
    "    dbscan = DBSCAN(\n",
    "    eps=eps, \n",
    "    min_samples=3, \n",
    "    metric='euclidean', \n",
    ").fit(scaled_df.iloc[:, :-1]) # uses the original data\n",
    "    \n",
    "    labels = dbscan.labels_\n",
    "    return labels\n",
    "\n",
    "tsne_df = TSNE( # bear in mind you can exchange this for UMAP as well\n",
    "    n_components=2, \n",
    "    perplexity=10,\n",
    "    random_state=CID,\n",
    ").fit_transform(scaled_df.iloc[:, :-1]) # uses the original data\n",
    "\n",
    "eps_values = np.arange(5, 30, 3)\n",
    "\n",
    "# creating subplots \n",
    "plt.figure(figsize=(10,7))\n",
    "for n, eps in enumerate(eps_values):\n",
    "    ax = plt.subplot(3, 3, n+1)\n",
    "    tdb = test_dbscan(eps)\n",
    "    sns.scatterplot(tsne_df[:, 0], tsne_df[:, 1], \n",
    "                    c=tdb, s=30, cmap='rainbow', ax=ax)\n",
    "    ax.set_title('eps='+str(eps)+' n_clusters='+str(len(set(tdb))))\n",
    "    if (n % 3) == 0:\n",
    "        ax.set_ylabel(\"t-SNE 2\")\n",
    "    if n>5:\n",
    "        ax.set_xlabel(\"t-SNE 1\")\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "383abfae",
   "metadata": {},
   "outputs": [],
   "source": [
    "pca = PCA(n_components=2)\n",
    "pca_df = pca.fit_transform(X_train)\n",
    "\n",
    "def test_dbscan2(eps):\n",
    "    dbscan = DBSCAN(\n",
    "    eps=eps, \n",
    "    min_samples=3, \n",
    "    metric='euclidean', \n",
    "    ).fit(pca_df)\n",
    "    \n",
    "    labels = dbscan.labels_\n",
    "    return labels\n",
    "\n",
    "eps_values = np.arange(5, 30, 3)\n",
    "\n",
    "# creating subplots \n",
    "plt.figure(figsize=(10,7))\n",
    "for n, eps in enumerate(eps_values):\n",
    "    ax = plt.subplot(3, 3, n+1)\n",
    "    tdb = test_dbscan2(eps)\n",
    "    sns.scatterplot(pca_df[:, 0], pca_df[:, 1],\n",
    "                    c=tdb, s=30, cmap='rainbow', ax=ax)\n",
    "    ax.set_title('eps='+str(eps)+' n_clusters='+str(len(set(tdb))))\n",
    "    if (n % 3) == 0:\n",
    "        ax.set_ylabel(\"PC2\")\n",
    "    if n>5:\n",
    "        ax.set_xlabel(\"PC1\")\n",
    "plt.tight_layout()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b9ac79e",
   "metadata": {},
   "outputs": [],
   "source": [
    "labels = test_dbscan2(5)\n",
    "\n",
    "ax = sns.scatterplot(pca_df[:, 0], pca_df[:, 1], c=labels, s=30, cmap='rainbow')\n",
    "ax.set_xlabel(\"PC1\")\n",
    "ax.set_ylabel(\"PC2\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2dfb06db",
   "metadata": {},
   "outputs": [],
   "source": [
    "pca = PCA(n_components=2)\n",
    "pca_df = pca.fit_transform(X_train)\n",
    "\n",
    "def test_dbscan2(eps):\n",
    "    dbscan = DBSCAN(\n",
    "    eps=eps, \n",
    "    min_samples=3, \n",
    "    metric='euclidean', \n",
    "    ).fit(pca_df)\n",
    "    \n",
    "    labels = dbscan.labels_\n",
    "    return labels\n",
    "\n",
    "eps_values = np.arange(5, 30, 3)\n",
    "\n",
    "# creating subplots \n",
    "plt.figure(figsize=(10,7))\n",
    "for n, eps in enumerate(eps_values):\n",
    "    ax = plt.subplot(3, 3, n+1)\n",
    "    sns.scatterplot(pca_df[:, 0], pca_df[:, 1], c=test_dbscan(eps), s=30, cmap='rainbow', ax=ax)\n",
    "    ax.set_title('eps='+str(eps)+' n_clusters='+str(len(set(test_dbscan(eps)[1]))))\n",
    "    ax.set_xlabel(\"PC1\")\n",
    "    ax.set_ylabel(\"PC2\")\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "375522cd",
   "metadata": {},
   "source": [
    "# Optional materials\n",
    "Please view the [Youtube video](https://youtu.be/xlqsYUXDuOk) online about OPTICS before you start this part of the tutorial."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee2e6fb9-3658-4f7b-a418-b3ce5ad8a0a0",
   "metadata": {},
   "source": [
    "## Ordering Points To Identify the Clustering Structure (OPTICS)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3105425-5f91-4cda-8126-7c1f1b42ebc6",
   "metadata": {},
   "source": [
    "OPTICS, closely related to DBSCAN, finds core sample of high density and expands clusters from them.  Unlike DBSCAN, OPTICS keeps cluster hierarchy for a variable neighbourhood radius. The sklearn [OPTICS](https://scikit-learn.org/stable/modules/generated/sklearn.cluster.OPTICS.html) function is better suited for usage on large datasets than the current sklearn implementation of DBSCAN. Clusters are then extracted using a DBSCAN-like method (```cluster_method``` = ‘dbscan’) or an automatic technique (```cluster_method``` = ‘xi’). This implementation deviates from the original OPTICS by first performing k-nearest-neighbourhood (kNN, BIDS 5) searches on all points to identify core sizes, then computing only the distances to unprocessed points when constructing the cluster order. Note that we do not employ a heap to manage the expansion candidates, so the time complexity will be O(n^2). The original methodology paper can be found [here](https://dl.acm.org/doi/10.1145/304181.304187). Below are some parameters we will be changing today: \n",
    "\n",
    "- `min_samples`: The number of samples in a neighbourhood for a point to be considered as a core point.\n",
    "- `max_eps`: The maximum distance between two samples for one to be considered as in the neighbourhood of the other. Default value of `np.inf` will identify clusters across all scales; reducing `max_eps` will result in shorter run times.\n",
    "- `metric`: Metric to use for distance computation, default here is '_minkowski_' as opposed to Euclidean in DBSCAN.\n",
    "- `p`: Parameter for the Minkowski metric from pairwise_distances. When p = 1, this is equivalent to using manhattan_distance (l1), and euclidean_distance (l2) for p = 2. For arbitrary p, minkowski_distance (l_p) is used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1285ef95-00ce-4ed2-95d8-f9a5ec811d55",
   "metadata": {},
   "outputs": [],
   "source": [
    "optics = OPTICS(\n",
    "    min_samples=3, \n",
    "    max_eps=np.inf,\n",
    "    p=4\n",
    ").fit(scaled_df.iloc[:, :-1])\n",
    "\n",
    "labels = optics.labels_\n",
    "n_clusters_ = len(set(labels))\n",
    "n_noise_ = list(labels).count(-1)\n",
    "\n",
    "print(\"Estimated number of clusters: %d\" % n_clusters_)\n",
    "print(\"Estimated number of noise points: %d\" % n_noise_)\n",
    "print(\"Homogeneity: %0.3f\" % metrics.homogeneity_score(scaled_df.iloc[:, -1], labels))\n",
    "print(\"Completeness: %0.3f\" % metrics.completeness_score(scaled_df.iloc[:, -1], labels))\n",
    "print(\"Silhouette Coefficient: %0.3f\" % metrics.silhouette_score(scaled_df.iloc[:, :-1], labels))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc7ead2a-c88a-42aa-8f9d-fbcebbe68500",
   "metadata": {},
   "outputs": [],
   "source": [
    "def dim_reduction(method):\n",
    "    df = method.fit_transform(scaled_df.iloc[:, :-1])\n",
    "    sns.scatterplot(df[:, 0], df[:, 1], c=labels, s=30, cmap='Dark2')\n",
    "    plt.title(method)\n",
    "    plt.show()\n",
    "\n",
    "dim_reduction(MDS(n_components=2, metric=True))\n",
    "dim_reduction(TSNE(n_components=2, perplexity=30))\n",
    "dim_reduction(umap.UMAP(n_components=2, n_neighbors=10, min_dist=0.5))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}