BIDS_02_DimensionReduction_1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "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",
   "metadata": {},
   "source": [
    "# Dimensionality Reduction 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Learning objectives\n",
    "1. Apply PCA to different data sets and interpret the output using the popular python library [sklearn](https://scikit-learn.org/stable/)\n",
    "2. Learn how to visualise the model output\n",
    "3. Interpret the results to learn about the data structure and potential outliers\n",
    "4. Code your own function to perform scaling (centering and auto-scaling) using only the numpy.mean and numpy.std functions and two datasets as input (training, test) - use generative AI to help write the code\n",
    "5. Investigate the effect of scaling the data on the output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import specific packages and functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd \n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from sklearn.decomposition import PCA, SparsePCA, KernelPCA\n",
    "from sklearn.preprocessing import StandardScaler, PowerTransformer, RobustScaler"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in the omics datasets using the pandas [read_excel()](https://pandas.pydata.org/docs/reference/api/pandas.read_excel.html) function.\n",
    "For this example we will be using some COVID19 proteomics data.\n",
    "\n",
    "**It is important that the Data folder is placed in the same location as the folders of each of the notebooks.** Otherwise you need to edit the data path in each notebook. Here it goes back to the parent directory of this Notebook and then into the Data-main folder to load the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "covid_proteomics = pd.read_excel(\"../Data-main/COVID19_proteomics.xlsx\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "covid_proteomics.head(20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see the data has three metadata colums: COVID19 (disease status), sample_time (when the blood draw was taken), and sample_id."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Principal Component Analysis (PCA)\n",
    "Read more about PCA in the sklearn [documentation](https://scikit-learn.org/stable/modules/decomposition.html#pca)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the sklearn [PCA()](https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html?highlight=pca#sklearn.decomposition.PCA) function to initialise a PCA object. What is the maximum number of components you can have for a PCA model?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca_res = PCA(n_components=4)\n",
    "# run PCA with 4 components"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apply the `fit_transform()` function to the PCA object to perform PCA dimensionality reduction on the data and project the data to the latent space. Note: we must exclude any sample metadata columns beforehand (remember Python starts counting at 0). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# covid_proteomics.iloc[:, 3:] dataframe without the first 3 metadata columns\n",
    "pca_covid = pca_res.fit_transform(covid_proteomics.iloc[:, 3:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By appling `fit_transform` we project the proteomics data into the latent subspace captured by PCA. The results, also known as **PCA scores** are stored in the results of `fit_transform`, in our case the variable `pca_covid`. Let's look at the PCA results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca_covid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(pca_covid.shape)\n",
    "print(type(pca_covid))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results are returned as a numpy array, in this case of size 382 rows (number of samples) by 4 columns (number of components we selected)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also obtain the PCA components, **also known as eigenvectors**, which represent the influence of each variable (in this case each protein) within each principal component. We do so using the `components_` attribute of the PCA results. There are 4 rows (components) and 450 columns (proteins)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(pca_res.components_)\n",
    "print()\n",
    "print(pca_res.components_.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PCA visualisation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualise two sets of PCA scores against each other with a PCA biplot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set plotting parameters\n",
    "sns.set_style(\"ticks\")\n",
    "sns.set_context(\"notebook\")\n",
    "plt.figure(figsize=(8, 8))\n",
    "\n",
    "# plot a scatterplot using seaborn\n",
    "# the x axis will contain the first column of the pca scores x=pca_covid[:, 0]\n",
    "p = sns.scatterplot(x=pca_covid[:, 0], y=pca_covid[:, 1])\n",
    "\n",
    "p.set_xlabel(\"PC1\")\n",
    "p.set_ylabel(\"PC2\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Colour the scatterpoints on the scatterplot by some metadata - here we will use COVID19 status"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_style(\"ticks\")\n",
    "sns.set_context(\"notebook\")\n",
    "plt.figure(figsize=(8, 8))\n",
    "\n",
    "p = sns.scatterplot(x=pca_covid[:, 0], y=pca_covid[:, 1], hue=covid_proteomics[\"COVID19\"])\n",
    "\n",
    "p.set_xlabel(\"PC1\")\n",
    "p.set_ylabel(\"PC2\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualise multiple components against each other. Here we make use of the seaborn [pairplot()](https://seaborn.pydata.org/generated/seaborn.pairplot.html) function to do so"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# make a dataframe to store the PCA scores for each component alongside the sample metadata\n",
    "pca_df = pd.DataFrame(pca_covid, columns=[\"PC\"+str(i) for i in range(1, pca_covid.shape[1]+1)])\n",
    "pca_df[\"COVID_status\"] = covid_proteomics[\"COVID19\"]\n",
    "\n",
    "pca_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.pairplot(data=pca_df, hue=\"COVID_status\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Let's plot the first two components and add the variance explained to the axes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_style(\"ticks\")\n",
    "sns.set_context(\"notebook\")\n",
    "plt.figure(figsize=(8, 8))\n",
    "\n",
    "p = sns.scatterplot(x=pca_covid[:, 0], y=pca_covid[:, 1], hue=covid_proteomics[\"COVID19\"])\n",
    "\n",
    "# proportions of variance explained by axes\n",
    "pvars = pca_res.explained_variance_ratio_[:2] * 100\n",
    "\n",
    "p.set_xlabel((\"PC1: \" \"{:.2f}%\".format(pvars[0])))\n",
    "p.set_ylabel((\"PC2: \" \"{:.2f}%\".format(pvars[1])))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scree plots show the percentage of the variance in the data explained by each principal component (eigenvalues)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Scree plot\n",
    "\n",
    "# perform PCA with 20 components\n",
    "pca_covid = PCA(n_components=20).fit(covid_proteomics.iloc[:, 3:])\n",
    "\n",
    "# use the attribute .explained_variance_ratio_ to get the eigenvalues\n",
    "variance_per_component = pca_covid.explained_variance_ratio_\n",
    "\n",
    "# sum the eigenvalues to get the cumulative variance explained for each component\n",
    "cumulative_variance = np.cumsum(variance_per_component)\n",
    "components = list(range(1, 21))\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(10, 10))\n",
    "sns.barplot(x=components, y=variance_per_component, palette=\"viridis\", ax=ax)\n",
    "\n",
    "# show the cumulative variance with a blue line\n",
    "sns.pointplot(x=components, y=cumulative_variance, ax=ax, color=\"blue\", label=\"Cumulative variance\")\n",
    "\n",
    "plt.xlabel(\"Components\")\n",
    "plt.ylabel(\"Eigenvalue (% variance explained)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How many components would be required to explain 75% of the variance in the dataset?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpretation of PCA results\n",
    "\n",
    "Looking at the PCA scoreplot and the PCA pairplot, which PCA component(s) show the most clear separation between COVID19 and non-COVID samples?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# extracting the scores and loadings (transposing the latter in the process) of the first 2 components (note: this can be different in other datasets)\n",
    "scores = pca_covid[:, :2]\n",
    "loadings = pca_res.components_[:2].T\n",
    "features = list(covid_proteomics.columns[3:])\n",
    "\n",
    "loadingVector = loadings * np.abs(scores).max(axis=0)\n",
    "\n",
    "sns.set_style(\"ticks\")\n",
    "sns.set_context(\"notebook\")\n",
    "plt.figure(figsize=(8, 8))\n",
    "\n",
    "p = sns.scatterplot(x=pca_covid[:, 0], y=pca_covid[:, 1], hue=covid_proteomics[\"COVID19\"])\n",
    "\n",
    "p.set_xlabel((\"PC1: \" \"{:.2f}%\".format(pvars[0])))\n",
    "p.set_ylabel((\"PC2: \" \"{:.2f}%\".format(pvars[1])))\n",
    "\n",
    "# features as arrows\n",
    "for i, arrow in enumerate(loadingVector):\n",
    "    plt.arrow(0, 0, *arrow, color='k', alpha=0.5, width=0.1, ec='none', length_includes_head=True)\n",
    "    #plt.text(*(arrow * 1.05), features[i], ha='center', va='center') # this is commented, see what happens when you uncomment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly we have too many variables to visualise in a biplot (hence we have switched off labelling them), let's visualise the most important features in a loading plot instead. (You could only visualise the top few variables above, but how do you know how many to select?)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12, 8))\n",
    "\n",
    "lp = sns.lineplot(data=loadings, palette=\"rocket\")\n",
    "\n",
    "lp.set_xlabel(\"Feature\")\n",
    "lp.set_ylabel(\"Loading weights\")\n",
    "lp.set_xticks(range(len(features)), labels=features, rotation=90)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example we might want to select some features that are positive for component 1 and negative for component 2 and those that are negative for component 1 and positive for component 2. As these vectors are in the rough direction of the class separation, be aware that supervised methods (BIDS 7-10) aim to do this as they find the best way to use the features to predict the classes. There can be other reasons to inspect this in an unsupervised analysis such as possible features that are outliers or those that cause groupings that are unrelated to any class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loading_vector = loadings[:,0]*loadings[:,1]\n",
    "n = 20\n",
    "ind = loading_vector.argsort()[:n]\n",
    "features = np.array(features)\n",
    "top_features = features[ind]\n",
    "top_features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loadingVector_n = loadingVector[ind,:]\n",
    "\n",
    "sns.set_style(\"ticks\")\n",
    "sns.set_context(\"notebook\")\n",
    "plt.figure(figsize=(8, 8))\n",
    "\n",
    "p = sns.scatterplot(x=pca_covid[:, 0], y=pca_covid[:, 1], hue=covid_proteomics[\"COVID19\"])\n",
    "\n",
    "p.set_xlabel((\"PC1: \" \"{:.2f}%\".format(pvars[0])))\n",
    "p.set_ylabel((\"PC2: \" \"{:.2f}%\".format(pvars[1])))\n",
    "\n",
    "# features as arrows\n",
    "for i, arrow in enumerate(loadingVector_n):\n",
    "    plt.arrow(0, 0, *arrow, color='k', alpha=0.5, width=0.1, ec='none', length_includes_head=True)\n",
    "    plt.text(*(arrow * 1.05), top_features[i], ha='center', va='center')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12, 8))\n",
    "\n",
    "lp = sns.lineplot(data=loadings[ind,:], palette=\"rocket\")\n",
    "\n",
    "lp.set_xlabel(\"Feature\")\n",
    "lp.set_ylabel(\"Loading weights\")\n",
    "lp.set_xticks(range(len(top_features)), labels=top_features, rotation=90)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Identifying outliers\n",
    "Samples can be labelled in order to identify those that are outlying."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca_covid = PCA(n_components=2).fit_transform(covid_proteomics.iloc[:, 3:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,8))\n",
    "sns.scatterplot(x=pca_covid[:, 0], y=pca_covid[:, 1], hue=covid_proteomics[\"COVID19\"])\n",
    "\n",
    "# for loop to add labels to each x, y pair along with the corresponding sample ID\n",
    "for i in range(pca_covid.shape[0]):\n",
    "    plt.text(x=pca_covid[:, 0][i]+0.3, y=pca_covid[:, 1][i]+0.3, s=covid_proteomics[\"sample_id\"][i], \n",
    "          fontdict=dict(color='black',size=8))\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Are the any obvious outliers in this dataset as seen from this PCA biplot?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scaling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "The standard score of a sample $x$ is defined as:\n",
    "\n",
    "\n",
    "$$z = \\frac{(x-\\mu)}{\\sigma}$$\n",
    "\n",
    "Where:\n",
    "- $\\mu$ is the feature mean\n",
    "- $\\sigma$ is the feature standard deviation\n",
    "\n",
    "\n",
    "\n",
    "Code a function to scale the data such that each feature has a:\n",
    "- Mean of 0 \n",
    "- Standard deviation of 1\n",
    "\n",
    "Use only the numpy [mean()](https://numpy.org/doc/stable/reference/generated/numpy.mean.html) and numpy [std()](https://numpy.org/doc/stable/reference/generated/numpy.std.html) functions\n",
    "\n",
    "You can check your answer by comparing it to the result achieved using the sklearn [StandardScaler()](https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.StandardScaler.html#:~:text=Standardize%20features%20by%20removing%20the%20mean%20and%20scaling%20to%20unit%20variance.&text=where%20u%20is%20the%20mean,or%20one%20if%20with_std%3DFalse%20.) function\n",
    "\n",
    "Use a generative AI ([ChatGPT](https://chat.openai.com), [Code-LLaMa](https://poe.com/Code-Llama-34b), etc) of your choice and prompt it step by step to obtain the code. Ensure you only allow it to use numpy functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def my_scaler(data):\n",
    "    \n",
    "    # get mean of each column\n",
    "\n",
    "    # subtract the column mean from each value in the column\n",
    "\n",
    "    # check the mean of each column is now 0\n",
    "\n",
    "    # get the standard deviations of each column\n",
    "\n",
    "    # divide the mean centered data by the standard deviation to scale by unit variance\n",
    "\n",
    "    return scaled_data\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# apply the scaler function to the data\n",
    "covid_proteomics_scaled = my_scaler(covid_proteomics.iloc[:, 3:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# perform PCA using the scaled data\n",
    "pca_covid_scaled = PCA(n_components=2).fit_transform(covid_proteomics_scaled)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How does scaling affect the PCA? visualise the output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5))\n",
    "\n",
    "sns.scatterplot(x=pca_covid[:, 0], y=pca_covid[:, 1], hue=covid_proteomics[\"COVID19\"], ax=ax1)\n",
    "ax1.set_title(\"No scaling\")\n",
    "\n",
    "sns.scatterplot(x=pca_covid_scaled[:, 0], y=pca_covid_scaled[:, 1], hue=covid_proteomics[\"COVID19\"], ax=ax2)\n",
    "ax2.set_title(\"Standard scaling\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "sklearn also has a number of inbuilt functions for scaling. Look into the following functions and apply them on the dataset:\n",
    "- [sklearn.StandardScaler()](https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.StandardScaler.html#:~:text=Standardize%20features%20by%20removing%20the%20mean%20and%20scaling%20to%20unit%20variance.&text=where%20u%20is%20the%20mean,or%20one%20if%20with_std%3DFalse%20.)\n",
    "- [sklearn.RobustScaler()](https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.RobustScaler.html?highlight=robust%20scaler#sklearn.preprocessing.RobustScaler)\n",
    "- [sklearn.PowerTransformer()](https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.PowerTransformer.html?highlight=power%20transformer#sklearn.preprocessing.PowerTransformer)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's begin with appling the PowerTransformer():"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# initialise PowerTransformer object\n",
    "pt = PowerTransformer()\n",
    "\n",
    "# Apply the power transformation to the data\n",
    "covid_proteomics_power_transform = pt.fit_transform(covid_proteomics.iloc[:, 3:])\n",
    "\n",
    "# Apply PCA to the power transformed data\n",
    "pca_covid_power_transform = PCA(n_components=2).fit_transform(covid_proteomics_power_transform)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Complete the following code for the Standard Scaler and Robust Scaler:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Standard scaler\n",
    "ss = StandardScaler()\n",
    "\n",
    "covid_proteomics_standard_scaled = \n",
    "\n",
    "pca_covid_standard_scaled = \n",
    "\n",
    "# Robust scaler\n",
    "rs = \n",
    "\n",
    "covid_proteomics_robust_scaled = \n",
    "\n",
    "pca_covid_robust_scaled = "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualise the results of the scaling using PCA biplots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 5))\n",
    "\n",
    "sns.scatterplot(x=pca_covid_power_transform[:, 0], y=pca_covid_power_transform[:, 1], hue=covid_proteomics[\"COVID19\"], ax=ax1)\n",
    "ax1.set_title(\"Power transformer\")\n",
    "\n",
    "sns.scatterplot(x=pca_covid_standard_scaled[:, 0], y=pca_covid_standard_scaled[:, 1], hue=covid_proteomics[\"COVID19\"], ax=ax2)\n",
    "ax2.set_title(\"Standard scaling\")\n",
    "\n",
    "sns.scatterplot(x=pca_covid_robust_scaled[:, 0], y=pca_covid_robust_scaled[:, 1], hue=covid_proteomics[\"COVID19\"], ax=ax3)\n",
    "ax3.set_title(\"Robust scaling\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Outlier detection on scaled data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now use a generative AI to create the code for calculating Hotelling's T<sup>2</sup> statistic and visualise the outliers (if any) using the [pairplot()](https://seaborn.pydata.org/generated/seaborn.pairplot.html) function as before.\n",
    "\n",
    "Example prompt:\n",
    "\n",
    "_Given these packages:_ <br>\n",
    "_import pandas as pd_ <br>\n",
    "_import numpy as np_ <br>\n",
    "_import matplotlib.pyplot as plt_ <br>\n",
    "_import seaborn as sns_ <br>\n",
    "_from sklearn.decomposition import PCA, SparsePCA, KernelPCA_ <br>\n",
    "_from sklearn.preprocessing import StandardScaler, PowerTransformer, RobustScaler_ <br>\n",
    "_Can you write the code to calculate the Hotelling's T2 statistic to find outliers based on the model below:_ <br>\n",
    "_pca_covid = PCA(n_components=7).fit_transform(covid_proteomics_robust_scaled)_ <br>\n",
    "_Then visualise any outliers in a pair plot_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the above is correct, now see if this can be rewritten in a function you can reuse later on:\n",
    "\n",
    "Example prompt: _Can you rewrite this as a function with the pca_df and percentile (default as 95) as inputs and outliers as output?_\n",
    "\n",
    "Sometimes they will get it right the first time, but other times it helps to break it down in parts first before continuing with other steps (visualising, write as function, etc)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add the function here, rewrite it in your preferred style\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# visualise it and compare that it produced the same output as above when it was not a function\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optional materials\n",
    "\n",
    "## PCA extensions: sparse PCA and kernel PCA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sparse PCA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sparse PCA is performed using the [SparsePCA()](https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.SparsePCA.html) function. The level of sparsity (the proportion of input variables contributing to the principal components) is controllable by the coefficient of the L1 penalty, given by the parameter `alpha`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some algorithms make use of a random process, e.g. to initialise weights it select initial starting points from a random number generator. This means that each time you use the algorithm, the answer might be slightly (sometimes this is unnoticeable) different. Also, this may mean the output you get is different from your neighbour.\n",
    "\n",
    "How do you get the same answer each time (as long as you do not change any other parameters)? You can set the random state to a specific number.\n",
    "\n",
    "Random number generators use the computer time and other factors to start at a random place in a huge sequence of random numbers. By setting the ```random_state``` property of an algorithm you will ensure it will start at a standard place. You can pick any positive integer for this. If you look at the documentation of [PCA](https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html?highlight=pca#sklearn.decomposition.PCA) and [SparsePCA](https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.SparsePCA.html#sklearn.decomposition.SparsePCA) you will see both of these methods have this argument as input, with a default value, however using other default parameters PCA will give the same output whereas the output from SparsePCA might change."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# enter your CID here, or date of birth, or another number of your choosing to use as random state\n",
    "# remember to check the documentation of each algorithm if setting the random_state is needed\n",
    "CID = 0 # do not add the leading 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sparse_pca = SparsePCA(n_components=2, alpha=1, random_state=CID)\n",
    "sparse_pca_covid = sparse_pca.fit_transform(covid_proteomics.iloc[:, 3:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look at the mean number of 0 values across the components (sparsity level). How does this change when you change the alpha parameter? \n",
    "\n",
    "If `alpha` = 0 there is no sparsity constraint, and all input variables will contribute to the principal components. The higher alpha is, the less variables will contribute to the principal components. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(sparse_pca.components_ == 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sparse_pca_1 = SparsePCA(n_components=2, alpha=1).fit_transform(covid_proteomics.iloc[:, 3:])\n",
    "sparse_pca_10 = SparsePCA(n_components=2, alpha=10).fit_transform(covid_proteomics.iloc[:, 3:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5))\n",
    "\n",
    "sns.scatterplot(x=sparse_pca_1[:, 0], y=sparse_pca_1[:, 1], hue=covid_proteomics[\"COVID19\"], ax=ax1)\n",
    "ax1.set_title(\"alpha=1\")\n",
    "\n",
    "sns.scatterplot(x=sparse_pca_10[:, 0], y=sparse_pca_10[:, 1], hue=covid_proteomics[\"COVID19\"], ax=ax2)\n",
    "ax2.set_title(\"alpha=10\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Kernel PCA\n",
    "We will be discussing more about kernels in BIDS 8, so you can revisit this section after that session.\n",
    "\n",
    "[Kernel PCA](https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.KernelPCA.html#sklearn.decomposition.KernelPCA) is a form of non-linear dimensionality reduction using kernels. There are several hyperparameters that can be tuned for this model, the main being the type of kernel used. The sklearn KernelPCA() function supports the following kernels: 'linear’, ‘poly’, ‘rbf’, ‘sigmoid’, ‘cosine’, ‘precomputed’, with the default being ’linear’."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# no kernel has been specified so it uses linear by default\n",
    "kernel_pca = KernelPCA(n_components=2)\n",
    "kernel_pca_covid = kernel_pca.fit_transform(covid_proteomics.iloc[:, 3:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now try to use another kernel type:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kernel_pca_rbf = KernelPCA(n_components=2, kernel=\"rbf\").fit_transform(covid_proteomics.iloc[:, 3:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Experiment with the other parameters listed on the reference page, such as Gamma - the kernel bandwidth parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5))\n",
    "\n",
    "sns.scatterplot(x=kernel_pca_covid[:, 0], y=kernel_pca_covid[:, 1], hue=covid_proteomics[\"COVID19\"], ax=ax1)\n",
    "ax1.set_title(\"Linear kernel\")\n",
    "\n",
    "sns.scatterplot(x=kernel_pca_rbf[:, 0], y=kernel_pca_rbf[:, 1], hue=covid_proteomics[\"COVID19\"], ax=ax2)\n",
    "ax2.set_title(\"Radial basis function kernel\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which kernel type do you think provides the best separation for this dataset?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualise the results of the COVID19 dataset using different PCA types:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 5))\n",
    "\n",
    "sns.scatterplot(x=pca_covid[:, 0], y=pca_covid[:, 1], hue=covid_proteomics[\"COVID19\"], ax=ax1)\n",
    "ax1.set_title(\"PCA\")\n",
    "\n",
    "sns.scatterplot(x=sparse_pca_covid[:, 0], y=sparse_pca_covid[:, 1], hue=covid_proteomics[\"COVID19\"], ax=ax2)\n",
    "ax2.set_title(\"Sparse PCA\")\n",
    "\n",
    "# if we use linear kernel PCA this will be the same result as standard PCA\n",
    "sns.scatterplot(x=kernel_pca_covid[:, 0], y=kernel_pca_covid[:, 1], hue=covid_proteomics[\"COVID19\"], ax=ax3)\n",
    "ax3.set_title(\"Kernel PCA\")\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Your turn\n",
    "Select another dataset from the `Data` folder and import it using the pandas `read_excel()` function as above. Scale the data, and then apply PCA (standard, sparse, or kernel) and visualise the results. Reuse the codes above for a new dataset, copy and paste them below. Reusing existing code (and sometimes slightly adapting it) is common practice, and you will do this as well for the final assignment where you will reuse the materials from the tutorials for the assignment.\n",
    "\n",
    "What can you interpret from the PCA results? Are there any outliers?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import datasets..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform scaling..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform PCA..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualise PCA results and detect outliers"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "b36f82fc6ec7da3c99303f5c6c1c64ad749da77112fc70d0beec9ca6a90b9819"
  },
  "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": 2
}