Multivariate Analysis PCA-NMR.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "n8itkL-uC2jG"
      },
      "source": [
        "# Multivariate Analysis - PCA"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "aVRBt5QBC2jH"
      },
      "source": [
        "In this notebook we will perform an unsupervised *multivariate* exploratory analysis with principal component analysis (PCA) of 4 racks from the NMR profiling dataset of the *Dementia research cohort*. For details of the study see the Metabolights Study [MTBLS719](https://www.ebi.ac.uk/metabolights/MTBLS719).\n",
        "\n",
        "The notebook is divided in the following steps:\n",
        "\n",
        "1) Model fitting basics: Fit PCA models to the dataset with different scaling options.\n",
        "\n",
        "2) Model Cross-validation and component selection: Describe model cross-validation routines, and best practices for model performance benchmarking and parameter selection.\n",
        "\n",
        "3) Outlier detection and model interpretation: Use PCA to explore the main trends in the dataset and detect potential outliers.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rgAwkmFbC2jI"
      },
      "source": [
        "## Code import\n",
        "\n",
        "Import all the packages and configure notebook plotting mode."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "9JiGMbjq38j7"
      },
      "outputs": [],
      "source": [
        "!pip install -q ipympl\n",
        "!pip install -q kneed"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Jsc30bf7C9x6"
      },
      "outputs": [],
      "source": [
        "!git clone https://github.com/IPTC-DataAnalysisCourse/chemometrics-tutorials.git"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "TCUagVyxFDhc"
      },
      "outputs": [],
      "source": [
        "%cd chemometrics-tutorials"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "UodohlppC2jI"
      },
      "outputs": [],
      "source": [
        "# Import the required python packages including\n",
        "# the custom Chemometric Model objects\n",
        "import numpy as np\n",
        "from sklearn import preprocessing\n",
        "import pandas as pds\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "from pyChemometrics.ChemometricsPCA import ChemometricsPCA\n",
        "from pyChemometrics.ChemometricsScaler import ChemometricsScaler\n",
        "from pyChemometrics.plotting_utils import plotLoadings\n",
        "\n",
        "# Use to obtain same values as in the text\n",
        "np.random.seed(350)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "vs4qMJwUaL-j"
      },
      "outputs": [],
      "source": [
        "import sklearn\n",
        "print('The scikit-learn version is {}.'.format(sklearn.__version__))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "njODRlceC2jJ"
      },
      "source": [
        "The next cell sets up the figure display mode. The *notebook* and *ipympl* modes allows interactive plotting. Another option is to select *inline*, to obtain static plots in a notebook cell."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "7nY_p-QMaL-k"
      },
      "outputs": [],
      "source": [
        "# Set the plot backend to support interactive plotting\n",
        "#%matplotlib notebook # run this line when running in jupyter notebook"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "H0CuKOELC2jJ"
      },
      "outputs": [],
      "source": [
        "# Set the plot backend to support interactive plotting - Run this when running in Google Colab\n",
        "%matplotlib ipympl\n",
        "\n",
        "import plotly.io as pio\n",
        "pio.renderers.default = \"colab\"\n",
        "\n",
        "from google.colab import output\n",
        "output.enable_custom_widget_manager()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2JtQ2G0qC2jJ"
      },
      "source": [
        "## Data import"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_7tTlTwhC2jK"
      },
      "source": [
        "We will now import the NMR data and the metadata (Y variable).\n",
        "\n",
        "X - NMR data matrix\n",
        "\n",
        "Y - Matrix with one metadata outcome\n",
        "\n",
        "ppm - Chemical shift axis for the NMR data in H $\\delta$ppm.\n",
        "\n",
        "#### Metadata\n",
        "Y - represents the sex of the individuals (1: Male, 0: Female, in original Y data matrix)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "2U7cC1_zC2jK"
      },
      "outputs": [],
      "source": [
        "# Load the dataset\n",
        "dementia_nmr_dataset = pds.read_csv(\"./Data/Dementia U NMR_PQN_normalised_Data.csv\",delimiter=',')\n",
        "\n",
        "#dementia_nmr_dataset = pds.read_csv(\"./Data/Dementia U NMR_Data.csv\",delimiter=',')\n",
        "\n",
        "# Delete samples where outcome variable is unknown - Study Samples in standard NPC pipeline\n",
        "dementia_nmr_dataset = dementia_nmr_dataset[~dementia_nmr_dataset['Gender'].isnull()]\n",
        "\n",
        "# Create the X matrix\n",
        "X = dementia_nmr_dataset.iloc[:, 5::].values\n",
        "\n",
        "# Use pandas Categorical type to generate the dummy enconding of the Y vector (0 and 1)\n",
        "Y = pds.Categorical(dementia_nmr_dataset['Gender']).codes\n",
        "\n",
        "ppm = np.array(dementia_nmr_dataset.columns[5::], dtype=np.float64)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7YNe7koVC2jL"
      },
      "source": [
        "**Note**: To apply the analyses exemplified in this notebook to any other dataset, just modify the cell above to import the data matrices and vectors X and Y from any other source file.\n",
        "\n",
        "The expected data types and formatting for **X** and **Y** are:\n",
        "\n",
        "   **X**: Any data matrix with n rows (observations/samples) and p columns (variables/features). The matrix should be provided as a [numpy.ndarray](https://numpy.org/doc/stable/reference/generated/numpy.ndarray.html) object, with 2 dimensions, and with shape = (n, p). We recommend using the *numpy* function [numpy.genfromtxt](https://numpy.org/devdocs/reference/generated/numpy.genfromtxt.html) or the *pandas* [pandas.read_csv](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.read_csv.html) function to read the data from a text file. When using the *pandas.read_csv* function, extract the data matrix as a *numpy.ndarray* from the pandas.DataFrame object using the `.values` attribute.\n",
        "```\n",
        "X_DataFrame = pds.read_csv(\"./data/X_spectra.csv\")\n",
        "X = X_DataFrame.values\n",
        "```\n",
        "   \n",
        "   **Y** vectors: Each **Y** vector should be a 1-dimensional [numpy.ndarray](https://numpy.org/doc/stable/reference/generated/numpy.ndarray.html) object, with a number and ordering of elements matching the rows in **X**. For continuous variables, any regular *numpy.ndarray* with a data type of `int` (integers only) or `float` can be used.\n",
        "   ```\n",
        "   Y_continuous = numpy.ndarray([23.4, 24, 0.3, -1.23], dtype='float')\n",
        "   ```\n",
        "To encode binary class labels, a *numpy.ndarray* of dtype `int`, with 0 and 1 as labels (e.g., 0 = Control, 1 = Case) must be used. The way in which classes are encoded will affect the model interpretation: the class labeled as 1 is used as the \"positive/case\" class by the *pyChemometrics* objects.\n",
        "   \n",
        "   In the example above, we used the *pandas* [Categorical](https://pandas.pydata.org/pandas-docs/stable/user_guide/categorical.html) datatype to handle the conversion of the original numerical values (1, 2) to the required (0, 1) labels. After converting a column to a `Categorical` datatype, the `.codes` attribute returns a vector with the same length of the original Y, but where each value is replaced by their integer (`int`) code. The correspondence between code and category can be inspected with the `categories` attribute. The order of the labels in `.codes` is the same as the order of the `categories` attribute (i.e. 0 is the first element in `categories`, 1 the second and so on).\n",
        "   ```\n",
        "   Y = pds.Categorical(Y.iloc[:, 1])\n",
        "   Y.codes # The numerical label\n",
        "   Y.categories # Original text or numerical description of the category\n",
        "   ```\n",
        "   [get_dummies](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.get_dummies.html) is another helpful function to perform dummy (0-1) encoding of variables."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "h11JHelyC2jL"
      },
      "source": [
        "Plot the spectra in the dataset.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "_DdLx_vyC2jL"
      },
      "outputs": [],
      "source": [
        "# Plot the spectra in the dataset\n",
        "plt.figure()\n",
        "plt.plot(ppm, X.T)\n",
        "plt.title(\"X matrix of spectra\")\n",
        "plt.xlabel(\"$\\delta$ppm 1H\")\n",
        "plt.gca().invert_xaxis()\n",
        "plt.ylabel(\"Intensity\")\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NbdDhZ18C2jM"
      },
      "source": [
        "# Exploratory Unsupervised analysis using PCA"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BZfIkaj3C2jM"
      },
      "source": [
        "In this tutorial we will fit a PCA model to explore the general trends in the dataset and assess if there are any potential outliers (cause by instrumental issues during data acquisition, for example). We start by describing the model syntax for fitting the PCA models."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PYHQbaTxC2jM"
      },
      "source": [
        "## 1) PCA model fitting and scaling"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "SPZdny1qC2jM"
      },
      "source": [
        "### Scaling options and preliminary model fitting\n",
        "\n",
        "We will start by calculating a series of PCA models with 2 components, each with one of the 3 common scaling choices in chemometrics - mean centring (MC), Unit Variance (UV) and Pareto (Par) scaling. The choice of components to use in the modeling will be addressed properly in the next section."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "5v-4CmVpC2jM"
      },
      "outputs": [],
      "source": [
        "# Select the scaling options:\n",
        "# Here we are generating 3 scaling objects to explore the effect of scaling in PCA:\n",
        "\n",
        "# Unit-Variance (UV) scaling:\n",
        "scaling_object_uv = ChemometricsScaler(scale_power=1)\n",
        "\n",
        "# Mean Centering (MC):\n",
        "scaling_object_mc = ChemometricsScaler(scale_power=0)\n",
        "\n",
        "# Pareto scaling (Par):\n",
        "scaling_object_par = ChemometricsScaler(scale_power=0.5)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "u0HvDLpfC2jM"
      },
      "source": [
        "The scaling object will store the vector of column means and standard deviations as it was estimated from the dataset passed to the *fit* method (i.e., during \"training\" of the classifier)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BzgPETuXC2jM"
      },
      "source": [
        "Pass each scaling object to the PCA method."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "l8LKDbhRC2jN"
      },
      "outputs": [],
      "source": [
        "# Create and fit the PCA model - starting with UV\n",
        "PCA_model_uv = ChemometricsPCA(ncomps=2, scaler=scaling_object_uv)\n",
        "PCA_model_uv.fit(X)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "v6DKj00TC2jN"
      },
      "outputs": [],
      "source": [
        "# Create and fit the PCA model - MC\n",
        "PCA_model_mc = ChemometricsPCA(ncomps=2, scaler=scaling_object_mc)\n",
        "PCA_model_mc.fit(X)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "-Kr_7fW8C2jN"
      },
      "outputs": [],
      "source": [
        "# Create and fit the PCA model - Par\n",
        "PCA_model_par = ChemometricsPCA(ncomps=2, scaler=scaling_object_par)\n",
        "PCA_model_par.fit(X)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "i0beuFycC2jN"
      },
      "source": [
        "### Effect of scaling on PCA Score plots"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Z8sqt_OKC2jN"
      },
      "source": [
        "These plots show the effect that different scaling parameters have on the PCA scores. The scores plots are a usefull summary of the multivariate similarity between samples, and will be used to inspect the main trends in the data."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "H2kzqCkFC2jN",
        "scrolled": false
      },
      "outputs": [],
      "source": [
        "# PCA score plot for the mean centered model\n",
        "PCA_model_mc.plot_scores(comps=[0, 1], label_outliers=True, plot_title='Mean centering')\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "tqwu-xUoC2jO"
      },
      "outputs": [],
      "source": [
        "# Score plot for the Pareto scaled model\n",
        "PCA_model_par.plot_scores(comps=[0, 1], label_outliers=True, plot_title='Pareto scaling')\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "71NJVZPrC2jO"
      },
      "outputs": [],
      "source": [
        "# PCA score plot for UV scaled model\n",
        "PCA_model_uv.plot_scores(comps=[0, 1],label_outliers=True, plot_title='UV scaling')\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "J3ns4NmaC2jO"
      },
      "source": [
        "### Effect of scaling on PCA loadings"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "y2fgVcO2C2jO"
      },
      "source": [
        "These plots show the effect that different scaling parameters have on the PCA loadings (also designated as $p$, when refering to a single loading vector or $P$ to the matrix containing the loading vectors for all components).\n",
        "\n",
        "Although different types of scaling can be used to investigate the structure of the dataset, each type of scaling results in different models, potentially with different interpretations. The *plot_scores* function automatically draws a 95% confidence Hotelling $T^{2}$ ellipse, and flags the points as potential outliers."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "2BWfXHtZC2jO"
      },
      "outputs": [],
      "source": [
        "# Plot of first principal component loadings of mean centering model\n",
        "ax = PCA_model_mc.plot_loadings(ncomp=1, xaxis=ppm)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "zjFOvzy7C2jO"
      },
      "outputs": [],
      "source": [
        "# Plot of first principal component loadings of Pareto scaled model\n",
        "ax = PCA_model_par.plot_loadings(ncomp=1, xaxis=ppm)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "f5cr4iqgC2jO"
      },
      "outputs": [],
      "source": [
        "# Plot of first principal component loadings of Unit Variance scaled model\n",
        "ax = PCA_model_uv.plot_loadings(ncomp=1, xaxis=ppm)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "p4e3d9juC2jO"
      },
      "source": [
        "This example highlights the effects of different scaling approaches on PCA scores, loadings, their interpretation and the structure recovered by the model.\n",
        "\n",
        "We will now proceed with the exploratory analysis of this dataset using UV scaling."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0YfFvDqZC2jP"
      },
      "source": [
        "Direct plots of the loading vectors can be hard to interpret, especially when using UV scaling. A better way to visualize these quantities is to overlay them on a representative spectrum (e.g. median spectrum) with a colourscale."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Ivc8vu7VC2jP"
      },
      "outputs": [],
      "source": [
        "# Plot of first principal component loadings overlaid on the median spectrum - UV\n",
        "plotLoadings(PCA_model_uv.loadings[0, :], ppm, spectra=X, cbarlabels ='Loadings')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "I2-mSeJvC2jP"
      },
      "outputs": [],
      "source": [
        "# Plot of first principal component loadings overlaid on the median spectrum - Mean centring\n",
        "plotLoadings(PCA_model_mc.loadings[0, :], ppm, spectra=X, cbarlabels = \"Loadings\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "f6sLbjjnC2jP"
      },
      "source": [
        "## 2) Model cross-validation and component selection"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8ifd_CCCC2jP"
      },
      "source": [
        "When generating a PCA model, the number of components is the main parameter that must be chosen.\n",
        "\n",
        "Ideally, we want to select enough components to capture as much structured variation in the data as possible, but not so  many that random noise starts to be also incorporated in the PCA model.\n",
        "\n",
        "A sugestion to select the number of components is to use the $Q^{2}X$ measure, and pick the number of components after this metric reaches a plateau (e.g. less than 5% increase compared to previous number of components)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "8dcQik8KC2jP"
      },
      "outputs": [],
      "source": [
        "PCA_model_uv.scree_plot(X, total_comps=10)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3cA6lQqTC2jW"
      },
      "source": [
        "This suggestion should suffice for exploratory data analysis with PCA. However, the $Q^{2}X$ measure obtained for K-Fold cross validation is sensitive to row permutation of the X matrix. A more robust alternative is to use repeated cross-validation, shuffle the rows in the X matrix each time, and see the distribution of $Q^{2}X$ values per component. This should give a more comprehensive overview of each component's robustness."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qQ1wsNPCC2jW"
      },
      "source": [
        "**Note**: In an exploratory PCA analysis, the selection of the number of components is not as critical as in a PLS-DA model, and the user can simply select a set of components to explore the data and interactively adjust. Nevertheless, its still important to benchmark the cross-validated performance of the model, to avoid interpreting or infering biological conclusions from non-robust principal components."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "uyC6mtZlC2jW"
      },
      "outputs": [],
      "source": [
        "rep_cv = PCA_model_uv.repeated_cv(X, repeats=5)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "npF9CbXsC2jW"
      },
      "source": [
        "Refit the model for further exploration with the selected number of components."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Q3s8kEEFC2jX"
      },
      "outputs": [],
      "source": [
        "# Create and fit the PCA model - UV scaling\n",
        "PCA_model_uv = ChemometricsPCA(ncomps=2, scaler=scaling_object_uv)\n",
        "PCA_model_uv.fit(X)\n",
        "PCA_model_uv.cross_validation(X)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zlrP5s4HC2jX"
      },
      "source": [
        "If the *cross_validation* method has been called before for a PCA model object, the plot_loadings method will also show the estimated confidence bands (in light red) for the loading parameters. These were estimated during the cross validation procedure. The default value of $sigma = 2$ means that the plot uses the average $\\pm2\\sigma$ to draw the confidence bands."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "TC_XQyBcaL-w"
      },
      "outputs": [],
      "source": [
        "ax = PCA_model_uv.plot_loadings(ncomp=1, sigma=2, xaxis=ppm, instrument = 'nmr')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PDVGxYadC2jX"
      },
      "source": [
        "## 3) Outlier detection and model interpretation"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VNw-vI0gC2jX"
      },
      "source": [
        "### Outlier detection\n",
        "\n",
        "PCA can be used to detect potential outliers, and screen samples where a potential problems might have occurred during data acquisition. The main outlier detection tool is the Hotelling $T^{2}$ statistic, a multivariate generalization of the Student's t-distribution. The *plot_scores* function automatically draws a 95% confidence ellipse for $T^{2}$ in the score plot. Samples outside the ellipse are candidate outliers and warrant further investigation.\n",
        "\n",
        "**Note**: Bear in mind that outlier exclusion can affect the model performance, and therefore every outlier removed during analysis should be justified and always recorded.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "QL4gcFuuC2jX"
      },
      "outputs": [],
      "source": [
        "PCA_model_uv.plot_scores(comps=[0, 1], label_outliers=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hKlOp0jeC2jX"
      },
      "source": [
        "The PCA score plot highlighted a set of candidate outliers for the 2 principal components model. There are some outliers in component 1 and 2.\n",
        "\n",
        "The *.outlier* function can be used to automatically obtain the indices for candidate outlying samples using all or only a set of selected components. This synthax can be used to obtain outliers matching any given visualization or alongside a single PC component, or to obtain an assessment of which samples are potential outliers using the whole model.\n",
        "\n",
        "If using Hotelling T2, the significance level can be adjusted to control for the proportion of false positives in outlier detection."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "696bfSyZC2jY"
      },
      "outputs": [],
      "source": [
        "outlier_idx = PCA_model_uv.outlier(X)\n",
        "print(\"Outliers for the full 2 component model : {0}\".format(outlier_idx))\n",
        "outlier_idx = PCA_model_uv.outlier(X, comps=[0])\n",
        "print(\"Outliers for the 1st principal component : {0}\".format(outlier_idx))\n",
        "outlier_idx = PCA_model_uv.outlier(X, comps=[1])\n",
        "print(\"Outliers for the 2nd principal component : {0}\".format(outlier_idx))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ttiADcCrC2jY"
      },
      "source": [
        "The next step is to identify the reason for the outlying scores. This can be done by\n",
        "inspection of the raw data (or a subset of it), analysis of the loading vectors, and by comparing the model predictions/data reconstruction performed by the model. One of the outliers has a very large peak close to 2.15 ppm as well some other large peaks in the aliphatic (0.8-4 ppm) and aromatic regions (6.75-8 ppm)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Xf1oSVq_aL-y"
      },
      "outputs": [],
      "source": [
        "# plot the mean spectrum calculated from the raw data (blue) and the outlying spectra (red)\n",
        "plt.figure()\n",
        "plt.plot(ppm, X[outlier_idx, :].T, 'r')\n",
        "plt.plot(ppm, np.mean(X, axis=0), 'b')\n",
        "plt.gca().invert_xaxis()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wNyEv5V8aL-y"
      },
      "source": [
        "As all the highligted samples are not extreme outliers therefore these need to examined separately by plotting only the samples which appear as extreme outliers in the PCA scores plot"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "djgl7dCUaL-y"
      },
      "outputs": [],
      "source": [
        "outlier_idx_extreme = np.array([169, 263, 283])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "NBX8OE4EC2jY"
      },
      "outputs": [],
      "source": [
        "# plot the mean spectrum calculated from the raw data (blue) and the outlying spectra (red)\n",
        "plt.figure()\n",
        "plt.plot(ppm, X[outlier_idx_extreme, :].T, 'r')\n",
        "plt.plot(ppm, np.mean(X, axis=0), 'b')\n",
        "plt.gca().invert_xaxis()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "cetd5KhFC2jY"
      },
      "source": [
        "The PCA score plot should be interpreted using the loading vectors for the corresponding components."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "YkrrGANoC2jY"
      },
      "outputs": [],
      "source": [
        "ax = PCA_model_uv.plot_loadings(ncomp=1, xaxis=ppm, instrument = 'nmr')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-wskRtLEC2jY"
      },
      "source": [
        "Direct inspection of the loadings plot obtained from the UV scaled model is not very straightforward...\n",
        "\n",
        "Another way to investigate outliers is to use score values (from actual samples or artificially created) to generate reconstructed spectra. We can explore the PCA model space in this way, and compare spectra representative of any desired region of the score plot.\n",
        "We backtransform the reconstructed samples obtained to the original dataspace (reverse the scaling and mean centering) for plotting and comparison."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "699Bea-sC2jY"
      },
      "outputs": [],
      "source": [
        "# For all outliers in PC2\n",
        "# Use the center of the model as control\n",
        "model_center_sample = PCA_model_uv.inverse_transform([0, 0])\n",
        "# Reconstruct spectra from all the outliers in PC2\n",
        "outlier_idx = PCA_model_uv.outlier(X, comps=[1])\n",
        "out_scores = PCA_model_uv.scores[outlier_idx, :]\n",
        "outliers = PCA_model_uv.inverse_transform(out_scores)\n",
        "# Reconstruct a spectrum for the \"mean\" of these outliers\n",
        "mean_outlier = PCA_model_uv.inverse_transform(out_scores.mean(axis=0))\n",
        "\n",
        "plt.figure()\n",
        "# \"center\" representative \"normal\" sample plotted in blue\n",
        "plt.plot(ppm, model_center_sample, 'b')\n",
        "# The outliers plotted in dashed red line\n",
        "plt.plot(ppm, outliers.T, 'r--',)\n",
        "# The mean outlier plotted in green\n",
        "plt.plot(ppm, mean_outlier, 'g')\n",
        "\n",
        "plt.gca().invert_xaxis()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "42G9z1y9aL-z"
      },
      "outputs": [],
      "source": [
        "# For extreme outliers from PC2 and the overall model\n",
        "# Use the center of the model as control\n",
        "model_center_sample = PCA_model_uv.inverse_transform([0, 0])\n",
        "# Reconstruct spectra from the three outliers in PC2\n",
        "outlier_idx = PCA_model_uv.outlier(X, comps=[1])\n",
        "out_scores = PCA_model_uv.scores[outlier_idx_extreme, :]\n",
        "outliers = PCA_model_uv.inverse_transform(out_scores)\n",
        "# Reconstruct a spectrum for the \"mean\" of these outliers\n",
        "mean_outlier = PCA_model_uv.inverse_transform(out_scores.mean(axis=0))\n",
        "\n",
        "plt.figure()\n",
        "# \"center\" representative \"normal\" sample plotted in blue\n",
        "plt.plot(ppm, model_center_sample, 'b')\n",
        "# The outliers plotted in dashed red line\n",
        "plt.plot(ppm, outliers.T, 'r--',)\n",
        "# The mean outlier plotted in green\n",
        "plt.plot(ppm, mean_outlier, 'g')\n",
        "\n",
        "plt.gca().invert_xaxis()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PZ4yMEfAC2jZ"
      },
      "source": [
        "Another useful measure for outlier detection is the DmodX (Distance to model X) measure. It is calculated based on the model residuals, and can be used to screen for samples which are poorly modelled by the PCA model.\n",
        "A DmodX plot shows the DmodX value for each sample. Higher values mean more unexplained variation (residual) for that particular sample. An exclusion criterion is delineated by using the F-statistic to discern if each sample seems to have a higher amount of unexplained variation than overall in the population of samples modelled.\n",
        "\n",
        "The DmodX outliers are more useful to decide wether or not a sample should is overall well \"explained\" by the model. In this case we have quite a few samples above the critical line, but only 4 which seem to be strongly outlying in DmodX.\n",
        "\n",
        "It is expected that in a variable biological population many observations will contain features not incorporated in a multivariate model fitted on the entire dataset. Adding more principal components could in theory add the residual variability not modelled which is responsible for the DmodX values. In this particular instance we do not take any decision based on the DmodX metric."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "NT8WBop6C2jZ"
      },
      "outputs": [],
      "source": [
        "# The DmodX plot\n",
        "PCA_model_uv.plot_dmodx(X, label_outliers=True)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "pFFb18peC2jZ"
      },
      "outputs": [],
      "source": [
        "# The outlier function can also be used to obtain the DmodX measure and outliers detected with it\n",
        "outlier_idx = PCA_model_uv.outlier(X, measure='DmodX', alpha=0.05)\n",
        "print(outlier_idx)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WfArp__-C2jZ"
      },
      "source": [
        "## Model interpretation\n",
        "\n",
        "We now refit the model without the outliers, and re-assess the findings. The *outlier* method returns the indexes of outlying obersavtions. In the next cell we call it with the *comps* argument = [1], to estimate outliers in the 2nd Principal Component.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "s3lE9hXjC2jZ"
      },
      "outputs": [],
      "source": [
        "outlier_idx = PCA_model_uv.outlier(X, comps=[0]) # indices for all outliers\n",
        "outlier_idx = PCA_model_uv.outlier(X, comps=[1]) # indices for all outliers\n",
        "outlier_idx_extreme = np.array([169, 263, 283])\n",
        "print(\"The following samples (row index) have been detected as outliers in PC1: {0}\".format(outlier_idx))\n",
        "print(\"The following samples (row index) have been detected as outliers PC2: {0}\".format(outlier_idx))\n",
        "print(\"The following samples (row index) have been detected as extreme outliers: {0}\".format(outlier_idx_extreme))\n",
        "#Delete the outlier observations (rows)\n",
        "X_rem = np.delete(X, outlier_idx_extreme, axis=0)\n",
        "Y_rem = np.delete(Y, outlier_idx_extreme, axis=0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fKFPpahTC2ja"
      },
      "source": [
        "After removing outliers, it is recommened to re-assess the model performance using cross validation, and to check whether a model with the same number of components as chosen before is still reliable."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "LUTn7FacC2ja"
      },
      "outputs": [],
      "source": [
        "# Create and fit the PCA model - UV scaling\n",
        "PCA_model_uv = ChemometricsPCA(ncomps=2, scaler=scaling_object_uv)\n",
        "PCA_model_uv.fit(X_rem)\n",
        "PCA_model_uv.scree_plot(X_rem, total_comps=10)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "0fVqjm2tC2ja"
      },
      "outputs": [],
      "source": [
        "rep_cv = PCA_model_uv.repeated_cv(X_rem, repeats=5, total_comps=10)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "1m9rBN_sC2ja"
      },
      "outputs": [],
      "source": [
        "PCA_model_uv.cross_validation(X_rem)\n",
        "print(\"The estimated Q2X from the model is {0}\".format(PCA_model_uv.cvParameters['Q2']))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7AuJMJmtC2ja"
      },
      "source": [
        "Applying similar criteria as before, we now refit a PCA model with a total of 2 components."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "G14QB1QlC2ja"
      },
      "outputs": [],
      "source": [
        "PCA_model_uv.plot_scores(label_outliers=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XZobmTjIC2ja"
      },
      "source": [
        "**Note**: Excluding outliers and re-fitting the model can uncover new candidate outliers. It is not the purpose of this exploratory analysis to obtain a completely *outlier* free dataset for subsequent modelling. We recommend using PCA to screen mainly for large outliers associated with the main Principal components (those who explain a large proportion of the dataset variance) and investigate whether these could be a potential problem in other analyses. If that is the case, further actions might be suggested by the model interpretation, for example, applying some type of batch effect correction or repeating data-preprocessing steps. Although some samples might be outliers due to biological reasons, we do not recommend their automatic exclusion from any further statistical analysis without a strong rationale behind it."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_lBGQLmCC2jb"
      },
      "source": [
        "### Exploring the trends in score plots\n",
        "\n",
        "The *plot_scores* method can use the values of a covariate (discrete or continuous) to colour the scores for each observation. In the next plots we will use the gender information to see if there are any biological trends detected in the first PCA components."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "HGed-YAUC2jb"
      },
      "outputs": [],
      "source": [
        "# Age seems to be one of the main driving forces of variation in the dataset, judging from component 1.\n",
        "PCA_model_uv.plot_scores(color=Y_rem, discrete=True)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "MqWfB2fzC2jb"
      },
      "outputs": [],
      "source": [
        "# The loadings for component number 1\n",
        "ax = PCA_model_uv.plot_loadings(ncomp=1, xaxis=ppm)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "KOkuEu1JaL-3"
      },
      "outputs": [],
      "source": [
        "plotLoadings(PCA_model_uv.loadings[0, :], ppm, spectra=X, cbarlabels = \"Loadings\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "amCDHv_eC2jb"
      },
      "outputs": [],
      "source": [
        "#PCA_model_uv.plot_scores(color=Y_rem, discrete=True, comps=[1, 2]) only run if a model has more than 2 components and plots PC2 vs PC3"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jdqTBLfNC2jb"
      },
      "source": [
        "#### PCA is a very usefull exploratory data analysis tool, especially valuable to visualise the main trends in complex multivariate datasets. It can be very usefull for outlier detection and for preliminary data quality assessement and presence of batch or run-order effects.\n",
        "\n",
        "**Note**: Always investigate as thoroughly as possible why an observation is an outlier, and record all samples that were excluded from an analysis.\n",
        "\n",
        "Although it can also be used to investigate biological differences, supervised methods are more apropriate for that purpose, because they can specifically measure the \"strength\" and effect size of metabolome/phenotype associations.\n",
        "\n",
        "In the next notebook, *Multivariate Analysis - PLS-DA* we will use a supervised model to explicitly investigate metabolic profile differences according to gender, and discuss PLS-DA model interpretation."
      ]
    }
  ],
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python [conda env:npcproject]",
      "language": "python",
      "name": "conda-env-npcproject-py"
    },
    "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.18"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}