week09a-PCA-EOF.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Principal Component Analysis (PCA) and Empirical Orthogonal Functions (EOFs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The goal of *principal component analysis (PCA)* is to efficiently represent variability that is common among multiple variables. One application is to create an \"index\" that describes how different variables co-vary, or how a single variable co-varies in different places. *Empirical orthogonal functions (EOFs)* are a very similar type of analysis that involves the same basic steps (often seen in Physical Oceanography or Meteorology literature).\n",
    "\n",
    "## Visualizing in two dimensions (two variables)\n",
    "\n",
    "PCA is easiest to visualize when there are only two variables. The same data can be visualized on different axes if you simply rotate the data into a new coordinate system with new axes. \n",
    "\n",
    "![images/PCA_1.png](images/PCA_1.png)\n",
    "\n",
    "*Image source*: Joliffe (1990) Weather\n",
    "\n",
    "In the figure above, the data are rotated into a new coordinate system that maximizes the variance, or spread, in the data along the x-axis (the major axis) and minimizes the variance in the y-axis (the minor axis). This particular type of PCA in two dimensions is often referred to as __principal axis__ analysis. This analysis, which is described in detail below, creates a new coordinate frame that maximizes variability along an axis (PC1) and corresponding orthogonal axis (PC2) where PC1 PC2 are uncorrelated.\n",
    "    \n",
    "A common application of this type of analysis is to rotate data from a East/North coordinate system to *natural coordinates* primary oriented with the coastline or topography.\n",
    "\n",
    "#### Example - Wind vectors\n",
    "\n",
    "![images/PCA_2.png](images/PCA_2.png)\n",
    "\n",
    "Source: Patrick Daniel, MS 263, Spring 2016\n",
    "\n",
    "#### Example from Emery and Thompson - Current meter data\n",
    "\n",
    "![images/PCA_3.png](images/PCA_3.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Principal component analysis: steps\n",
    "1. Create data matrix (size: N x M, N rows and M columns)\n",
    "   - standardize if necessary (if data is a combination of different variables)\n",
    "2. Form the covariance matrix (M x M)\n",
    "3. Extract eigenvalues and eigenvectors from covariance matrix\n",
    "4. The eigenvectors are the principal components, the eigenvalues are their magnitudes\n",
    "\n",
    "Here these steps are decribed for a two dimensional example, focusing on the current meter example above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1. Data matrix, __X__\n",
    "\n",
    "Creating a data matrix, $X$\n",
    "\n",
    "N rows, M columns\n",
    "\n",
    "$$\\begin{bmatrix} u_1 & v_1 \n",
    "\\\\ u_2 & v_2 \n",
    "\\\\ u_3 & v_3\n",
    "\\\\ \\vdots & \\vdots\n",
    "\\\\ u_N & v_N \\end{bmatrix}$$\n",
    "                \n",
    "N - number of concurrent samples \n",
    "\n",
    "M - number of variables\n",
    "\n",
    "For the current meter example, $u$ is eastward velocity and $v$ is northward velocity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2. Covariance matrix, C  (or correlation matrix, R)\n",
    "\n",
    "The covariance matrix __C__ has  M rows and M columns (M x M). It contains the covariance between all possible combinations of the M different variables that are stored in the different columns of __X__.\n",
    "    \n",
    "$$\\begin{bmatrix}\n",
    "c_{11} & c_{12} & c_{13} & \\cdots & c_{1M}\\\\\n",
    "c_{21} & c_{22} & c_{23} & \\cdots & c_{2M}\\\\\n",
    "c_{31} & c_{32} & c_{33} & \\cdots & c_{3M}\\\\\n",
    "\\vdots & \\vdots & \\vdots & \\ddots &       \\\\\n",
    "c_{M1} & c_{M2} & c_{M3} &        & c_{MM}\\\\\n",
    "\\end{bmatrix}$$\n",
    "        \n",
    "Alternatively, the correlation matrix, __R__, can also be used. Recall that the correlation between two variables is the covariance scaled by the standard deviation of the variables, $r_{uv} = \\frac{c_{uv}}{s_u s_v}$. \n",
    "   \n",
    "$$\\begin{bmatrix}\n",
    "1      & r_{12} & r_{13} & \\cdots & r_{1M}\\\\\n",
    "r_{21} & 1      & r_{23} & \\cdots & r_{2M}\\\\\n",
    "r_{31} & r_{32} & 1      & \\cdots & r_{3M}\\\\\n",
    "\\vdots & \\vdots & \\vdots & \\ddots &       \\\\\n",
    "r_{M1} & r_{M2} & r_{M3} &        & 1     \\\\\n",
    "\\end{bmatrix}$$\n",
    "\n",
    "The diagonal elements of the correlation matrix are all equal to 1. \n",
    "\n",
    "The covariance matrix is generally used when all variables have the same units. Variables or locations with greater variance are weighted higher. The correlation matrix is used when the variables have different units, or if the overall magnitudes do not matter. Use of the correlation matrix standardizes the data so that they are all weighted equally. It is equivalent to transforming the data into z-scores first, and then computing the covariance matrix.\n",
    "   \n",
    "##### Current meter example: covariance matrix \n",
    "\n",
    "For the case of finding the principal component of wind or current meter data, the variables are the eastward component ($u$) and the northward component ($v$). The covariance matrix has two rows and two columns.\n",
    "\n",
    "$$\\begin{bmatrix}\n",
    "c_{uu} & c_{uv}\\\\\n",
    "c_{vu} & c_{vv}\\\\\n",
    "\\end{bmatrix}$$\n",
    "\n",
    "$c_{uu}$  variance of $u$ <br>\n",
    "\n",
    "$c_{vv}$  variance of $v$<br>\n",
    "\n",
    "$c_{uv}$  covariance of $u$ and $v$ <br>\n",
    "\n",
    "The values for this example are given in Emery and Thomson,\n",
    "\n",
    "$$\\begin{bmatrix}\n",
    "441 & -302.4\\\\\n",
    "-302.4 & 324\\\\\n",
    "\\end{bmatrix}$$\n",
    "\n",
    "where units are in $\\frac{cm^2}{s^2}$\n",
    "\n",
    "Calculating the covariance matrix can also be done in Python with the numpy.cov() function. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3. Extracting the eigenvalues and eigenvectors from the covariance matrix\n",
    "\n",
    "The eigenvectors (contained in the matrix $V$) and eigenvalues (contained in the matrix $\\Lambda$) of the $C$ matrix are defined by the relationship:\n",
    "\n",
    "$$CV =  V \\Lambda$$\n",
    "\n",
    "$V$ =  eigenvectors (M x M matrix)\n",
    "\n",
    "$$\\begin{bmatrix}\n",
    "v_{11} & v_{12} & v_{13} & \\cdots\\\\\n",
    "v_{21} & v_{22} & v_{23} & \\cdots\\\\\n",
    "v_{31} & v_{32} & v_{33} & \\cdots\\\\\n",
    "\\vdots & \\vdots & \\vdots & \\ddots\n",
    "\\end{bmatrix}$$\n",
    "\n",
    "Each column is one eigenvector\n",
    "\n",
    "$\\Lambda$ = eigenvalues (M x M)\n",
    "\n",
    "$$\\begin{bmatrix}\n",
    "\\lambda_{1} & 0 & 0 & \\cdots \\\\\n",
    "0 & \\lambda_{2} & 0 & \\cdots \\\\\n",
    "0 & 0 & \\lambda_{3} & \\cdots \\\\\n",
    "\\vdots & \\vdots & \\vdots & \\ddots\n",
    "\\end{bmatrix}$$\n",
    "\n",
    "where $\\lambda _1$ - eigenvalue # 1, $\\lambda _2$ - eigenvalue # 2, etc.\n",
    "\n",
    "#### 4. Physical meaning of the eigenvalues and eigenvectors\n",
    "\n",
    "The *eigenvectors* describe directions in M-dimensional space. These are called the principal components, which define the orientation of new axes for the data. All eigenvectors are orthogonal (or perpendicular, for the 2D case). The vectors that define the axes can be scaled any magnitude, but usually they are scaled to have a value of one (unit vectors).\n",
    "\n",
    "The *eigenvalues* describe how each eigenvector is scaled. They describe the variance in the data along each axis defined by the orientations of the eigenvectors.\n",
    "\n",
    "#### Back to the current meter example\n",
    "\n",
    "To compute the eigenvectors and eigenvalues of a matrix in Python, use:\n",
    "```python\n",
    "from scipy import linalg\n",
    "linalg.eig(C)\n",
    "```\n",
    "<br>\n",
    "\n",
    "For the covariance matrix $C$ given above,\n",
    "\n",
    "$$V = \\begin{bmatrix}\n",
    ".771 & 0.636\\\\\n",
    "-.636 & .771\n",
    "\\end{bmatrix}$$\n",
    "\n",
    "The columns of the V matrix define the orientation of the new axes (or principal components) in $u,v$ space. The first column defining PC 1 has a negative slope ($u$ = 0.771, $v$ = -0.636) because the variance in the data is primarily oriented in the northwest/southeast axis in $u,v$ space.\n",
    "\n",
    "$$\\Lambda = \\begin{bmatrix}\n",
    "690.5 & 0\\\\\n",
    "0 & 74.5\n",
    "\\end{bmatrix}$$\n",
    "\n",
    "\n",
    "$\\lambda _1$ = 690.6 $\\frac{cm^2}{s^2}$ = variance along PC 1 <br>\n",
    "$\\lambda _2$ = 74.5 $\\frac{cm^2}{s^2}$ = variance along PC 2  <br>\n",
    "\n",
    "Note that there is much more variance along PC 1 (the major axis) than PC 2 (the minor axis). The data have been rotated so that the variance is maximized along PC 1 and minimized along PC 2.\n",
    "\n",
    "To recover the data in the rotated coordinates,\n",
    "\n",
    "$\\tau _1 = X V_i$\n",
    "\n",
    "$\\tau _1$ = u(0.771) + v(-0.636) <br>\n",
    "$\\tau _2$ = u(0.636) + v(0.771)\n",
    "\n",
    "For this example, this would create new variables for alongshore and cross-shore velocity. This is the same data, just rotated into a new set of axes, where alongshore is defined by the axis of maximum variance (the major axis).\n",
    "\n",
    "Factor loading is another common terminology encountered in EOFs, and refers to the larger spread of data in the first eigenvector (original example v') than in the second (u')and thereby quantifying their relative importance and are defined by:\n",
    "\n",
    "$A =  V \\sqrt{\\Lambda}$\n",
    "\n",
    "Factor loading can be visualized as the length of a PC vectors, whereby PC1 would have a longer vector compared to PC2 (see Emery and Thompson PCA figure above). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 3-D visualization of principal components\n",
    "\n",
    "For three variables, defined by three axes, three principal components are obtained. The first two define a plane in three-dimensional space. These are the two axes that define a a plane where the data would be most spread out if you looked at the data from a perspective perpendicular to the plane.\n",
    "\n",
    "<img src=\"images/PCA_4.png\" width=\"700\">\n",
    "[source](http://ordination.okstate.edu/PCA.htm)\n",
    "\n",
    "For an interactive visualization of PCA in three dimensions, see this link:\n",
    "\n",
    "http://setosa.io/ev/principal-component-analysis/\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: the Pacific Decadal Oscillation\n",
    "\n",
    "![images/PCA_PDO.png](images/PCA_PDO.png)\n",
    "\n",
    "[source](http://research.jisao.washington.edu/pdo/)\n",
    "\n",
    "The Pacific Decadal Oscillation (PDO) is defined as the first principal component of SST in the North Pacific. The origin of this index comes from Davis (1976).\n",
    "\n",
    "This study used binned SST data from 1947-1974 in 82 different sites across the Pacific Ocean.\n",
    "\n",
    "![images/PCA_PDO_Bins.png](images/PCA_PDO_Bins.png)\n",
    "\n",
    "Davis, R.E. (1976) Predictability of sea surface temperature and sea level pressure anomalies over the North Pacific Ocean. _J. Phys. Oceanogr._, 6, 249-266.\n",
    "\n",
    "![images/PCA_PDO_STD.png](images/PCA_PDO_STD.png)\n",
    "\n",
    "Standard deviation of temperature anomaly [deg C] \n",
    "\n",
    "Temperature anomaly = temperature – seasonal monthly mean temperature\n",
    "\n",
    "#### Data matrix:\n",
    "\n",
    "(336 X 82) -> 28 years * 12 month\n",
    "    \n",
    "#### Create a covariance matrix (what size?) \n",
    "\n",
    "(82 x 82)\n",
    "    \n",
    "#### First principal component/empirical orthogonal function of North Pacific SST\n",
    "\n",
    "![images/PCA_PDO_First.png](images/PCA_PDO_First.png)\n",
    "\n",
    "This is the dominant \"mode\" of temperature variability in the North Pacific, based on how temperature at different locations co-vary together (i.e. the covariance matrix). \n",
    "\n",
    "This PC is a vector defined by 82 different numbers. It is the first column of an (82 x 82) eigenvector matrix (V). It can be visualized by mapping each value back to its corresponding spatial location.\n",
    "\n",
    "#### Other principal components of North Pacific SST\n",
    "\n",
    "![images/PCA_PDO_All.png](images/PCA_PDO_All.png)\n",
    "\n",
    "#### Fraction of Variance that is accounted for by the fist M empirical orthogonal functions\n",
    "\n",
    "![images/PCA_PDO_Var.png](images/PCA_PDO_Var.png)\n",
    "\n",
    "#### Principal component scores\n",
    "\n",
    "The first PC can be used to create an index describing the time variability of the dominant mode. This is the PDO index, which decribes a mode of temperature variability that can be compared with ecosystem parameters (e.g. zooplankton abundance).\n",
    "\n",
    "![images/PCA_PDO_Index.png](images/PCA_PDO_Index.png)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### PCA assumptions\n",
    "- Relationships between variables are linear\n",
    "- Mean and variance are sufficient statistics (variables have Gaussian/normal distribution)\n",
    "- Large variances have import dynamics (high signal to noise ration)\n",
    "- The principal components are orthogonal\n",
    "\n",
    "##### PCA Issues\n",
    "\n",
    "- Does not take into account phase lags between correlated variables, can detect \"standing\" features but not propagating waves. The PDO oscillates back and forth between cool and warm phases like a see-saw, it does not propagate across the Pacific like a wave. If there are wave-like features in the data, PCA will split up the variance between multiple modes, confusing interpration.\n",
    "\n",
    "- Domain dependence: can get a different answer by changing the spatial extent of the data\n",
    "\n",
    "Shlens, J., A tutorial on principal component analysis: Derivation, discussion and singular value decomposition.  https://www.cs.princeton.edu/picasso/mats/PCA-Tutorial-Intuition_jp.pdf \n",
    "\n",
    "Example:\n",
    "![images/PCA_Issue.png](images/PCA_Issue.png)\n",
    "\n",
    "If no dominant pattern exists, PCA will pick out a pattern that may not be robust."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}