R语言源码那些事（二） PCA分析

R语言PCA主成分分析全流程及源码追溯：factoextra vs 基础方法对比解析

一、PCA分析流程回顾

library(ggplot2)
library(factoextra)
library(patchwork)
 
data(iris)
df <- iris[, 1:4]
 
# 标准化数据并进行 PCA
pca_res <- prcomp(df, scale. = TRUE)
 
# 可视化：主成分方差解释比例
fviz_eig(pca_res)
 
# 可视化：样本个体图，按物种着色
p1 <- fviz_pca_ind(pca_res,
                  geom.ind = "point",
                  col.ind = iris$Species,
                  palette = "npg",
                  addEllipses = TRUE,
                  legend.title = "Species")
 
# 可视化：变量贡献度图
p2 <- fviz_pca_var(pca_res,
                  col.var = "contrib",
                  gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
                  repel = TRUE)
 
# 可视化：双标图，综合个体与变量信息
p3 <- fviz_pca_biplot(pca_res,
                     col.ind = iris$Species,
                     palette = "jco",
                     col.var = "black",
                     addEllipses = TRUE,
                     legend.title = "Species")
 
# 三图组合展示
p1 + p2 + p3

---

二、factoextra 包对变量贡献度计算的源码追溯

factoextra::get_pca_var() 函数其实基于 prcomp 的输出，先计算变量与主成分的相关系数矩阵，然后基于相关性的平方（cos²）计算贡献度，核心代码如下：

if (inherits(res.pca, "prcomp")) {
    # 计算变量载荷乘以主成分标准差，得到相关系数矩阵
    var_cor_func <- function(var.loadings, comp.sdev) {
        var.loadings * comp.sdev
    }
    var.cor <- t(apply(res.pca$rotation, 1, var_cor_func, res.pca$sdev))
    
    # 计算贡献度、相关性等指标
    var <- .get_pca_var_results(var.cor)
}

其中，计算贡献度的函数定义为：

.get_pca_var_results <- function(var.coord) {
    var.cor <- var.coord
    var.cos2 <- var.cor^2                  # 相关性平方，即cos²
    comp.cos2 <- apply(var.cos2, 2, sum)  # 每个主成分的总cos²
    
    contrib <- function(var.cos2, comp.cos2) {
        var.cos2 * 100 / comp.cos2         # 贡献度百分比
    }
    var.contrib <- t(apply(var.cos2, 1, contrib, comp.cos2))
    
    # 命名列名
    colnames(var.coord) <- colnames(var.cor) <- colnames(var.cos2) <- colnames(var.contrib) <- paste0("Dim.", 1:ncol(var.coord))
    
    list(coord = var.coord, cor = var.cor, cos2 = var.cos2, contrib = var.contrib)
}

---

三、提取主成分得分和载荷

为方便对比，我们用基础 prcomp 输出，按上述方法手动计算变量的相关系数和贡献度：

#1 提取变量
var_res <- factoextra::get_pca_var(pca_res)
	#载荷
		var_res$coord 
		var.cor <- pca_res$rotation %*% diag(pca_res$sdev) %>% as.data.frame() #base提取法
	#变量对主成分的贡献
		var_res$contrib
	#变量在主成分方向带相关性平方
		var_res$cos2
 
#2 提取样本相关信息
	#样本主成分得分
		smp_res <- factoextra::get_pca_ind(pca_res)
		smp_res$coord 
		pca_res$x #base提取法
	# Contributions to the PCs
		smp_res$contrib
	# Quality of representation
		smp_res$cos2 
 
 
#3 主成分的贡献
	summary_res <- summary(pca_res)
	summary_res$importance  ##base提取法
	factoextra::get_eigenvalue(pca_res)