Introduction

Principal Component Analysis (PCA) is one of the most fundamental and widely-used techniques in statistics and data science for dimensionality reduction and exploratory data analysis. In our increasingly data-rich world, we often encounter datasets with many variables (features), making them difficult to visualize, analyze, or even store efficiently. PCA provides an elegant solution by transforming high-dimensional data into a lower-dimensional representation while preserving as much of the original information as possible (James et al. 2013).

What is PCA?

PCA is an unsupervised learning technique that finds the directions (principal components) along which the data varies the most. These directions are linear combinations of the original variables, and they form a new coordinate system that captures the maximum variance in the data with fewer dimensions.

Think of it this way: imagine you’re looking at a 3D object and trying to take the most informative 2D photograph of it. You would position your camera to capture the angle that shows the most detail and variation. PCA does something similar mathematically – it finds the “best angles” to view your high-dimensional data in lower dimensions.

Why Do We Need PCA?

There are several compelling reasons to use PCA:

  1. Dimensionality Reduction: Reduce storage space and computational complexity
  2. Visualization: Plot high-dimensional data in 2D or 3D
  3. Noise Reduction: Focus on the most important patterns in the data
  4. Feature Extraction: Create new features that capture the essence of the original data
  5. Multicollinearity: Handle highly correlated variables in regression analysis

Learning Objectives

By the end of this tutorial, you will be able to:

  • Understand the mathematical foundations of PCA
  • Implement PCA in R using both manual calculations and built-in functions
  • Interpret PCA results including loadings and scores
  • Use explained variance to choose the optimal number of components
  • Apply PCA to both synthetic and real-world datasets
  • Create effective visualizations of PCA results

Mathematical Foundation

The Mathematics Behind PCA

PCA is fundamentally about finding the eigenvectors and eigenvalues of a covariance matrix. Let’s break this down step by step.

Step 1: Data Preparation

Given a dataset \(X\) with \(n\) observations and \(p\) variables, we first center the data by subtracting the mean from each variable:

\[X_{centered} = X - \bar{X}\]

where \(\bar{X}\) is the vector of variable means.

Step 2: Covariance Matrix

We compute the covariance matrix \(C\):

\[C = \frac{1}{n-1}X_{centered}^T X_{centered}\]

Step 3: Eigendecomposition

The principal components are the eigenvectors of the covariance matrix. We solve:

\[C \mathbf{v} = \lambda \mathbf{v}\]

where \(\mathbf{v}\) are the eigenvectors (principal components) and \(\lambda\) are the eigenvalues.

Step 4: Sorting and Selection

We sort the eigenvectors by their corresponding eigenvalues in descending order. The first principal component corresponds to the largest eigenvalue and captures the most variance.

Step 5: Transformation

Finally, we project the original data onto the selected principal components:

\[Y = X_{centered} \mathbf{V}\]

where \(\mathbf{V}\) is the matrix of selected eigenvectors and \(Y\) is the transformed data.

Geometric Interpretation

Geometrically, PCA finds a new coordinate system where:

  1. The first axis (PC1) points in the direction of maximum variance
  2. The second axis (PC2) is orthogonal to PC1 and points in the direction of maximum remaining variance
  3. This process continues for all dimensions

Let’s visualize this concept with a simple 2D example:

# Create synthetic 2D data with correlation
set.seed(123)
n <- 100
x1 <- rnorm(n)
x2 <- 0.7 * x1 + 0.5 * rnorm(n)  # Create correlation
data_2d <- data.frame(x1 = x1, x2 = x2)

# Perform PCA
pca_2d <- prcomp(data_2d, center = TRUE, scale. = FALSE)

# Extract PC directions for plotting
pc1_slope <- pca_2d$rotation[2,1] / pca_2d$rotation[1,1]
pc2_slope <- pca_2d$rotation[2,2] / pca_2d$rotation[1,2]

# Create the plot
ggplot(data_2d, aes(x = x1, y = x2)) +
  geom_point(alpha = 0.6, size = 2) +
  geom_abline(slope = pc1_slope, intercept = 0, 
              color = "red", size = 1.2, 
              linetype = "dashed") +
  geom_abline(slope = pc2_slope, intercept = 0, 
              color = "blue", size = 1.2, 
              linetype = "dashed") +
  annotate("text", x = 2, y = 1.5, label = "PC1", 
           color = "red", size = 4, fontface = "bold") +
  annotate("text", x = -1.5, y = 1, label = "PC2", 
           color = "blue", size = 4, fontface = "bold") +
  labs(title = "Geometric Interpretation of PCA",
       subtitle = "Red line shows PC1 (maximum variance), blue line shows PC2",
       x = "X1", y = "X2") +
  coord_fixed() +
  theme_minimal()
Geometric interpretation of PCA showing original data points and principal component directions

Geometric interpretation of PCA showing original data points and principal component directions

Synthetic Data Example

Let’s start with a controlled synthetic example to understand how PCA works step by step.

Creating Synthetic Data

We’ll create a 3D dataset where we know the underlying structure:

# Set seed for reproducibility
set.seed(42)

# Create synthetic 3D data with known structure
n <- 200
t <- seq(0, 4*pi, length.out = n)

# Original data lies approximately on a 2D manifold in 3D space
x <- cos(t) + 0.1 * rnorm(n)
y <- sin(t) + 0.1 * rnorm(n)  
z <- t/4 + 0.05 * rnorm(n)

# Create additional correlated variables
x2 <- x + 0.5 * y + 0.2 * rnorm(n)
y2 <- y - 0.3 * x + 0.2 * rnorm(n)

synthetic_data <- data.frame(
  Var1 = x,
  Var2 = y, 
  Var3 = z,
  Var4 = x2,
  Var5 = y2
)

head(synthetic_data)
##       Var1        Var2        Var3      Var4        Var5
## 1 1.137096 -0.20009292  0.06674563 0.9873528 -0.40346013
## 2 0.941537  0.09648335 -0.02767669 1.0742428 -0.04096116
## 3 1.028348  0.24309222  0.03434814 1.3474250 -0.02193620
## 4 1.045396  0.39426559  0.04981404 1.4096420  0.04031557
## 5 1.008695  0.11222680  0.03422980 0.9327042 -0.46351971
## 6 0.939955  0.19543248  0.02899756 1.3504852 -0.14834156

Exploring the Data

Let’s examine the correlation structure and distributions:

# Correlation matrix
cor_matrix <- cor(synthetic_data)
corrplot(cor_matrix, method = "color", type = "upper", 
         order = "hclust", tl.cex = 0.8, tl.col = "black")
Correlation matrix and pairwise relationships in synthetic data

Correlation matrix and pairwise relationships in synthetic data 

# Summary statistics
summary(synthetic_data)
##       Var1                Var2                Var3               Var4         
##  Min.   :-1.166647   Min.   :-1.183752   Min.   :-0.02768   Min.   :-1.72432  
##  1st Qu.:-0.726299   1st Qu.:-0.709321   1st Qu.: 0.78040   1st Qu.:-0.77150  
##  Median : 0.005275   Median :-0.038220   Median : 1.58409   Median :-0.01383  
##  Mean   : 0.002252   Mean   : 0.001128   Mean   : 1.56791   Mean   :-0.02279  
##  3rd Qu.: 0.685007   3rd Qu.: 0.691908   3rd Qu.: 2.36661   3rd Qu.: 0.69961  
##  Max.   : 1.179530   Max.   : 1.180658   Max.   : 3.19294   Max.   : 1.46789  
##       Var5         
##  Min.   :-1.53951  
##  1st Qu.:-0.58726  
##  Median :-0.03125  
##  Mean   : 0.01502  
##  3rd Qu.: 0.63215  
##  Max.   : 1.30652

Performing PCA on Synthetic Data

Now let’s perform PCA and examine the results:

# Perform PCA
pca_synthetic <- prcomp(synthetic_data, center = TRUE, scale. = TRUE)

# Display the summary
summary(pca_synthetic)
## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5
## Standard deviation     1.5262 1.3485 0.8841 0.21248 0.15892
## Proportion of Variance 0.4659 0.3637 0.1563 0.00903 0.00505
## Cumulative Proportion  0.4659 0.8296 0.9859 0.99495 1.00000
# Principal component loadings (rotation matrix)
print("Principal Component Loadings:")
## [1] "Principal Component Loadings:"
print(round(pca_synthetic$rotation, 3))
##         PC1    PC2    PC3    PC4    PC5
## Var1 -0.168 -0.712 -0.001  0.263  0.629
## Var2 -0.608  0.192 -0.252 -0.651  0.326
## Var3  0.367 -0.102 -0.924 -0.020 -0.011
## Var4 -0.429 -0.550 -0.099 -0.075 -0.705
## Var5 -0.532  0.379 -0.268  0.708 -0.010

Visualizing PCA Results

# Create a comprehensive visualization
par(mfrow = c(2, 2))

# 1. Scree plot
var_explained <- (pca_synthetic$sdev^2) / sum(pca_synthetic$sdev^2)
barplot(var_explained, names.arg = paste0("PC", 1:5), 
        main = "Scree Plot - Proportion of Variance Explained",
        ylab = "Proportion of Variance", ylim = c(0, 0.6))

# 2. Cumulative variance plot  
cumvar <- cumsum(var_explained)
plot(1:5, cumvar, type = "b", pch = 19, 
     main = "Cumulative Variance Explained",
     xlab = "Number of Components", 
     ylab = "Cumulative Proportion of Variance",
     ylim = c(0, 1))
abline(h = 0.95, col = "red", lty = 2)
text(3, 0.97, "95% threshold", col = "red")

# 3. PC1 vs PC2 scores plot
plot(pca_synthetic$x[,1], pca_synthetic$x[,2], 
     pch = 19, col = rainbow(n)[order(synthetic_data$Var3)],
     main = "PC1 vs PC2 Scores", 
     xlab = paste0("PC1 (", round(var_explained[1]*100, 1), "%)"),
     ylab = paste0("PC2 (", round(var_explained[2]*100, 1), "%)"))

# 4. Biplot showing both scores and loadings
biplot(pca_synthetic, scale = 0, cex = 0.6,
       main = "PCA Biplot")
PCA results for synthetic data showing scores plots and explained variance

PCA results for synthetic data showing scores plots and explained variance 

par(mfrow = c(1, 1))

Interpreting the Results

From our synthetic data analysis:

  1. PC1 explains 46.6% of the variance
  2. PC1 + PC2 together explain 83% of the variance
  3. The first 3 components capture 98.6% of the total variance

This demonstrates how PCA can effectively reduce dimensionality from 5 variables to 2-3 principal components while retaining most of the information.

Real Dataset Example: The Iris Dataset

Now let’s apply PCA to a classic real-world dataset – the famous iris dataset (Anderson 2003).

Loading and Exploring the Iris Data

# Load the iris dataset
data(iris)

# Remove species column for PCA (we'll use it for coloring later)
iris_numeric <- iris[, 1:4]

# Explore the data
head(iris)
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa
summary(iris_numeric)
##   Sepal.Length    Sepal.Width     Petal.Length    Petal.Width   
##  Min.   :4.300   Min.   :2.000   Min.   :1.000   Min.   :0.100  
##  1st Qu.:5.100   1st Qu.:2.800   1st Qu.:1.600   1st Qu.:0.300  
##  Median :5.800   Median :3.000   Median :4.350   Median :1.300  
##  Mean   :5.843   Mean   :3.057   Mean   :3.758   Mean   :1.199  
##  3rd Qu.:6.400   3rd Qu.:3.300   3rd Qu.:5.100   3rd Qu.:1.800  
##  Max.   :7.900   Max.   :4.400   Max.   :6.900   Max.   :2.500
# Check for missing values
sum(is.na(iris_numeric))
## [1] 0

Correlation Analysis

# Correlation matrix
iris_cor <- cor(iris_numeric)
corrplot(iris_cor, method = "color", type = "upper", 
         order = "hclust", tl.cex = 0.8, tl.col = "black",
         title = "Iris Dataset - Correlation Matrix")
Correlation structure in the iris dataset

Correlation structure in the iris dataset 

# Pairwise scatter plots
GGally::ggpairs(iris_numeric, 
                title = "Pairwise Relationships in Iris Data") +
  theme_minimal()
Correlation structure in the iris dataset

Correlation structure in the iris dataset

Performing PCA on Iris Data

# Perform PCA (scaling is important here due to different units)
pca_iris <- prcomp(iris_numeric, center = TRUE, scale. = TRUE)

# Summary of PCA results
summary(pca_iris)
## Importance of components:
##                           PC1    PC2     PC3     PC4
## Standard deviation     1.7084 0.9560 0.38309 0.14393
## Proportion of Variance 0.7296 0.2285 0.03669 0.00518
## Cumulative Proportion  0.7296 0.9581 0.99482 1.00000
# Loadings
print("Principal Component Loadings:")
## [1] "Principal Component Loadings:"
kable(round(pca_iris$rotation, 3), caption = "PC Loadings for Iris Dataset")
PC Loadings for Iris Dataset
PC1 PC2 PC3 PC4
Sepal.Length 0.521 -0.377 0.720 0.261
Sepal.Width -0.269 -0.923 -0.244 -0.124
Petal.Length 0.580 -0.024 -0.142 -0.801
Petal.Width 0.565 -0.067 -0.634 0.524 
# Calculate variance explained
var_exp_iris <- (pca_iris$sdev^2) / sum(pca_iris$sdev^2)
print(paste("PC1 explains", round(var_exp_iris[1]*100, 1), "% of variance"))
## [1] "PC1 explains 73 % of variance"
print(paste("PC2 explains", round(var_exp_iris[2]*100, 1), "% of variance"))
## [1] "PC2 explains 22.9 % of variance"
print(paste("PC1+PC2 explain", round(sum(var_exp_iris[1:2])*100, 1), "% of variance"))
## [1] "PC1+PC2 explain 95.8 % of variance"

Visualizing Iris PCA Results

# Create data frame for plotting
iris_pca_df <- data.frame(
  PC1 = pca_iris$x[,1],
  PC2 = pca_iris$x[,2], 
  PC3 = pca_iris$x[,3],
  Species = iris$Species
)

# Create multiple plots
p1 <- ggplot(iris_pca_df, aes(x = PC1, y = PC2, color = Species)) +
  geom_point(size = 2, alpha = 0.7) +
  labs(title = "PCA of Iris Dataset",
       subtitle = "Clear separation of species in PC space",
       x = paste0("PC1 (", round(var_exp_iris[1]*100, 1), "%)"),
       y = paste0("PC2 (", round(var_exp_iris[2]*100, 1), "%)")) +
  theme_minimal() +
  stat_ellipse(aes(color = Species), level = 0.95, linetype = 2)

# Scree plot
scree_data <- data.frame(
  Component = factor(1:4, labels = paste0("PC", 1:4)),
  Variance = var_exp_iris,
  Cumulative = cumsum(var_exp_iris)
)

p2 <- ggplot(scree_data, aes(x = Component, y = Variance)) +
  geom_bar(stat = "identity", fill = "steelblue", alpha = 0.7) +
  geom_line(aes(group = 1), color = "red", size = 1) +
  geom_point(color = "red", size = 2) +
  labs(title = "Scree Plot - Iris Dataset",
       y = "Proportion of Variance Explained") +
  theme_minimal()

# Cumulative variance plot
p3 <- ggplot(scree_data, aes(x = Component, y = Cumulative)) +
  geom_bar(stat = "identity", fill = "lightgreen", alpha = 0.7) +
  geom_hline(yintercept = 0.95, color = "red", linetype = "dashed") +
  annotate("text", x = 2, y = 0.97, label = "95% threshold", color = "red") +
  labs(title = "Cumulative Variance Explained",
       y = "Cumulative Proportion") +
  theme_minimal()

# 3D plot (PC1 vs PC2 vs PC3)
p4 <- ggplot(iris_pca_df, aes(x = PC1, y = PC3, color = Species)) +
  geom_point(size = 2, alpha = 0.7) +
  labs(title = "PC1 vs PC3",
       x = paste0("PC1 (", round(var_exp_iris[1]*100, 1), "%)"),
       y = paste0("PC3 (", round(var_exp_iris[3]*100, 1), "%)")) +
  theme_minimal()

# Combine plots
grid.arrange(p1, p2, p3, p4, ncol = 2, nrow = 2)
Comprehensive PCA analysis of the iris dataset

Comprehensive PCA analysis of the iris dataset

Explained Variance and Choosing Dimensions

One of the most critical decisions in PCA is determining how many principal components to retain. This section explores various methods and criteria for making this choice.

Understanding Explained Variance

The explained variance tells us how much of the original data’s variability is captured by each principal component. The proportion of variance explained by the \(k\)-th component is:

\[\text{Prop. Var}_k = \frac{\lambda_k}{\sum_{i=1}^p \lambda_i}\]

where \(\lambda_k\) is the \(k\)-th eigenvalue and \(p\) is the total number of variables.

Methods for Choosing the Number of Components

1. Kaiser Criterion (Eigenvalue > 1)

Keep components with eigenvalues greater than 1 (for standardized data).

2. Proportion of Variance Explained

Common thresholds: - 80% for exploratory analysis - 90% for more comprehensive analysis
- 95% for minimal information loss

3. Scree Plot Analysis

Look for the “elbow” where the curve levels off.

Let’s implement and compare these methods:

# Function to analyze dimension selection for any PCA result
analyze_dimensions <- function(pca_result, dataset_name) {
  
  # Calculate variance explained
  var_explained <- (pca_result$sdev^2) / sum(pca_result$sdev^2)
  eigenvalues <- pca_result$sdev^2
  cumvar <- cumsum(var_explained)
  
  # Create results data frame
  results <- data.frame(
    Component = 1:length(var_explained),
    Eigenvalue = eigenvalues,
    Prop_Var = var_explained,
    Cumulative_Var = cumvar
  )
  
  # Kaiser criterion (eigenvalue > 1)
  kaiser_components <- sum(eigenvalues > 1)
  
  # Variance thresholds
  comp_80 <- which(cumvar >= 0.80)[1]
  comp_90 <- which(cumvar >= 0.90)[1]  
  comp_95 <- which(cumvar >= 0.95)[1]
  
  cat(paste("Analysis for", dataset_name, ":\n"))
  cat(paste("Kaiser criterion (eigenvalue > 1):", kaiser_components, "components\n"))
  cat(paste("80% variance explained:", comp_80, "components\n"))
  cat(paste("90% variance explained:", comp_90, "components\n"))
  cat(paste("95% variance explained:", comp_95, "components\n\n"))
  
  return(results)
}

# Analyze the iris dataset
iris_results <- analyze_dimensions(pca_iris, "Iris Dataset")
## Analysis for Iris Dataset :
## Kaiser criterion (eigenvalue > 1): 1 components
## 80% variance explained: 2 components
## 90% variance explained: 2 components
## 95% variance explained: 2 components
# Create comprehensive visualization for iris
par(mfrow = c(2, 2))

# Scree plot
plot(iris_results$Component, iris_results$Prop_Var, type = "b", pch = 19,
     main = "Iris - Scree Plot", xlab = "Component", ylab = "Proportion of Variance")
abline(h = 1/ncol(iris_numeric), col = "red", lty = 2)

# Cumulative variance
plot(iris_results$Component, iris_results$Cumulative_Var, type = "b", pch = 19,
     main = "Iris - Cumulative Variance", xlab = "Component", ylab = "Cumulative Variance")
abline(h = c(0.8, 0.9, 0.95), col = c("orange", "blue", "red"), lty = 2)
legend("bottomright", c("80%", "90%", "95%"), col = c("orange", "blue", "red"), lty = 2, cex = 0.8)

# Eigenvalues
plot(iris_results$Component, iris_results$Eigenvalue, type = "b", pch = 19,
     main = "Iris - Eigenvalues", xlab = "Component", ylab = "Eigenvalue")
abline(h = 1, col = "red", lty = 2)
text(2, 1.1, "Kaiser criterion", col = "red")

# Loadings heatmap
loadings_matrix <- pca_iris$rotation[, 1:3]
image(1:3, 1:4, t(loadings_matrix), 
      main = "Loadings Heatmap", 
      xlab = "Principal Component", 
      ylab = "Variable",
      axes = FALSE)
axis(1, at = 1:3, labels = paste0("PC", 1:3))
axis(2, at = 1:4, labels = rownames(loadings_matrix))
Multiple approaches for selecting the optimal number of principal components

Multiple approaches for selecting the optimal number of principal components 

par(mfrow = c(1, 1))

Practical Implementation in R

Manual PCA Implementation

To deepen understanding, let’s implement PCA manually:

# Manual PCA function
manual_pca <- function(data, scale = TRUE) {
  # Step 1: Center (and optionally scale) the data
  if(scale) {
    data_processed <- scale(data, center = TRUE, scale = TRUE)
  } else {
    data_processed <- scale(data, center = TRUE, scale = FALSE)
  }
  
  # Step 2: Calculate covariance matrix
  cov_matrix <- cov(data_processed)
  
  # Step 3: Calculate eigenvalues and eigenvectors
  eigen_result <- eigen(cov_matrix)
  eigenvalues <- eigen_result$values
  eigenvectors <- eigen_result$vectors
  
  # Step 4: Sort by eigenvalues (descending)
  order_indices <- order(eigenvalues, decreasing = TRUE)
  eigenvalues <- eigenvalues[order_indices]
  eigenvectors <- eigenvectors[, order_indices]
  
  # Step 5: Calculate principal components (scores)
  pc_scores <- data_processed %*% eigenvectors
  
  # Return results
  list(
    eigenvalues = eigenvalues,
    eigenvectors = eigenvectors, 
    scores = pc_scores,
    variance_explained = eigenvalues / sum(eigenvalues)
  )
}

# Test manual implementation on iris data
manual_result <- manual_pca(iris_numeric, scale = TRUE)

# Compare with prcomp results
cat("Comparison of eigenvalues:\n")
## Comparison of eigenvalues:
comparison <- data.frame(
  Manual = round(manual_result$eigenvalues, 4),
  prcomp = round(pca_iris$sdev^2, 4)
)
print(comparison)
##   Manual prcomp
## 1 2.9185 2.9185
## 2 0.9140 0.9140
## 3 0.1468 0.1468
## 4 0.0207 0.0207
# The results should be identical (or very close due to numerical precision)
all.equal(manual_result$eigenvalues, pca_iris$sdev^2, tolerance = 1e-10)
## [1] TRUE

Conclusion and Summary

What We’ve Learned

Throughout this tutorial, we have explored Principal Component Analysis from multiple perspectives:

Mathematical Understanding

  • PCA finds the directions of maximum variance through eigendecomposition of the covariance matrix
  • Principal components are orthogonal linear combinations of original variables
  • The transformation preserves distances while reducing dimensionality

Practical Implementation

  • Scaling matters: Always consider whether to scale your data based on variable units and ranges
  • R implementation: prcomp() is the preferred function for most applications
  • Interpretation: Loadings show how original variables contribute to each component

Real-World Applications

We successfully applied PCA to: - Synthetic data: Demonstrated the mathematical concepts with controlled examples - Iris dataset: Showed excellent species separation in reduced dimensions

Dimension Selection

We explored multiple criteria for choosing the number of components: - Kaiser criterion (eigenvalues > 1) - Variance explained thresholds (80%, 90%, 95%) - Scree plot elbow method

Key Takeaways

  1. PCA is a powerful tool for dimensionality reduction, visualization, and exploratory data analysis

  2. Preprocessing is crucial: Always center your data, and scale when variables have different units or vastly different variances

  3. No single “correct” number of components: The choice depends on your specific goals and acceptable information loss

  4. Interpretation requires domain knowledge: Statistical results must be interpreted in the context of your specific field and research questions

  5. PCA assumptions: Works best when relationships between variables are linear and when you want to maximize variance

When to Use PCA

Good applications: - High-dimensional data visualization - Preprocessing for machine learning algorithms - Noise reduction in data - Identifying patterns in complex datasets - Data compression

Consider alternatives when: - Relationships are highly non-linear (consider kernel PCA or other manifold learning techniques) - Interpretability of original features is crucial - Small sample sizes relative to number of variables - Categorical variables dominate the dataset

Final Recommendations

  1. Always visualize your data before and after PCA
  2. Validate results using multiple criteria for dimension selection
  3. Document your choices regarding scaling, number of components, and interpretation
  4. Consider the audience when presenting PCA results - focus on insights rather than technical details
  5. Practice with different datasets to build intuition

PCA remains one of the most valuable tools in the data scientist’s toolkit. By understanding both its mathematical foundations and practical implementation, you can effectively apply it to reveal hidden patterns and structures in high-dimensional data.

References

The mathematical foundations and applications discussed in this tutorial draw from several key sources in statistical learning and multivariate analysis (James et al. 2013; Hastie, Tibshirani, and Friedman 2009; Pearson 1901; Hotelling 1933; Anderson 2003). For deeper exploration of these topics, consult the references below.

This tutorial was created as an educational resource. All code examples are provided under open-source principles and can be freely used and modified for learning purposes.

Anderson, Theodore W. 2003. An Introduction to Multivariate Statistical Analysis. John Wiley & Sons.
Hastie, Trevor, Robert Tibshirani, and Jerome Friedman. 2009. The Elements of Statistical Learning: Data Mining, Inference, and Prediction. Springer Science & Business Media.
Hotelling, Harold. 1933. “Analysis of a Complex of Statistical Variables into Principal Components.” Journal of Educational Psychology 24 (6): 417.
James, Gareth, Daniela Witten, Trevor Hastie, and Robert Tibshirani. 2013. An Introduction to Statistical Learning. Vol. 112. Springer. https://www.statlearning.com/.
Pearson, Karl. 1901. “LIII. On Lines and Planes of Closest Fit to Systems of Points in Space.” The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 2 (11): 559–72.