Theory of PCA - Answer Key

Authors

Noor Sohail

Will Gammerdinger

Published

September 5, 2025

Exercise 1

There are a couple properties of variance and co-variance that we can verify:

cov(X,X) is equal to var(X). We can observe this is mathematically below:

\[ \operatorname{cov}(X, X) = \frac{1}{n-1} \sum_{i=1}^{n} (x_i - \bar{x})(x_i - \bar{x}) = \frac{1}{n-1} \sum_{i=1}^{n} (x_i - \bar{x})^{2} = \operatorname{var}(X) \]

As a result, you will sometimes see covariance matrices written as:

\[ \begin{bmatrix} \operatorname{var}(X) & \operatorname{cov}(X, Y) & \dots & \operatorname{cov}(X, Z) \\ \operatorname{cov}(Y, X) & \operatorname{var}(Y) & \dots & \operatorname{cov}(Y, Z) \\ \dots & \dots & \dots & \dots \\ \operatorname{cov}(Z, X) & \operatorname{cov}(Z, Y) & \dots & \operatorname{var}(Z) \end{bmatrix} \]

Confirm this property by estimating the variance for Gene A.

# Estimate the variance for Gene A
var(recentered_expression_matrix[, "Gene_A"])

[1] 799.3333

# Estimate the variance for Gene A
print(recentered_expression_df["Gene_A"].var())

799.3333333333334

Now estimate the covariance for Gene A and Gene A

# Estimate the covariance for Gene A and Gene A
cov(recentered_expression_matrix[, "Gene_A"], recentered_expression_matrix[, "Gene_A"])

[1] 799.3333

# Estimate the covariance for Gene A and Gene A
print(recentered_expression_df["Gene_A"].cov(recentered_expression_df["Gene_A"]))

799.3333333333333

Is the value the same? Does it match the value in the covariance matrix for Gene A and Gene A?

# Extract the covariance estimate of Gene A and Gene A from the covariance matrix
cov_matrix["Gene_A", "Gene_A"]

[1] 799.3333

# Extract the covariance estimate of Gene A and Gene A from the covariance matrix
print(cov_matrix.loc["Gene_A", "Gene_A"])

799.3333333333333

cov(X,Y) is equal to cov(Y,X). We can observe this is mathematically below:

\[ \operatorname{cov}(X, Y) = \frac{1}{n-1} \sum_{i=1}^{n} (x_i - \bar{x})(y_i - \bar{y}) = \frac{1}{n-1} \sum_{i=1}^{n} (y_i - \bar{y})(x_i - \bar{x}) = \operatorname{cov}(Y, X) \]

Estimate the covariance of Gene B and Gene A.

# Estimate the covariance of Gene B and Gene A
cov(recentered_expression_tibble[, "Gene_B"], recentered_expression_tibble[, "Gene_A"])

         Gene_A
Gene_B 584.6667

# Estimate the covariance of Gene B and Gene A
print(recentered_expression_df["Gene_B"].cov(recentered_expression_df["Gene_A"]))

584.6666666666666

How does this compare to the covariance that we estimated by hand?

It is the same.

Exercise 2

When looking at the percent explained by each principal component, the first principal component should explain the most and each of the following principal components should explain less than the previous principal component. Let’s have a look at our pct_var_explained object, are our results congruent with this expectation?

# Print the percent variance explained
pct_var_explained

    PC_1     PC_2 
96.16723  3.83277

# Print the percent variance explained
print(pct_var_explained)

[96.16723045  3.83276955]

Yes, PC_1 explains python pct_var_explained[0] and PC_2 explains python pct_var_explained[1]

Exercise 3

Create a plot of the Principal Components Analysis derived from prcomp()/PCA(). Is it the same as the plot we derived except only reflected over the x-axis (R)/rotated 180°(Python)?

# Create a tibble to hold the PC scores prcomp() found and also make the Cell IDs into a column
prcomp_pc_scores_tibble <- prcomp_PCA$x %>% 
  as.data.frame() %>% 
  rownames_to_column("cells") %>% 
  as_tibble()

# Plot the PC scores found by prcomp()
ggplot(prcomp_pc_scores_tibble, aes(x = PC1, y = PC2, label = cells)) +
  geom_point( color = "cornflowerblue") +
  geom_text(hjust = 0, vjust = -1) +
  theme_bw() +
  xlim(-50, 50) +
  ylim(-12, 12) +
  xlab(paste0("PC 1 (Variance Explained ", round(prcomp_eigenvalues["PC_1"]/sum(prcomp_eigenvalues) * 100, digits = 2),"%)")) +
  ylab(paste0("PC 2 (Variance Explained ", round(prcomp_eigenvalues["PC_2"]/sum(prcomp_eigenvalues) * 100, digits = 2),"%)")) +
  ggtitle("PCA of Expression Values from Four Cells") +
  theme(plot.title = element_text(hjust = 0.5))

Figure 1: Scatterplot of `PC_1` vs. `PC_2` from `prcomp()`.

Yes, it is the same plot just reflected over the x-axis.

sns.set(style="whitegrid")

# Initialize a plot with a specific size
plt.figure(figsize = (8, 6))

# Add a scatterplot layer to the plot, coloring points by genotype
PCA_sklearn_plot = sns.scatterplot(data = pc_scores_sklearn_df,
                x = "PC_1", 
                y = "PC_2")

# Set x-axis limits
plt.xlim(left = -55,
         right = 55)

# Set y-axis limits
plt.ylim(bottom = -15,
         top = 15)

# Add cell labels at PC coordinates
for cell_id, row in pc_scores_sklearn_df.iterrows():
    PCA_sklearn_plot.text(
        x = row["PC_1"] + 1,
        y = row["PC_2"] + 1,
        s = cell_id,
        fontsize = 9,
        ha = "left",
        va = "bottom"
    )

# Change the text of the x-axis label
plt.xlabel(xlabel = f"PC 1 (Variance Explained {pca_sklearn.explained_variance_ratio_[0].round(4) * 100}%)")

# Change the text of the y-axis label
plt.ylabel(ylabel = f"PC 2 (Variance Explained {pca_sklearn.explained_variance_ratio_[1].round(4) * 100}%)")

# Add plot title
plt.title(label = "PCA of Expression Values from Four Cells")

# Render the plot
plt.show()

Figure 2: Scatterplot of `PC_1` vs. `PC_2` from scikit-learn.

Yes, it is the same plot just rotated 180°.

--- title: "Theory of PCA - Answer Key" author: - Noor Sohail - Will Gammerdinger date: "2025-09-05" --- ```{r} #| label: load_data_R #| echo: false # Items to pre-load library(tidyverse) # Create a vector for Cell IDs cells <- c("Cell_1", "Cell_2", "Cell_3", "Cell_4") # Create a vector to hold expression values for Gene A across all of the cells Gene_A <- c(0, 12, 65, 23) # Create a vector to hold expression values for Gene B across all of the cells Gene_B <- c(4, 30, 57, 18) # Create a tibble to hold the cell names and expression values expression_tibble <- tibble(cells, Gene_A, Gene_B) # Determine the center of the data by: # Finding the average expression of gene A Gene_A_mean <- mean(expression_tibble$Gene_A) # Finding the average expression of gene B Gene_B_mean <- mean(expression_tibble$Gene_B) # Create a vector to hold the center of the data center_of_data <- c(Gene_A_mean, Gene_B_mean) # Assign names to the components of the vector names(center_of_data) <- c("Gene_A", "Gene_B") # Shift the data points so that they data is centered on the origin recentered_expression_tibble <- expression_tibble %>% mutate( Gene_A = Gene_A - Gene_A_mean, Gene_B = Gene_B - Gene_B_mean ) # Move the cell IDs to the rownames and convert the tibble to a matrix recentered_expression_matrix <- recentered_expression_tibble %>% column_to_rownames("cells") %>% as.matrix() # Create a covariance matrix cov_matrix <- cov(recentered_expression_matrix) # Find the eigenvalues and eigenvectors of the covariance matrix eig <- eigen(cov_matrix) # Transform the data into PC space by multiply the re-centered expression matrix by the eigenvectors pc_scores <- recentered_expression_matrix %*% eig$vectors # Name the columns in pc_scores object colnames(pc_scores) <- c("PC_1", "PC_2") # Calculate the percent of variance explained by each PC using the eigenvalues pct_var_explained <- (eig$values / sum(eig$values)) * 100 # Name the elements of the pct_var_explained by their PC names(pct_var_explained) <- c("PC_1", "PC_2") # Create a tibble to hold the PC scores we found and also make the Cell IDs into a column pc_scores_tibble <- pc_scores %>% as.data.frame() %>% rownames_to_column("cells") %>% as_tibble() # Run prcomp() on the expression tibble after moving the Cell IDs to be rownames prcomp_PCA <- expression_tibble %>% column_to_rownames("cells") %>% prcomp() # Print the eigenvalues found by prcomp() by squaring prcomp_PCA$sdev prcomp_eigenvalues <- prcomp_PCA$sdev ** 2 # Name the elements of the prcomp_eigenvalues by their PC names(prcomp_eigenvalues) <- c("PC_1", "PC_2") ``` ```{python} #| label: load_data_Python #| echo: false # Import libraries import numpy as np import pandas as pd import matplotlib.pyplot as plt import seaborn as sns from sklearn.decomposition import PCA # Create a list for Cell IDs cells = ["Cell_1", "Cell_2", "Cell_3", "Cell_4"] # Create a list to hold expression values for Gene A across all of the cells Gene_A = [0, 12, 65, 23] # Create a list to hold expression values for Gene B across all of the cells Gene_B = [4, 30, 57, 18] # Create a DataFrame to hold the cell names and expression values expression_df = pd.DataFrame({ "cells": cells, "Gene_A": Gene_A, "Gene_B": Gene_B }).set_index("cells") # Determine the center of the data by: # Finding the average expression of gene A Gene_A_mean = expression_df["Gene_A"].mean() # Finding the average expression of gene B Gene_B_mean = expression_df["Gene_B"].mean() center_of_data = { "Gene_A": Gene_A_mean, "Gene_B": Gene_B_mean, } # Shift the data so it is centered on the origin recentered_expression_df = expression_df.assign( Gene_A = expression_df["Gene_A"] - Gene_A_mean, Gene_B = expression_df["Gene_B"] - Gene_B_mean, ) # Create a covariance matrix cov_matrix = recentered_expression_df.cov() # Find the eigenvalues and eigenvectors of the covariance matrix eig_values, eig_vectors = np.linalg.eigh(cov_matrix) # Sort the eigenvalues and eigenvectors by descending variance explained # Indices for descending order idx = np.argsort(eig_values)[::-1] # Reordered eigenvalues eig_values = eig_values[idx] # Reordered eigenvectors eig_vectors = eig_vectors[:, idx] # Place the eigenvectors into a dataframe eig_vectors_df = pd.DataFrame( eig_vectors ) pc_scores = recentered_expression_df.values.dot(eig_vectors) pc_scores_df = pd.DataFrame( pc_scores, index = recentered_expression_df.index, columns = ["PC_1", "PC_2"] ) # Calculate the percent of variance explained by each PC using the eigenvalues pct_var_explained = (eig_values / eig_values.sum()) * 100 # Create PCA object pca_sklearn = PCA() # Run PCA on our expression data pc_scores_sklearn = pca_sklearn.fit_transform(expression_df) # Print the eigenvectors found by scikit-learn eigenvector_sklearn = pd.DataFrame( pca_sklearn.components_.T ) pc_scores_sklearn_df = pd.DataFrame( pc_scores_sklearn, index = expression_df.index, columns = ["PC_1", "PC_2"] ) ``` ## Exercise 1 There are a couple properties of variance and co-variance that we can verify: *cov(X,X)* is equal to *var(X)*. We can observe this is mathematically below: $$ \operatorname{cov}(X, X) = \frac{1}{n-1} \sum_{i=1}^{n} (x_i - \bar{x})(x_i - \bar{x}) = \frac{1}{n-1} \sum_{i=1}^{n} (x_i - \bar{x})^{2} = \operatorname{var}(X) $$ As a result, you will sometimes see covariance matrices written as: $$ \begin{bmatrix} \operatorname{var}(X) & \operatorname{cov}(X, Y) & \dots & \operatorname{cov}(X, Z) \\ \operatorname{cov}(Y, X) & \operatorname{var}(Y) & \dots & \operatorname{cov}(Y, Z) \\ \dots & \dots & \dots & \dots \\ \operatorname{cov}(Z, X) & \operatorname{cov}(Z, Y) & \dots & \operatorname{var}(Z) \end{bmatrix} $$ 1. Confirm this property by estimating the variance for Gene A. ::: {.panel-tabset group="language"} ### R ```{r} #| label: covariance_check_1_R # Estimate the variance for Gene A var(recentered_expression_matrix[, "Gene_A"]) ``` ### Python ```{python} #| label: covariance_check_1_Python # Estimate the variance for Gene A print(recentered_expression_df["Gene_A"].var()) ``` ::: 2. Now estimate the covariance for Gene A and Gene A ::: {.panel-tabset group="language"} ### R ```{r} #| label: covariance_check_2_R # Estimate the covariance for Gene A and Gene A cov(recentered_expression_matrix[, "Gene_A"], recentered_expression_matrix[, "Gene_A"]) ``` ### Python ```{python} #| label: covariance_check_2_Python # Estimate the covariance for Gene A and Gene A print(recentered_expression_df["Gene_A"].cov(recentered_expression_df["Gene_A"])) ``` ::: 3. Is the value the same? Does it match the value in the covariance matrix for Gene A and Gene A? ::: {.panel-tabset group="language"} ### R ```{r} #| label: covariance_matrix_check_R # Extract the covariance estimate of Gene A and Gene A from the covariance matrix cov_matrix["Gene_A", "Gene_A"] ``` ### Python ```{python} #| label: covariance_matrix_check_Python # Extract the covariance estimate of Gene A and Gene A from the covariance matrix print(cov_matrix.loc["Gene_A", "Gene_A"]) ``` ::: 4. *cov(X,Y)* is equal to *cov(Y,X)*. We can observe this is mathematically below: $$ \operatorname{cov}(X, Y) = \frac{1}{n-1} \sum_{i=1}^{n} (x_i - \bar{x})(y_i - \bar{y}) = \frac{1}{n-1} \sum_{i=1}^{n} (y_i - \bar{y})(x_i - \bar{x}) = \operatorname{cov}(Y, X) $$ Estimate the covariance of Gene B and Gene A. ::: {.panel-tabset group="language"} ### R ```{r} #| label: covariance_check_3_R # Estimate the covariance of Gene B and Gene A cov(recentered_expression_tibble[, "Gene_B"], recentered_expression_tibble[, "Gene_A"]) ``` ### Python ```{python} #| label: covariance_check_3_Python # Estimate the covariance of Gene B and Gene A print(recentered_expression_df["Gene_B"].cov(recentered_expression_df["Gene_A"])) ``` ::: 5. How does this compare to the covariance that we estimated by hand? It is the same. ## Exercise 2 When looking at the percent explained by each principal component, the first principal component should explain the most and each of the following principal components should explain less than the previous principal component. Let's have a look at our `pct_var_explained` object, are our results congruent with this expectation? ::: {.panel-tabset group="language"} ### R ```{r} #| label: pct_explained_R # Print the percent variance explained pct_var_explained ``` ### Python ```{python} #| label: pct_explained_Python # Print the percent variance explained print(pct_var_explained) ``` ::: Yes, `PC_1` explains `python pct_var_explained[0]` and `PC_2` explains `python pct_var_explained[1]` ## Exercise 3 Create a plot of the Principal Components Analysis derived from `prcomp()`/`PCA()`. Is it the same as the plot we derived except only reflected over the x-axis (R)/rotated 180°(Python)? ::: {.panel-tabset group="language"} ### R ```{r} #| label: fig-plotting_PCA_prcomp_R #| fig-cap: "Scatterplot of `PC_1` vs. `PC_2` from `prcomp()`." # Create a tibble to hold the PC scores prcomp() found and also make the Cell IDs into a column prcomp_pc_scores_tibble <- prcomp_PCA$x %>% as.data.frame() %>% rownames_to_column("cells") %>% as_tibble() # Plot the PC scores found by prcomp() ggplot(prcomp_pc_scores_tibble, aes(x = PC1, y = PC2, label = cells)) + geom_point( color = "cornflowerblue") + geom_text(hjust = 0, vjust = -1) + theme_bw() + xlim(-50, 50) + ylim(-12, 12) + xlab(paste0("PC 1 (Variance Explained ", round(prcomp_eigenvalues["PC_1"]/sum(prcomp_eigenvalues) * 100, digits = 2),"%)")) + ylab(paste0("PC 2 (Variance Explained ", round(prcomp_eigenvalues["PC_2"]/sum(prcomp_eigenvalues) * 100, digits = 2),"%)")) + ggtitle("PCA of Expression Values from Four Cells") + theme(plot.title = element_text(hjust = 0.5)) ``` Yes, it is the same plot just reflected over the x-axis. ### Python ```{python} #| label: plotting_PCA_PCA_Python_show #| eval: false #| echo: true sns.set(style="whitegrid") # Initialize a plot with a specific size plt.figure(figsize = (8, 6)) # Add a scatterplot layer to the plot, coloring points by genotype PCA_sklearn_plot = sns.scatterplot(data = pc_scores_sklearn_df, x = "PC_1", y = "PC_2") # Set x-axis limits plt.xlim(left = -55, right = 55) # Set y-axis limits plt.ylim(bottom = -15, top = 15) # Add cell labels at PC coordinates for cell_id, row in pc_scores_sklearn_df.iterrows(): PCA_sklearn_plot.text( x = row["PC_1"] + 1, y = row["PC_2"] + 1, s = cell_id, fontsize = 9, ha = "left", va = "bottom" ) # Change the text of the x-axis label plt.xlabel(xlabel = f"PC 1 (Variance Explained {pca_sklearn.explained_variance_ratio_[0].round(4) * 100}%)") # Change the text of the y-axis label plt.ylabel(ylabel = f"PC 2 (Variance Explained {pca_sklearn.explained_variance_ratio_[1].round(4) * 100}%)") # Add plot title plt.title(label = "PCA of Expression Values from Four Cells") # Render the plot plt.show() ``` ```{python} #| label: fig-plotting_PCA_PCA_Python_hidden #| fig-cap: "Scatterplot of `PC_1` vs. `PC_2` from scikit-learn." #| echo: false def plot_pca_sklearn(pc_scores_sklearn_df, pca_sklearn): sns.set(style="whitegrid") # Initialize a plot with a specific size fig, ax = plt.subplots(figsize=(8, 6)) # Add a scatterplot layer to the plot sns.scatterplot( data=pc_scores_sklearn_df, x="PC_1", y="PC_2", ax=ax ) # Set axis limits via the Axes object ax.set_xlim(-55, 55) ax.set_ylim(-15, 15) # Add cell labels at PC coordinates for cell_id, row in pc_scores_sklearn_df.iterrows(): ax.text( x=row["PC_1"] + 1, y=row["PC_2"] + 1, s=cell_id, fontsize=9, ha="left", va="bottom" ) # Axis labels with variance explained from sklearn PCA ax.set_xlabel( f"PC 1 (Variance Explained {pca_sklearn.explained_variance_ratio_[0].round(4) * 100}%)" ) ax.set_ylabel( f"PC 2 (Variance Explained {pca_sklearn.explained_variance_ratio_[1].round(4) * 100}%)" ) # Title ax.set_title("PCA of Expression Values from Four Cells") return fig, ax fig, ax = plot_pca_sklearn(pc_scores_sklearn_df, pca_sklearn) plt.show() ``` Yes, it is the same plot just rotated 180°. :::