Lab 18 - PCA in Python.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This lab on Principal Components Analysis is a python adaptation of p. 401-404,\n",
    "408-410 of \"Introduction to Statistical Learning with Applications in R\" by Gareth James,\n",
    "Daniela Witten, Trevor Hastie and Robert Tibshirani. Original adaptation by J. Warmenhoven, updated by R. Jordan Crouser at Smith College for\n",
    "SDS293: Machine Learning (Spring 2016)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 10.4: Principal Components Analysis\n",
    "\n",
    "In this lab, we perform PCA on the ${\\tt USArrests}$ data set. The rows of the data set contain the 50 states, in\n",
    "alphabetical order:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df = pd.read_csv('USArrests.csv', index_col=0)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The columns of the data set contain four variables relating to various crimes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start by taking a quick look at the column means of the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see right away the the data have **vastly** different means. We can also examine the variances of the four variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df.var()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not surprisingly, the variables also have vastly different variances: the\n",
    "${\\tt UrbanPop}$ variable measures the percentage of the population in each state\n",
    "living in an urban area, which is not a comparable number to the number\n",
    "of crimes committeed in each state per 100,000 individuals. If we failed to scale the\n",
    "variables before performing PCA, then most of the principal components\n",
    "that we observed would be driven by the ${\\tt Assault}$ variable, since it has by\n",
    "far the largest mean and variance. \n",
    "\n",
    "Thus, it is important to standardize the\n",
    "variables to have mean zero and standard deviation 1 before performing\n",
    "PCA. We can do this using the ${\\tt scale()}$ function from ${\\tt sklearn}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import scale\n",
    "X = pd.DataFrame(scale(df), index=df.index, columns=df.columns)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll use the ${\\tt PCA()}$ function from ${\\tt sklearn}$ to compute the loading vectors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sklearn.decomposition import PCA\n",
    "\n",
    "pca_loadings = pd.DataFrame(PCA().fit(X).components_.T, index=df.columns, columns=['V1', 'V2', 'V3', 'V4'])\n",
    "pca_loadings"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that there are four distinct principal components. This is to be\n",
    "expected because there are in general ${\\tt min(n − 1, p)}$ informative principal\n",
    "components in a data set with $n$ observations and $p$ variables.\n",
    "\n",
    "Using the ${\\tt fit_transform()}$ function, we can get the principal component scores of the original data. We'll take a look at the first few states:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Fit the PCA model and transform X to get the principal components\n",
    "pca = PCA()\n",
    "df_plot = pd.DataFrame(pca.fit_transform(X), columns=['PC1', 'PC2', 'PC3', 'PC4'], index=X.index)\n",
    "df_plot.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can construct a **biplot** of the first two principal components using our loading vectors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig , ax1 = plt.subplots(figsize=(9,7))\n",
    "\n",
    "ax1.set_xlim(-3.5,3.5)\n",
    "ax1.set_ylim(-3.5,3.5)\n",
    "\n",
    "# Plot Principal Components 1 and 2\n",
    "for i in df_plot.index:\n",
    "    ax1.annotate(i, (-df_plot.PC1.loc[i], -df_plot.PC2.loc[i]), ha='center')\n",
    "\n",
    "# Plot reference lines\n",
    "ax1.hlines(0,-3.5,3.5, linestyles='dotted', colors='grey')\n",
    "ax1.vlines(0,-3.5,3.5, linestyles='dotted', colors='grey')\n",
    "\n",
    "ax1.set_xlabel('First Principal Component')\n",
    "ax1.set_ylabel('Second Principal Component')\n",
    "    \n",
    "# Plot Principal Component loading vectors, using a second y-axis.\n",
    "ax2 = ax1.twinx().twiny() \n",
    "\n",
    "ax2.set_ylim(-1,1)\n",
    "ax2.set_xlim(-1,1)\n",
    "ax2.set_xlabel('Principal Component loading vectors', color='red')\n",
    "\n",
    "# Plot labels for vectors. Variable 'a' is a small offset parameter to separate arrow tip and text.\n",
    "a = 1.07  \n",
    "for i in pca_loadings[['V1', 'V2']].index:\n",
    "    ax2.annotate(i, (-pca_loadings.V1.loc[i]*a, -pca_loadings.V2.loc[i]*a), color='red')\n",
    "\n",
    "# Plot vectors\n",
    "ax2.arrow(0,0,-pca_loadings.V1[0], -pca_loadings.V2[0])\n",
    "ax2.arrow(0,0,-pca_loadings.V1[1], -pca_loadings.V2[1])\n",
    "ax2.arrow(0,0,-pca_loadings.V1[2], -pca_loadings.V2[2])\n",
    "ax2.arrow(0,0,-pca_loadings.V1[3], -pca_loadings.V2[3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ${\\tt PCA()}$ function also outputs the variance explained by of each principal\n",
    "component. We can access these values as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pca.explained_variance_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also get the proportion of variance explained:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pca.explained_variance_ratio_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the first principal component explains 62.0% of the variance\n",
    "in the data, the next principal component explains 24.7% of the variance,\n",
    "and so forth. We can plot the PVE explained by each component as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(7,5))\n",
    "plt.plot([1,2,3,4], pca.explained_variance_ratio_, '-o')\n",
    "plt.ylabel('Proportion of Variance Explained')\n",
    "plt.xlabel('Principal Component')\n",
    "plt.xlim(0.75,4.25)\n",
    "plt.ylim(0,1.05)\n",
    "plt.xticks([1,2,3,4])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also use the function ${\\tt cumsum()}$, which computes the cumulative sum of the elements of a numeric vector, to plot the cumulative PVE:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(7,5))\n",
    "plt.plot([1,2,3,4], np.cumsum(pca.explained_variance_ratio_), '-s')\n",
    "plt.ylabel('Proportion of Variance Explained')\n",
    "plt.xlabel('Principal Component')\n",
    "plt.xlim(0.75,4.25)\n",
    "plt.ylim(0,1.05)\n",
    "plt.xticks([1,2,3,4])\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 10.6: NCI60 Data Example\n",
    "\n",
    "Let's return to the ${\\tt NCI60}$ cancer cell line microarray data, which\n",
    "consists of 6,830 gene expression measurements on 64 cancer cell lines:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df2 = pd.read_csv('NCI60.csv').drop('Unnamed: 0', axis=1)\n",
    "df2.columns = np.arange(df2.columns.size)\n",
    "df2.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Read in the labels to check our work later\n",
    "y = pd.read_csv('NCI60_y.csv', usecols=[1], skiprows=1, names=['type'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 10.6.1 PCA on the NCI60 Data\n",
    "\n",
    "We first perform PCA on the data after scaling the variables (genes) to\n",
    "have standard deviation one, although one could reasonably argue that it\n",
    "is better not to scale the genes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Scale the data\n",
    "X = pd.DataFrame(scale(df2))\n",
    "X.shape\n",
    "\n",
    "# Fit the PCA model and transform X to get the principal components\n",
    "pca2 = PCA()\n",
    "df2_plot = pd.DataFrame(pca2.fit_transform(X))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now plot the first few principal component score vectors, in order to\n",
    "visualize the data. The observations (cell lines) corresponding to a given\n",
    "cancer type will be plotted in the same color, so that we can see to what\n",
    "extent the observations within a cancer type are similar to each other:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1,2, figsize=(15,6))\n",
    "\n",
    "color_idx = pd.factorize(y.type)[0]\n",
    "cmap = mpl.cm.hsv\n",
    "\n",
    "# Left plot\n",
    "ax1.scatter(df2_plot.iloc[:,0], df2_plot.iloc[:,1], c=color_idx, cmap=cmap, alpha=0.5, s=50)\n",
    "ax1.set_ylabel('Principal Component 2')\n",
    "\n",
    "# Right plot\n",
    "ax2.scatter(df2_plot.iloc[:,0], df2_plot.iloc[:,2], c=color_idx, cmap=cmap, alpha=0.5, s=50)\n",
    "ax2.set_ylabel('Principal Component 3')\n",
    "\n",
    "# Custom legend for the classes (y) since we do not create scatter plots per class (which could have their own labels).\n",
    "handles = []\n",
    "labels = pd.factorize(y.type.unique())\n",
    "norm = mpl.colors.Normalize(vmin=0.0, vmax=14.0)\n",
    "\n",
    "for i, v in zip(labels[0], labels[1]):\n",
    "    handles.append(mpl.patches.Patch(color=cmap(norm(i)), label=v, alpha=0.5))\n",
    "\n",
    "ax2.legend(handles=handles, bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
    "\n",
    "# xlabel for both plots\n",
    "for ax in fig.axes:\n",
    "    ax.set_xlabel('Principal Component 1') "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On the whole, cell lines corresponding to a single cancer type do tend to have similar values on the\n",
    "first few principal component score vectors. This indicates that cell lines\n",
    "from the same cancer type tend to have pretty similar gene expression\n",
    "levels.\n",
    "\n",
    "We can generate a summary of the proportion of variance explained (PVE)\n",
    "of the first few principal components:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pd.DataFrame([df2_plot.iloc[:,:5].std(axis=0, ddof=0).as_matrix(),\n",
    "              pca2.explained_variance_ratio_[:5],\n",
    "              np.cumsum(pca2.explained_variance_ratio_[:5])],\n",
    "             index=['Standard Deviation', 'Proportion of Variance', 'Cumulative Proportion'],\n",
    "             columns=['PC1', 'PC2', 'PC3', 'PC4', 'PC5'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the ${\\tt plot()}$ function, we can also plot the variance explained by the\n",
    "first few principal components:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df2_plot.iloc[:,:10].var(axis=0, ddof=0).plot(kind='bar', rot=0)\n",
    "plt.ylabel('Variances')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, it is generally more informative to\n",
    "plot the PVE of each principal component (i.e. a **scree plot**) and the cumulative\n",
    "PVE of each principal component. This can be done with just a\n",
    "little tweaking:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig , (ax1,ax2) = plt.subplots(1,2, figsize=(15,5))\n",
    "\n",
    "# Left plot\n",
    "ax1.plot(pca2.explained_variance_ratio_, '-o')\n",
    "ax1.set_ylabel('Proportion of Variance Explained')\n",
    "ax1.set_ylim(ymin=-0.01)\n",
    "\n",
    "# Right plot\n",
    "ax2.plot(np.cumsum(pca2.explained_variance_ratio_), '-ro')\n",
    "ax2.set_ylabel('Cumulative Proportion of Variance Explained')\n",
    "ax2.set_ylim(ymax=1.05)\n",
    "\n",
    "for ax in fig.axes:\n",
    "    ax.set_xlabel('Principal Component')\n",
    "    ax.set_xlim(-1,65)    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that together, the first seven principal components\n",
    "explain around 40% of the variance in the data. This is not a huge amount\n",
    "of the variance. However, looking at the scree plot, we see that while each\n",
    "of the first seven principal components explain a substantial amount of\n",
    "variance, there is a marked decrease in the variance explained by further\n",
    "principal components. That is, there is an **elbow** in the plot after approximately\n",
    "the seventh principal component. This suggests that there may\n",
    "be little benefit to examining more than seven or so principal components\n",
    "(phew! even examining seven principal components may be difficult)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}