PCA_II.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Copy of 2021_Spring_BUMK744_PCA_II.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MeA9MTo0qb7t"
      },
      "source": [
        "# **Principal Component Analysis (PCA) - II**\n",
        "\n",
        "\n",
        "**objectives:**\n",
        "\n",
        "In this Jupyter Notebook we continue with **Principal Component Analysis (PCA)**. The general learning goal is twofold:\n",
        "\n",
        "1. Deeper understanding of principal component scores and principal directions\n",
        "2. Understand how much variation each principal component explains, and how to reduce the dimensionality based on this measure\n",
        "\n",
        "> - [A] Importing\n",
        "> - [B] Centering (demeaning) the variables in your data\n",
        "> - [C] PCA on the 2 length variables in the Iris dataset\n",
        "> - [D] Properties of the principal component scores\n",
        "> - [E] PCA on the full Iris data set\n",
        "> - [F] Explained variance by each principal component\n",
        "> - [G] [G] Dimensionality reduction using PCA\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kROvkX_4SGiY"
      },
      "source": [
        "---\n",
        "## **[A] Importing**"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "N9fyPyZJsc7j"
      },
      "source": [
        "We import the same packages and data as in the previous flie:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "E5rxwtMjsl5l"
      },
      "source": [
        "import numpy as np\n",
        "import pandas as pd\n",
        "\n",
        "import matplotlib.pyplot as plt\n",
        "import seaborn as sns\n",
        "sns.set()\n",
        "\n",
        "from sklearn.decomposition import PCA\n",
        "from sklearn import datasets\n",
        "\n",
        "iris = datasets.load_iris()\n",
        "\n",
        "# Load X, y, variable labels for the Iris data set\n",
        "X = iris.data\n",
        "y = iris.target\n",
        "variable_labels = iris.feature_names\n",
        "\n",
        "# Create a dataframe for X\n",
        "X = pd.DataFrame(X, columns=variable_labels)\n",
        "X_only_length = X.loc[:, ['sepal length (cm)', 'petal length (cm)']]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "POxF2oI3UraR"
      },
      "source": [
        "---\n",
        "## **[B] Centering (demeaning) the variables in your data**"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MCU6p_4SUv70"
      },
      "source": [
        "In the previous flie we have seen that the PCA functionality in scikit-learn automatically centers (demeans) the variables in your data, such that each column has zero mean.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8qREmM-QVFkG"
      },
      "source": [
        "X_only_length_centered = X_only_length - X_only_length.mean()\n",
        "X_only_length_centered.mean().round(6)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9BeH9DOsVIUD"
      },
      "source": [
        "A benefit of centering your data is that it facilitates interpretation.\n",
        "- A **positive value** for an (observation, variable) pair indicates that this observation has an **above average** value\n",
        "- A **negative value** for an (observation, variable) pair indicates that this observation has a **below average** value\n",
        "\n",
        "For example, consider the first observation (index 0) after centering:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "HFOw2TGSVMiD"
      },
      "source": [
        "X_only_length_centered.iloc[0]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Mb7rhO_6VTKP"
      },
      "source": [
        "Both values are negative, implying that the 1st observation in our data has a below average sepal length and petal length. We can manually confirm this:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "4gQMjQwTVrda"
      },
      "source": [
        "X_only_length.iloc[0]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qdDBV12FVtkv"
      },
      "source": [
        "Comparing that against the average values:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "UZiHeaYaVxTz"
      },
      "source": [
        "X_only_length.mean()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kbEgADg3L8yq"
      },
      "source": [
        "Indeed this shows that both values are below average.\n",
        "\n",
        "Note that centering your data does **not** change the variance in your variables: "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "jGeiecGcMCNs"
      },
      "source": [
        "X_only_length.var()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "FxubEbQQMEhh"
      },
      "source": [
        "X_only_length_centered.var()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8tlNLNzXVzru"
      },
      "source": [
        "Because working with centered data is so useful, and to avoid any confusion, we will **overwrite** our `X` and `X_only_length` variables to be centered as well in the rest of this notebook:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "d5HZ1F8w6FNs"
      },
      "source": [
        "X = X - X.mean()\n",
        "X_only_length = X_only_length - X_only_length.mean()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "njgjlcLKsZK3"
      },
      "source": [
        "---\n",
        "## **[C] PCA on the 2 length variables in the Iris dataset**"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2HO15dK9Xkzo"
      },
      "source": [
        "We continue were we left of in the previous notebook, by running PCA on the dataset `X_only_length` that only contains the two length variables from the Iris dataset:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Ni5yYc6xEnuh"
      },
      "source": [
        "pca_only_length = PCA()\n",
        "pca_only_length = pca_only_length.fit(X_only_length)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LdcpJ5W-CSNg"
      },
      "source": [
        "The **principal component directions**:\n",
        "- Each direction is expressed as a linear combination of the **original variables**\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7pqngLzVCTO4"
      },
      "source": [
        "pc_directions = pd.DataFrame(\n",
        "    pca_only_length.components_,\n",
        "    index=['PC1', 'PC2'],\n",
        "    columns = X_only_length.columns\n",
        ")\n",
        "pc_directions"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "trpHwysmFULb"
      },
      "source": [
        "The **principal component scores**, for each of the 150 flowers in the Iris dataset:\n",
        "- The principal component scores for each observation is now expressed in terms of the two **principal components**"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "b86faQbTEp6-"
      },
      "source": [
        "pc_scores = pd.DataFrame(\n",
        "    pca_only_length.transform(X_only_length),\n",
        "    columns=['PC1', 'PC2'],\n",
        ")\n",
        "pc_scores"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PCtSeTPfGbiW"
      },
      "source": [
        "---\n",
        "## **[D] Properties of the principal component scores**"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rEPFnXPHGhBb"
      },
      "source": [
        "The principal component scores are designed such that they have a number of different properties:"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uu2dR7lQDfLi"
      },
      "source": [
        "- Principal component scores are **centered** (have zero mean)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "FVZljekl1DsM"
      },
      "source": [
        "pc_scores.mean().round(6)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "o4aSlWmkFP6Z"
      },
      "source": [
        "- Principal component scores are **sorted in descending order of variance**"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "36QTKtVZFS58"
      },
      "source": [
        "pc_scores.var()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eSyxcEmCE05e"
      },
      "source": [
        "- Principal component scores are **uncorrelated**\n",
        " - By design, they are **efficient** in the sense that PC1 does not contain any information about PC2"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "WsbZ5mZ1Dr_t"
      },
      "source": [
        "pc_scores.cov().round(6)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "SKNNf13ZE3EA"
      },
      "source": [
        "- Capture all the variation in the two length variables\n",
        " - Note: This only holds if we keep **all** principal components\n",
        " - That is, if $L = P$ (we will get to this later)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Sj0y92ZRE8fL"
      },
      "source": [
        "np.sum(X_only_length.var())"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9ynaOkFzE_4_"
      },
      "source": [
        "np.sum(pc_scores.var())"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "JrFpk6mtGNSk"
      },
      "source": [
        "assert np.isclose(np.sum(X_only_length.var()), np.sum(pc_scores.var()))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "10yrgQRSUe7W"
      },
      "source": [
        "### Computing principal component scores"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tcX0g1yWJaYO"
      },
      "source": [
        "Each principal component score is computed as follows.\n",
        "\n",
        "For each observation and each principal component: The score for the principal component is obtained by taking the observation's values of the **original variables** and **multiplying** them with the **coefficients** of the principal component direction.\n",
        "\n",
        "> $\\text{score}_{i,d} = x_{i,1} \\text{direction}_{d,1} + x_{i,2} \\text{direction}_{d, 2}.$\n",
        "\n",
        "For example, consider the 7th observation (at index 6) and the 2nd principal component (at index 1):\n",
        "\n",
        "> $\\text{score}_{7,2} = x_{7,1} \\text{direction}_{2,1} + x_{7,2} \\text{direction}_{2, 2}.$\n",
        "\n",
        "The principal component score:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-4G9qpP0KDUR"
      },
      "source": [
        "np.sum(X_only_length.iloc[6] * pc_directions.iloc[1])"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RcIucBweHlMD"
      },
      "source": [
        "Let's verify that this is equal to score of the 2nd principal component for the 7th observation:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "McRTQHpyLGNK"
      },
      "source": [
        "pc_scores.iloc[6, 1]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "q_UdWlh6LLNf"
      },
      "source": [
        "The above can be summarized using a **vector-vector** multiplication (dot product) between the vector of values for the 7th observation and the direction vector for the 2nd principal component:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "lrzo2ZiTLjBV"
      },
      "source": [
        "X_only_length.iloc[6, :] @ pc_directions.iloc[1, :]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8dN8OZMML1Ta"
      },
      "source": [
        "Similarly, we can get **all** principal component scores for the 7th observation using a **matrix-vector** multiplication between the **`directions`** matrix for **all** principal components and the vector of values for the 7th observation:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "PSswL2N2MViu"
      },
      "source": [
        "pc_directions @ X_only_length.iloc[6, :]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "djdSfdW9MYi_"
      },
      "source": [
        "Again, let's verify:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "mvYUwrLZMZpW"
      },
      "source": [
        "pc_scores.iloc[6]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CY90tMVxMeEJ"
      },
      "source": [
        "Finally, we can compute the principal component scores for **all** observations and **all** principal components using a **matrix-matrix** multiplication between the data matrix and the **transposed** directions matrix:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_3_g2_lXMojr"
      },
      "source": [
        "X_only_length @ pc_directions.T"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gZXUNWE3Myym"
      },
      "source": [
        "Again, let's verify:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "IQ9nzSpnMvbt"
      },
      "source": [
        "pc_scores"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "O7W6De9jXCVp"
      },
      "source": [
        "### Interpreting the principal component scores"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xRLDFJl4YZqU"
      },
      "source": [
        "The next step is to provide an interpretation to the principal component scores. For example:\n",
        "\n",
        "> What does it mean if an observation has a high PC1 score?\n",
        "\n",
        "We will rephrase this question slightly:\n",
        "\n",
        "> **Question**: When do we **expect** an observation to have a large positive PC1 score?\n",
        "\n",
        "To answer this question, we first need to examine the **direction** corresponding to PC1:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "R_whyrqyXK0R"
      },
      "source": [
        "pc_directions.iloc[0]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gqFTiNgxXcsS"
      },
      "source": [
        "Note that both variables **positively** contribute to PC1.\n",
        "\n",
        "We **expect** an observation to have a large positive PC1 score if:\n",
        "- The value corresponding to petal length is **positive**\n",
        "- The value corresponding to sepal length is **positive**\n",
        "\n",
        "However, remember that we have **centered** our data:\n",
        "- Positive values correspond to **above average** values\n",
        "\n",
        "This allows us to provide a **more intuitive** interpretation.\n",
        "\n",
        "Restating the above, we **expect** an observation to have a large positive PC1 score if:\n",
        "- The value corresponding to petal length (cm) is **above average**\n",
        "- The value corresponding to sepal length (cm) is **above average**\n",
        "\n",
        "Also, because the coefficient for petal length is larger (0.919279) than the coefficient for sepal length (0.393606), intuitively petal length will have a **larger influence** on the score for PC1."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6J0N0NO_Y8J0"
      },
      "source": [
        "Let's verify this intuition by creating a scatterplot of the PC1 and PC2 scores.\n",
        "\n",
        "We have done something similar in the previous lecture, but now we will add **index labels** for each observation to the plot as well:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "qkiCLWWGYYJy"
      },
      "source": [
        "fig, ax = plt.subplots(nrows=1, ncols=2, figsize=[20, 10])\n",
        "\n",
        "# Scatterplot between petal length (horizontal) and sepal length (vertical)\n",
        "ax[0].set_title('Scatter plot between sepal length and petal length')\n",
        "ax[0].set_xlabel('Sepal length (centered)')\n",
        "ax[0].set_ylabel('Petal length (centered)')\n",
        "\n",
        "ax[0].axvline(0)\n",
        "ax[0].axhline(0)\n",
        "\n",
        "ax[0] = sns.scatterplot(\n",
        "    x=X_only_length.loc[:, 'sepal length (cm)'],\n",
        "    y=X_only_length.loc[:, 'petal length (cm)'],\n",
        "    ax=ax[0],\n",
        ")\n",
        "ax[0].axis('equal')\n",
        "\n",
        "# Add index label for each observation\n",
        "for row_idx, (row_x, row_y) in X_only_length.iterrows():\n",
        "    ax[0].text(row_x, row_y, row_idx)\n",
        "\n",
        "# Scatterplot between PC1 (horizontal) and PC2 (vertical)\n",
        "plt.title('Scatter plot between PC1 and PC2')\n",
        "plt.xlabel('PC1')\n",
        "plt.ylabel('PC2')\n",
        "\n",
        "ax[1].axvline(0)\n",
        "ax[1].axhline(0)\n",
        "\n",
        "ax[1] = sns.scatterplot(\n",
        "    x=pc_scores.loc[:, 'PC1'],\n",
        "    y=pc_scores.loc[:, 'PC2'],\n",
        "    ax=ax[1],\n",
        ")\n",
        "ax[1].axis('equal')\n",
        "\n",
        "# Add index label for each observation\n",
        "for row_idx, (row_x, row_y) in pc_scores.iterrows():\n",
        "    ax[1].text(row_x, row_y, row_idx)\n",
        "\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "UOCV0PrYZOtt"
      },
      "source": [
        "From the two plot it follows that:\n",
        "- Observations that have an above average sepal length and petal length (for example: observations 118, 122, 105, 131) also have a relatively **large positive** PC1 score\n",
        "\n",
        "- Observations that have a below average sepal length and petal length (for example: observations 13, 22) also have a relatively **large negative** PC2 score\n",
        "\n",
        "In other words, if we want to **describe** PC1, we could state that it measures the **length** of a plant in both the sepal and petal dimension, with the caveat that petal length has a larger influence on this score (in other words, it does not represent just the unweighted combination  of sepal and petal length)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KIzdCxoFdsvG"
      },
      "source": [
        "---\n",
        "\n",
        "We could ask a similar question for PC2:\n",
        "\n",
        "> **Question**: When do we **expect** an observation to have a high PC2 score?\n",
        "\n",
        "Examining the **direction** corresponding to PC2:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Hd8D0539UQKu"
      },
      "source": [
        "pc_directions.iloc[1]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Ucp91I44eIKe"
      },
      "source": [
        "Note that **sepal length** has a **negative** contribution to PC2, while **petal length** has a **positive** contribution to PC2.\n",
        "\n",
        "We **expect** an observation to have a large positive PC2 score if:\n",
        "- The value corresponding to petal length is **above average** (positive)\n",
        "- The value corresponding to sepal length is **below average** (negative)\n",
        "\n",
        "Also, because the coefficient for sepal length is larger in an absolute sense (-0.919279) than the coefficient for petal length (0.393606), intuitively sepal length will have a **larger influence** on the score for PC2.\n",
        "\n",
        "We can again confirm this intuition by examining the scatter plots above."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3yp0Kd3B992H"
      },
      "source": [
        "---\n",
        "## [E] PCA on the full Iris data set"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "--FsiwYx-Baa"
      },
      "source": [
        "Now that we have an understanding of **principal directions** and **principal component scores**, we can move on and apply PCA to the complete Iris data set, containing 4 variables. We will store the results in the `pca_full` object:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ruum_1Nz-Dvf"
      },
      "source": [
        "pca_full = PCA()\n",
        "pca_full.fit(X)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "RJfRrc4t-G4x"
      },
      "source": [
        "pc_full_directions = pd.DataFrame(\n",
        "    pca_full.components_,\n",
        "    index=['PC' + str(i + 1) for i in range(len(X.columns))],\n",
        "    columns = X.columns\n",
        ")\n",
        "pc_full_scores = pd.DataFrame(\n",
        "    pca_full.transform(X),\n",
        "    columns=['PC' + str(i + 1) for i in range(len(X.columns))],\n",
        ")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0HAM7wM2D0Nl"
      },
      "source": [
        "Now that we have a larger matrix with PC directions we should spend some time on how to interpret these numbers, which represent **coefficients**."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ow5LZDmYD043"
      },
      "source": [
        "pc_full_directions"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jBbIri8e_jM5"
      },
      "source": [
        "In thie interpretation we will focus per principal direction, that is, per row:\n",
        "We will focus on the interpretation per principal direction, i.e. per row:\n",
        "\n",
        "**Sign of the coefficients per row**\n",
        "> - If you want to interpret the **relative importance** of a variable for a principal direction, you should ignore the **sign** of the coefficient and only focus on its **absolute value**\n",
        "\n",
        "**Maximum (absolute) value of the coefficient per row**\n",
        "> - The largest coefficient that we can find per row is equal to 1 or -1 (again, ignore the sign here)\n",
        "> - The closer one of the coefficients is to 1 or -1, the more influential that variable is\n",
        "> - If one of the coefficients is exactly equal to 1 or -1, all the other coefficients for that direction should be equal to 0\n",
        "\n",
        "**To interpret a principal direction (each row)**\n",
        "> - Look for variables that have a large coefficient (in absolute sense, so close to 1 or -1), because they have a large influence on the score of the principal component\n",
        "> - Variables with smaller coefficients closer to 0 should be ignored, or receive less emphasis, as these do not contribute as much to the score for that principal direction"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Ncu0HuG5_ixR"
      },
      "source": [
        "pc_full_directions"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Yvr0gTM7BkBe"
      },
      "source": [
        "To provide an interpretation for the principal directions above:\n",
        "> - PC1 is primarily about petal length. A high score for PC1 likely corresponds with a petal length that is above arage. Conversely, a low score for PC1 likely corresponds with a petal length that is below average\n",
        "> - PC2 describes the size of the sepal. A high score for PC2 likely corresponds with a relatively large sepal (above average), while a low score for PC2 likely corresponds with a relatively small sepal (below average).\n",
        "> - PC3 is a little more work to interpret. It does not load high on petal length so we can ignore that variable. It places approximately the same weight (in absolute value) on the other 3 variables in the data. Note that both width variables have a positive sign, while the sepal length variable has a negative sign. This implies that the PC3 score seems to correspond to the ratio of general flower width to sepal length. A high value for PC3 implies that a flower likely has an above average width and a below average sepal length. A low value for PC3 implies that a flower likely has a below average width and an above average sepal length"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hMdmyiYosQjF"
      },
      "source": [
        "---\n",
        "## [F] Explained variance by each principal component"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NPM8NiioGCP9"
      },
      "source": [
        "### The `explained_variance_` attribute\n",
        "\n",
        "The `explained_variance_` attribute of the PCA object gives the **eigenvalue** for each principal component:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "FvSSSeb-DmrM"
      },
      "source": [
        "eigenvalues = pca_full.explained_variance_\n",
        "eigenvalues"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WtyNj-b4Edoq"
      },
      "source": [
        "Each eigenvalue describes **how much variation** is explained by each principal component.\n",
        "\n",
        "- Remember, the principal components are constructed such that the first principal component explains the most variance, the second principal component the second most variance, and so and so forth\n",
        "- This also shows from the eigenvalues, as those are sorted in descending order"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5VDMkzSNE2ZQ"
      },
      "source": [
        "Let's look at the **first eigenvalue**, which corresponds to the **first principal component**"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6kVrhuCqD7j0"
      },
      "source": [
        "eigenvalues[0]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jsn42y_iFJB_"
      },
      "source": [
        "Compare this against the variance of the first principal component score:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "O3cmcEh3FM23"
      },
      "source": [
        "pc_full_scores.loc[:, 'PC1'].var()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LBDqMgp2FTA0"
      },
      "source": [
        "Indeed, the first eigenvalue is identical to the variance of the first principal component score.\n",
        "\n",
        "Similarly, for all eigenvalues and all principal component scores:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "EU2DmoSwFecj"
      },
      "source": [
        "eigenvalues"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "YdxCq9F2FfpX"
      },
      "source": [
        "pc_full_scores.var()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JIjUoizXFosl"
      },
      "source": [
        "In fact, the sum of the eigenvalues is equal to the sum of the variance in the data set:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "zAu_2g3KFuLK"
      },
      "source": [
        "np.sum(eigenvalues)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NGod4YaMD_5c"
      },
      "source": [
        "and:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "gsBCZUarEEB7"
      },
      "source": [
        "np.sum(X.var())"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eEAHXtZ9EGMf"
      },
      "source": [
        "and naturally also identical to the sum of the variance in the principal component scores:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Rrwa32gSERL_"
      },
      "source": [
        "np.sum(pc_full_scores.var())"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FQp0kNzVGTxG"
      },
      "source": [
        "### The `explained_variance_ratio_` attribute\n",
        "\n",
        "The `explained_variance_ratio_` attribute gives the ratio of total variance explained by each principal component:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "qC0dWsfZGf1s"
      },
      "source": [
        "pca_full.explained_variance_ratio_"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OIq3dfmsG_v6"
      },
      "source": [
        "Each of these numbers tell us the ratio of the **total variation** in the data that is explained by each principal component.\n",
        "\n",
        "In this case:\n",
        "- PC1 accounts for approximately 92.46% of all the variation in the data\n",
        "- PC2 accounts for another 5.31%\n",
        "- The remaining variation (about 2%) is explained by PC3 and PC4\n",
        "\n",
        "In other words:\n",
        "- PC3 and PC4 do not add a lot of information\n",
        "- The combination of PC1 and PC2 is able to summarize the variation of **all 4 variables in the data** very well\n",
        "\n",
        "This is what we typically see in PCA solutions. We only need a small number of principal components to summarize a lot of the variation the data, even if the data consists of many (hundreds) of variables."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nm4LFnbIGjr7"
      },
      "source": [
        "Note that the `explained_variance_ratio_` adds up to 1:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "F6F4966SGk4w"
      },
      "source": [
        "np.sum(pca_full.explained_variance_ratio_)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "SJsk4Ry6GnZ8"
      },
      "source": [
        "Note that the `explained_variance_ratio_` attribute is just here for convenience, each element is equal to the corresponding eigenvalue divided by its sum:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "KXTh4GJzGwGT"
      },
      "source": [
        "eigenvalues / np.sum(eigenvalues)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TWYRJtjFIR1v"
      },
      "source": [
        "---\n",
        "## [G] Dimensionality reduction using PCA"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "t1cipDT3IfUf"
      },
      "source": [
        "The **number of principal components that we retain** is given by the number $L$. Setting $L < P$ will lead to **dimensionality reduction**.\n",
        "\n",
        "The **eigenvalues** can directly be used for some (heuristic) tests to determine the number of principal components we would like to retain in our analysis. The intuition here is that we only want to **keep those principal components that explain a lot of variation** in the data.\n",
        "\n",
        "We will discuss 3 heuristic tests to determine $L$:\n",
        "- Kaiser's rule\n",
        "- Scree plot\n",
        "- Proportion of variance explained\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ls0Ppr7UJHBp"
      },
      "source": [
        "### Kaiser's rule"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eFvmLd3KJSKQ"
      },
      "source": [
        "**Kaiser's rule**: Only keep the principal components that have an eigenvalue larger than 1.\n",
        "- Only works well on data that has been standardized (such that all variables have a variance of 1)\n",
        "\n",
        "> Intuition: We only keep principal components that explain the variance of **at least one variable**.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "T9FT5HKTJoaH"
      },
      "source": [
        "eigenvalues"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Cy5Yz8DMJqPL"
      },
      "source": [
        "In this case, using Kaiser's rule we would keep only the first principal component"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vj2LYd6TJxHw"
      },
      "source": [
        "### Scree plot"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EMyLbs0kJ0No"
      },
      "source": [
        "**Scree plot**: Create a **line plot** of the eigenvalues and only retain the principal components **before** the eigenvalues start to level off\n",
        "> Often, the scree plot has the shape of an **elbow**, sometimes the scree plot is also referred to as the \"elbow plot\"\n",
        "\n",
        "> The scree test then corresponds to finding the principal component index that corresponds to the \"bend of the elbow\". This is the point where the eigenvalues start to level off.\n",
        "\n",
        "> The scree test then says that we should select all the principal components **before** the bend of the elbow."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "sRjZezOMJ78o"
      },
      "source": [
        "plt.figure(figsize=(10, 10))\n",
        "\n",
        "plt.title('Scree plot of the eigenvalues')\n",
        "plt.xlabel('Principal component index')\n",
        "plt.ylabel('Eigenvalue')\n",
        "\n",
        "sns.lineplot(\n",
        "    x=np.arange(len(eigenvalues)) + 1,\n",
        "    y=eigenvalues,\n",
        "    marker='o',\n",
        "    markersize=12,\n",
        ")\n",
        "plt.xticks(np.arange(len(eigenvalues)) + 1)\n",
        "plt.axis('equal')\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mZ3LBJX1KHwD"
      },
      "source": [
        "In the scree plot above, the \"bending point of the elbow\" can be clearly defined as occuring **at** the second principal component:\n",
        "\n",
        "> This implies that starting from the second principal components, the eigenvalues level off, and each principal component more or less explains the same amount of variance\n",
        "\n",
        "In this case, using the scree plot, we would only keep the first principal component.\n",
        "\n",
        "**Gotcha #1**: You should be aware that the scree plot does not always represent an elbow! In particular this is the case if PCA is not able to summarize information in just a few components\n",
        "\n",
        "**Gotcha #2**: In this case the Kaiser rule and the scree plot resulted in the same number of selected components. In general, this does not have to be the case!"
      ]
    }
  ]
}