Multivariate Analysis - PCA.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "1E02NPscV97D"
   },
   "source": [
    "# Multivariate Analysis - PCA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "jzQxU7A9V97F"
   },
   "source": [
    "In this notebook we will perform an unsupervised *multivariate* exploratory analysis with principal component analysis (PCA) of the [AddNeuroMed](https://nyaspubs.onlinelibrary.wiley.com/doi/10.1111/j.1749-6632.2009.05064.x) dataset of Alzheimer's disease.\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",
    "4) **Model Interpretation**\n",
    "\n",
    "\n",
    "<br>\n",
    "\n",
    "*Credits: This tutorial was originally created by Gonçalo Correia and was adapted by Lukas Kopecky and Frederico Soares in November/December 2023.*"
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Preliminary Steps"
   ],
   "metadata": {
    "collapsed": false,
    "id": "v9IeUI_hV97F"
   }
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "5-upppX_V97H"
   },
   "source": [
    "### Import Packages\n",
    "\n",
    "Import all the packages and configure notebook plotting mode."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "sRoSV1tbV97I",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701685675791,
     "user_tz": 0,
     "elapsed": 3425,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    }
   },
   "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 pd\n",
    "import matplotlib.pyplot as plt\n",
    "import warnings\n",
    "from sklearn.exceptions import DataConversionWarning\n",
    "\n",
    "from pyChemometrics.ChemometricsPCA import ChemometricsPCA\n",
    "from pyChemometrics.ChemometricsScaler import ChemometricsScaler\n",
    "from pyChemometrics.plotting_utils import plotLoadings\n",
    "from pyChemometrics.plotting_utils import _scatterplots\n",
    "\n",
    "\n",
    "# Use to obtain same values as in the text\n",
    "np.random.seed(350)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Set the data conversion warnings to appear only once to avoid repetition during CV\n",
    "warnings.filterwarnings(\"ignore\", category=DataConversionWarning)\n",
    "warnings.filterwarnings(\"ignore\", category=FutureWarning)"
   ],
   "metadata": {
    "id": "7D2Acm92V97J",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701685675792,
     "user_tz": 0,
     "elapsed": 12,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    }
   }
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "zcjVOOGAV97J"
   },
   "source": [
    "The next cell sets up the figure display mode. The *notebook* mode allows interactive plotting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "-uMQJw7uV97K",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701685675792,
     "user_tz": 0,
     "elapsed": 9,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    }
   },
   "outputs": [],
   "source": [
    "# Set the plot backend to support interactive plotting\n",
    "# %matplotlib notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "efArnUJaV97K"
   },
   "source": [
    "### Data import"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "b57aOsrtV97K"
   },
   "source": [
    "We will now import the LC-MS data with the metadata (Y variables) and feature annotation for LC-MS.\n",
    "\n",
    "Then we split the data into two parts:\n",
    "\n",
    "**X** - LC-MS data matrix\n",
    "\n",
    "**Y** - Metadata that will help with the interpretation of the PCA model (sex/gender in this instance)\n",
    "\n",
    "**Y2** - Metadata that will help with the interpretation of the PCA model (age in this instance)\n",
    "\n",
    "Then we extract the feature annotations:\n",
    "\n",
    "**retention_times**, **mz_values** - annotation for the features of the rpos_x_matrix data\n",
    "\n",
    "<br>\n",
    "\n",
    "*NB - The dataset used in this tutorial has reduced features to speed up the model fitting proces. Full data available from [https://zenodo.org/doi/10.5281/zenodo.4053166](https://zenodo.org/doi/10.5281/zenodo.4053166). We recommend you to try to run the tutorial with the full dataset in your own time, or you can even try to run it using your own data.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Load the dataset\n",
    "dementia_rpos_dataset = pd.read_csv(\"./data/Dementia U RPOS_combinedData.csv\",delimiter=',')\n",
    "\n",
    "# Inspect the dataset\n",
    "dementia_rpos_dataset.head(10)"
   ],
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 603
    },
    "id": "VyvLyVI5V97K",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701685680025,
     "user_tz": 0,
     "elapsed": 4241,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    },
    "outputId": "3a4cd249-8a67-4a43-9f88-97dabff0554e"
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "#### Data Filtering"
   ],
   "metadata": {
    "collapsed": false,
    "id": "lKoFrJJuV97L"
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "zIDOMp_HV97L",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701685680026,
     "user_tz": 0,
     "elapsed": 6,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "# Delete samples where outcome variable is unknown (In this example QC samples)\n",
    "dementia_rpos_dataset = dementia_rpos_dataset[~dementia_rpos_dataset['Gender'].isnull()]\n",
    "\n",
    "X = dementia_rpos_dataset.iloc[:, 5::].values\n",
    "\n",
    "variable_names = dementia_rpos_dataset.columns[5::]\n",
    "\n",
    "# Use pandas Categorical type\n",
    "Y = pd.Categorical(dementia_rpos_dataset['Gender']).codes\n",
    "Y2 = dementia_rpos_dataset['Age'].values\n",
    "\n",
    "# Extract the retention times and m/z to use in 2D plots of the dataset\n",
    "retention_times = np.array([x.split('_')[0] for x in variable_names], dtype='float')/60\n",
    "mz_values = np.array([x.split('_')[1][0:-3] for x in variable_names], dtype='float')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "qhPeekHyV97L"
   },
   "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 = pd.read_csv(\"./data/X_mass_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",
    "   Y1 = pd.Categorical(Y.iloc[:, 1])\n",
    "   Y1.codes # The numerical label\n",
    "   Y1.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": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "print(dementia_rpos_dataset['Gender'].value_counts())\n",
    "print(pd.DataFrame(Y).value_counts())"
   ],
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "3btDOqH0V97M",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701685680026,
     "user_tz": 0,
     "elapsed": 6,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    },
    "outputId": "a4b2de86-7913-4d6b-fbe6-a32077a357ea"
   }
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "NgkNBegNV97M"
   },
   "source": [
    "Plot the log intensities of the features in the dataset. Each point in the scatterplot represents a feature (m/z and retention time combination) and the color represents the log intensity of that feature in the sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 449
    },
    "id": "JOuy9tyIV97M",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701685686836,
     "user_tz": 0,
     "elapsed": 6813,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    },
    "outputId": "4ceff976-75ff-41db-b78e-fe21ae5dfdf2"
   },
   "outputs": [],
   "source": [
    "_scatterplots(np.log1p(X).mean(axis=0), xaxis=retention_times, yaxis=mz_values, marker_size=3, alpha=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "wlUxwI3mV97N"
   },
   "source": [
    "## 1) PCA model fitting\n",
    "\n",
    "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": "sfseQxVJV97N"
   },
   "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": "AQ9sRYDNV97N",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701685687290,
     "user_tz": 0,
     "elapsed": 465,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    }
   },
   "outputs": [],
   "source": [
    "# Start with log scaling of the data\n",
    "X = np.log1p(X)\n",
    "\n",
    "# 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": "1tJnKCWHV97N"
   },
   "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": "JtAXVTbQV97O"
   },
   "source": [
    "Pass each scaling object to the PCA method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "rRDs2AMiV97O",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701685690096,
     "user_tz": 0,
     "elapsed": 2821,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    }
   },
   "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": "p63tA_TaV97O",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701685693353,
     "user_tz": 0,
     "elapsed": 3278,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    }
   },
   "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": "-uiDFHTCV97O",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701685695166,
     "user_tz": 0,
     "elapsed": 1837,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    }
   },
   "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": "v-DxJa5WV97O"
   },
   "source": [
    "### Effect of scaling on PCA Score plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Zy7gjs6cV97O"
   },
   "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": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 516
    },
    "id": "KtdbsLL8V97P",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701685695881,
     "user_tz": 0,
     "elapsed": 721,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    },
    "outputId": "394075cd-9c70-41b2-daec-54a49c6192d3"
   },
   "outputs": [],
   "source": [
    "# PCA score plot for the mean centered model\n",
    "PCA_model_mc.plot_scores(comps=[0, 1], plot_title='Mean centering')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 516
    },
    "id": "PnI9sT5nV97P",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701685696446,
     "user_tz": 0,
     "elapsed": 582,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    },
    "outputId": "5b655ac7-e18e-467b-a6a7-7865fa4084ae"
   },
   "outputs": [],
   "source": [
    "# Score plot for the Pareto scaled model\n",
    "PCA_model_par.plot_scores(comps=[0, 1], plot_title='Pareto scaling')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 516
    },
    "id": "s4bG3nK-V97P",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701685697032,
     "user_tz": 0,
     "elapsed": 600,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    },
    "outputId": "152138c5-620e-4472-da8f-e80035c95e76"
   },
   "outputs": [],
   "source": [
    "# PCA score plot for UV scaled model\n",
    "PCA_model_uv.plot_scores(comps=[0, 1], plot_title='UV scaling')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "MFevEYV_V97Q"
   },
   "source": [
    "### Effect of scaling on PCA loadings"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "FZq1fs4KV97Q"
   },
   "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": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 455
    },
    "id": "NwwQEttFV97Q",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701685701360,
     "user_tz": 0,
     "elapsed": 4343,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    },
    "outputId": "fc2540cc-5271-4386-9aea-48dd20e1c6f6"
   },
   "outputs": [],
   "source": [
    "# Plot of first principal component loadings of mean centering model\n",
    "PCA_model_mc.plot_loadings(ncomp=1,xaxis=retention_times, yaxis=mz_values,instrument='lcms', marker_size=5, alpha = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 452
    },
    "id": "QlVUWsRmV97Q",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701685705239,
     "user_tz": 0,
     "elapsed": 3901,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    },
    "outputId": "b86a4aa1-5eb3-4dff-d407-8429fef1e504"
   },
   "outputs": [],
   "source": [
    "# Plot of first principal component loadings of Pareto scaled model\n",
    "PCA_model_par.plot_loadings(ncomp=1,xaxis=retention_times, yaxis=mz_values,instrument='lcms', marker_size=5, alpha = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 449
    },
    "id": "3IfDB_k6V97Q",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701685708661,
     "user_tz": 0,
     "elapsed": 3431,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    },
    "outputId": "ab27d9d9-b9c2-46c1-d94d-c1c4ee591bd2"
   },
   "outputs": [],
   "source": [
    "# Plot of first principal component loadings of Unit Variance scaled model\n",
    "PCA_model_uv.plot_loadings(ncomp=1,xaxis=retention_times, yaxis=mz_values,instrument='lcms', marker_size=5, alpha = 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "EuF2iTJCV97R"
   },
   "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": "tlTb1VfgV97R"
   },
   "source": [
    "## 2) Model cross-validation and component selection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "MXIgOPsVV97S"
   },
   "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": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 468
    },
    "id": "g434iHCVV97S",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701685802394,
     "user_tz": 0,
     "elapsed": 93749,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    },
    "outputId": "f6370d93-f0f7-4204-84ad-b9a7851aad43"
   },
   "outputs": [],
   "source": [
    "PCA_model_uv.scree_plot(X, total_comps=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "8FqR0_18V97S"
   },
   "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": "uUvVweX4V97S"
   },
   "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, it's still important to benchmark the cross-validated performance of the model, to avoid interpreting or infering biological conclusions from non-robust principal components. **DO NOT RUN THIS PART IN GOOGLE COLABORATORY AS THIS WOULD TAKE UP TO 20 MINUTES TO RUN**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "v3W1A4tLV97S",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701685802397,
     "user_tz": 0,
     "elapsed": 31,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    }
   },
   "outputs": [],
   "source": [
    "# # DO NOT RUN IN GOOGLE COLAB\n",
    "# rep_cv = PCA_model_uv.repeated_cv(X, repeats=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "s-UYMofaV97S"
   },
   "source": [
    "Refit the model for further exploration with the selected number of components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "GZc9hawsV97T",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701685811448,
     "user_tz": 0,
     "elapsed": 9068,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    }
   },
   "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": "lCeQVNu_V97T"
   },
   "source": [
    "## 3) Outlier detection and model interpretation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "BxlLemn1V97T"
   },
   "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": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 516
    },
    "id": "VwstjBVQV97T",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701685811956,
     "user_tz": 0,
     "elapsed": 526,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    },
    "outputId": "c7e5d98f-1773-433f-b83a-0821b52e7938"
   },
   "outputs": [],
   "source": [
    "PCA_model_uv.plot_scores(comps=[0, 1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "LkOlhBv1V97U"
   },
   "source": [
    "The PCA score plot highlighted a set of candidate outliers. There is a cluster of 14 main 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": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "XuDpOzMaV97U",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701685811957,
     "user_tz": 0,
     "elapsed": 21,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    },
    "outputId": "0d870970-f6a8-424d-ed33-685735637df5"
   },
   "outputs": [],
   "source": [
    "outlier_idx = PCA_model_uv.outlier(X)\n",
    "print(\"Outliers for the full 4 component model : {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": "-6dj2CNiV97U"
   },
   "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 491
    },
    "id": "xmRRT6_YV97U",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701685819142,
     "user_tz": 0,
     "elapsed": 7198,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    },
    "outputId": "ccdb238e-9b86-4a64-8b2e-ff0ec52d8798"
   },
   "outputs": [],
   "source": [
    "# plot the mean spectrum calculated from the raw data (blue) and the outlying intensity data (red)\n",
    "plt.figure()\n",
    "plt.plot( X[outlier_idx, :].T, 'r')\n",
    "plt.plot( np.mean(X, axis=0), 'b')\n",
    "plt.gca()\n",
    "plt.title('Comparison of mean intensities with intensities of outliers')\n",
    "plt.ylabel('Intensity')\n",
    "plt.xlabel('Feature Index')\n",
    "\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "i9Mp6CdIV97U"
   },
   "source": [
    "The PCA score plot should be interpreted using the loading vectors for the corresponding components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 449
    },
    "id": "rM72m3V2V97U",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701685823531,
     "user_tz": 0,
     "elapsed": 4402,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    },
    "outputId": "a6da6b60-3a39-4a66-d682-43018d9a9f56"
   },
   "outputs": [],
   "source": [
    "PCA_model_uv.plot_loadings(ncomp=1,xaxis=retention_times, yaxis=mz_values,instrument='lcms', alpha=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "TAt_FJhGV97V"
   },
   "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 peak intensities. We can explore the PCA model space in this way, and compare peak intensities representative of any desired region of the score plot.\n",
    "\n",
    "The following plots show the main features, and if an outlier exhibits a value markedly higher than the raw data (depicted in grey), it is visualized in accordance with the colour-bar."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Xl-MYPxfV97V",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 449
    },
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701685824445,
     "user_tz": 0,
     "elapsed": 936,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    },
    "outputId": "45977934-8955-41de-aac0-7ac644f224b0"
   },
   "outputs": [],
   "source": [
    "PCA_model_uv.plot_outliers(X, outlier_idx, instrument='lcms', xaxis=retention_times, yaxis=mz_values, sigma = 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "L8oLc53nV97V"
   },
   "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-Statistics to discern if each sample seems to have a higher amount of unexplained variation than overall in the population of samples modelled. F-statistic which assumes Normal distribution of residuals. We also expect 5% of samples beyond the 95% CI.\n",
    "\n",
    "The DmodX outliers are more useful to decide whether a sample is overall well \"explained\" by the model. In this case we have quite a few samples above the critical line, but only 19 which seem to be strongly outlying in both DmodX and Hotelling T2.\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": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 468
    },
    "id": "yhIkOs0FV97V",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701685825350,
     "user_tz": 0,
     "elapsed": 913,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    },
    "outputId": "78e24ef5-146f-494d-b4e1-acaa7a6ce841"
   },
   "outputs": [],
   "source": [
    "# The DmodX plot\n",
    "PCA_model_uv.plot_dmodx(X, label_outliers=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "R5F-dS1cV97V",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701685825350,
     "user_tz": 0,
     "elapsed": 21,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    },
    "outputId": "52c27e50-1d51-48c1-e9e6-a5933325a6e5"
   },
   "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": "hHCQeaotV97V"
   },
   "source": [
    "## 4) 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 = [0], to estimate outliers in the 1st Principal Component.\n",
    "We seleceted the 1st principal component based on the last scores plot, where the outliers are in the positive value of PC1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "mYGX4k0fV97V",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701685825807,
     "user_tz": 0,
     "elapsed": 465,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    },
    "outputId": "9c1c4d07-c806-4770-ae3e-f76a7b31dc5c"
   },
   "outputs": [],
   "source": [
    "outlier_idx = PCA_model_uv.outlier(X, comps=[0])\n",
    "\n",
    "print(\"The following samples (row index) have been detected as outliers: {0}\".format(outlier_idx))\n",
    "#Delete the outlier observations (rows)\n",
    "X_rem = np.delete(X, outlier_idx, axis=0)\n",
    "Y1_rem = np.delete(Y, outlier_idx, axis=0)\n",
    "Y2_rem = np.delete(Y2, outlier_idx, axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "upGfoSz9V97V"
   },
   "source": [
    "After removing outliers, it is recommended 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": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 468
    },
    "id": "2Zg3Jl_fV97W",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701685915880,
     "user_tz": 0,
     "elapsed": 90089,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    },
    "outputId": "60281381-4778-44eb-8c56-bf0786b7abf9"
   },
   "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": "_mrXgmYoV97W",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701685915881,
     "user_tz": 0,
     "elapsed": 24,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    }
   },
   "outputs": [],
   "source": [
    "# # DO NOT RUN IN GOOGLE COLABORATORY\n",
    "# rep_cv = PCA_model_uv.repeated_cv(X_rem, repeats=5, total_comps=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "a5FLXyeLV97W",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701685923605,
     "user_tz": 0,
     "elapsed": 7738,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    },
    "outputId": "6ddd4188-0d1f-4fc3-8089-db2009ba6658"
   },
   "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": "z3k-He6YV97W"
   },
   "source": [
    "Applying similar criteria as before, we now refit a PCA model with a total of 4 components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "PCA_model_uv = ChemometricsPCA(ncomps=4, scaler=scaling_object_uv)\n",
    "PCA_model_uv.fit(X_rem)"
   ],
   "metadata": {
    "id": "rVzKZBJyV97W",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701686425670,
     "user_tz": 0,
     "elapsed": 2031,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 516
    },
    "id": "L_sy3hL8V97W",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701686427144,
     "user_tz": 0,
     "elapsed": 10,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    },
    "outputId": "37f5d0f7-a525-4cca-edc6-aa7e1260ad10"
   },
   "outputs": [],
   "source": [
    "PCA_model_uv.plot_scores()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ArBaLMaWV97W"
   },
   "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": "NrkkfRfwV97W"
   },
   "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 Age and Genotype information to see if there are any biological trends detected in the first PCA components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 516
    },
    "id": "nfTHvolnV97W",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701686433331,
     "user_tz": 0,
     "elapsed": 1238,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    },
    "outputId": "8971e32f-3623-4216-d8fa-59f1f93ababd"
   },
   "outputs": [],
   "source": [
    "# Age does not seem to be  main driving forces of variation in the dataset, judging from component 2 and 3.\n",
    "PCA_model_uv.plot_scores(color=Y2_rem, discrete=False, comps=[2, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 516
    },
    "id": "x86jO4vzV97W",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701686434401,
     "user_tz": 0,
     "elapsed": 645,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    },
    "outputId": "2854ffb2-6fdf-43aa-ebfe-4e48d525bff8"
   },
   "outputs": [],
   "source": [
    "# Age does not seem to have any impact on the biological variation either.\n",
    "PCA_model_uv.plot_scores(color=Y1_rem, discrete=True, comps=[0, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "scores_df = pd.DataFrame(PCA_model_uv.scores, columns=range(1, PCA_model_uv.scores.shape[1]+1))\n",
    "scores_df['Gender'] = Y1_rem\n",
    "sns.pairplot(\n",
    "    data=scores_df,\n",
    "    hue='Gender'\n",
    ")\n",
    "\n",
    "plt.xlabel('Component')\n",
    "plt.ylabel('Component')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ],
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "4caT0-e4V97X",
    "executionInfo": {
     "status": "ok",
     "timestamp": 1701686444621,
     "user_tz": 0,
     "elapsed": 8025,
     "user": {
      "displayName": "Lukas Kopecky",
      "userId": "18051270114113256328"
     }
    },
    "outputId": "a50af1bf-fd00-448c-a8bd-c1aef9a31457"
   }
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "uZCv0A1lV97X"
   },
   "source": [
    "**PCA is a very useful exploratory data analysis tool, especially valuable to visualise the main trends in complex multivariate datasets. It can be very useful for outlier detection and for preliminary data quality assessment 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",
    "\n",
    "\n",
    "Although it can also be used to investigate biological differences or differences caused by analytical errors, supervised methods are more appropriate 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 sex, and discuss PLS-DA model interpretation."
   ]
  }
 ],
 "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.11.5"
  },
  "colab": {
   "provenance": []
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}