lvzhaopeng的数字花园

统计学高级 CA 代码比较

23分钟阅读

1️⃣ 计算逻辑的相同点和区别

相同点：基本 CCA 步骤

步骤上一个代码 (geigen + compu)rcc()
中心化scale(X, center=TRUE, scale=FALSE)X.aux = scale(X, center=TRUE, scale=FALSE)
协方差矩阵Cxx = var(X) + λICxx = var(X) + λI（ridge）或者 cov.shrink()（shrinkage）
典型相关系数 & 权向量res <- geigen(Cxy, Cxx, Cyy) → xcoef/ycoef使用 Cholesky 白化 + SVD 计算 xcoef/ycoef
典型变量xscores = X.aux %*% xcoefU = X.aux %*% xcoef
载荷cor(X, xscores)loadings = xcoef / ycoef（权向量本身作为载荷，也可算 correlation 载荷）

✅ 基本数学思想完全相同：

  1. 中心化 X、Y
  2. 计算协方差矩阵（带正则化）
  3. 求广义特征值或通过 SVD 求典型向量
  4. 得到典型变量（U、V）
  5. 计算载荷（或直接使用权向量）

您提出了一个非常关键且容易混淆的点。您的直觉是正确的，cor(X, xscores) 的结果不是一个单一的标量值，而是一个矩阵。这个矩阵包含了丰富的信息，它精确地反映了原始变量集 X 与提取出的新综合变量 xscores 之间的全方位相关关系。

下面我们来详细解释它如何反应，以及我们应该如何解读它。

1. 维度和结构

首先，明确各个元素的维度：

  • X：是一个 n x p 的矩阵（n 个样本，p 个原始变量）。
  • xscores：是一个 n x m 的矩阵（n 个样本，m 个主成分得分）。通常 m <= p。
  • cor(X, xscores)：计算这两个矩阵的列与列之间的相关系数。
    • 结果是一个 p x m 的矩阵。
    • 行 代表原始的 p 个变量。
    • 列 代表提取出的 m 个主成分（xscores）。

我们把这个矩阵称为 因子载荷矩阵，通常记为 Loadings。

2. 如何解读这个矩阵：它反映了什么？

这个 p x m 的载荷矩阵是主成分分析的核心解读工具，它从两个维度回答了您的问题“怎么能反应”：

视角一：按行看（观察原始变量）

每一行描述了 一个原始变量 与 所有主成分 的相关关系。

  • 目的：理解每个原始变量主要由哪些主成分所代表的新维度来解释。
  • 解读：
    • 如果一个变量在某个主成分上的载荷（相关系数）的绝对值很大（例如 > |0.5| 或 > |0.7|），说明这个变量高度相关于该主成分。
    • 通过观察一行中哪个载荷最大，可以判断这个变量最主要由哪个主成分所代表。
    • 举例：假设我们有变量：身高、体重、血压、心率。
      • 我们可能发现“身高”和“体重”在第一主成分上有很高的正载荷（如0.9和0.85）。
      • 这告诉我们，第一主成分主要捕捉的是“体型大小”这个综合信息。

视角二：按列看（观察主成分）

每一列描述了 一个主成分 与 所有原始变量 的相关关系。

  • 目的：为每个主成分（这个新的、抽象的综合指标）赋予实际的业务含义。
  • 解读：
    • 观察一列中哪些原始变量的载荷（相关系数）的绝对值很大。
    • 这些高载荷的变量共同定义了该主成分的实际意义。
    • 举例：承接上面的例子。
      • 第一主成分这一列，在“身高”、“体重”上载荷很高，而在“血压”、“心率”上载荷很低。我们就可以将第一主成分解释为“体型因子”。
      • 第二主成分可能发现在“血压”、“心率”上载荷很高，而在“身高”、“体重”上载荷很低。我们就可以将第二主成分解释为“心血管因子”。

3. 一个具体的数字例子

假设 X 有3个变量 (X1, X2, X3)，我们提取了2个主成分 (PC1, PC2)。cor(X, xscores) 的结果可能如下：

变量与 PC1 的相关系数与 PC2 的相关系数
X10.950.05
X20.870.35
X3-0.100.98

解读：

  • 从主成分看（按列）：

    • PC1：与 X1 和 X2 高度正相关。PC1 主要代表了 X1 和 X2 所共同蕴含的信息。我们可以将 PC1 命名为 “X1-X2 因子”。
    • PC2：与 X3 高度正相关，与 X1 和 X2 几乎无关。PC2 主要代表了 X3 的独立信息。我们可以将 PC2 命名为 “X3 因子”。

  • 从变量看（按行）：

    • X1：其信息几乎完全被 PC1 所捕获（相关系数 0.95）。
    • X2：其信息主要被 PC1 所捕获（相关系数 0.87），但也有一小部分体现在 PC2 上。
    • X3：其信息几乎完全独立，由 PC2 所代表（相关系数 0.98）。

4. 可视化工具：载荷图

为了更直观地“反应”这种关系，我们常常会绘制 载荷图。

  • 方法：以 PC1 为横轴，PC2 为纵轴，将每个变量作为一个点画在图中，点的坐标就是该变量在 (PC1, PC2) 上的载荷。
  • 解读：
    • 点离原点越远，说明它与该主成分的关系越强。
    • 在同一个方向上聚集的点，代表它们被同一个主成分所代表，可能具有相似的含义。
    • 点与坐标轴夹角的余弦值近似等于其与对应主成分的相关系数。

结论

所以，cor(X, xscores) 这个 p x m 的矩阵，绝不是一个简单的值。它是连接原始数据空间和降维后主成分空间的桥梁和罗塞塔石碑。

  • 它通过列视角告诉我们每个主成分的实际意义（为抽象主成分命名）。
  • 它通过行视角告诉我们每个原始变量的信息被分配到了哪些主成分上。

正是通过深入分析这个相关系数矩阵，我们才能理解主成分分析的结果，并赋予其可解释的现实意义。

您完全正确！您的理解是精确的，这正是计算相关矩阵的数学过程。

cor(X, xscores) 最终得到一个 p x m 的矩阵，是通过系统地、成对地计算 X 的每一列与 xscores 的每一列之间的相关系数来实现的。

让我们一步步拆解这个过程。

1. 输入数据的维度

  • X: 维度为 n x p

    • n 是观测样本数
    • p 是原始变量的个数
    • 可以表示为：X = [X₁, X₂, ..., Xₚ]，其中每个 Xⱼ 都是一个 n x 1 的列向量。

  • xscores: 维度为 n x m

    • n 是观测样本数（与 X 相同）
    • m 是主成分的个数（通常 m <= p）
    • 可以表示为：xscores = [PC₁, PC₂, ..., PCₘ]，其中每个 PCₖ 都是一个 n x 1 的列向量。

2. 计算过程：您的理解完全正确

对于 X 的每一列（每一个原始变量 Xⱼ）和 xscores 的每一列（每一个主成分 PCₖ），计算它们之间的皮尔逊相关系数 r_{jk}。

这个系数的计算公式，正如您所说，是：
r_{jk} = cor(Xⱼ, PCₖ) = cov(Xⱼ, PCₖ) / (std(Xⱼ) * std(PCₖ))

其中：

  • cov(Xⱼ, PCₖ)：是 Xⱼ 和 PCₖ 两个向量之间的协方差（一个标量）。
  • std(Xⱼ)：是变量 Xⱼ 的标准差（一个标量）。
  • std(PCₖ)：是主成分 PCₖ 的标准差（一个标量）。

因此，每一个 r_{jk} 都是一个标量。

3. 组装成 p x m 矩阵

现在，我们系统地执行这个计算：

  1. 固定主成分 PCₖ（即固定 xscores 的列 k）。

  2. 让 k=1（第一个主成分），然后遍历 X 的所有 p 个变量：

    • 计算 r₁₁ = cor(X₁, PC₁)
    • 计算 r₂₁ = cor(X₂, PC₁)
    • …
    • 计算 rₚ₁ = cor(Xₚ, PC₁)
      这 p 个计算结果形成了一个 p x 1 的列向量。这个向量代表了所有原始变量与第一个主成分的相关性。

  3. 移动到下一个主成分 PC₂，重复步骤2：

    • 计算 r₁₂ = cor(X₁, PC₂)
    • 计算 r₂₂ = cor(X₂, PC₂)
    • …
    • 计算 rₚ₂ = cor(Xₚ, PC₂)
      这形成了第二个 p x 1 的列向量。

  4. 对所有的 m 个主成分重复此过程。

  5. 最后，将得到的 m 个 p x 1 的列向量并排组合在一起，就形成了一个 p x m 的矩阵。

最终矩阵的结构：

[
\text{cor}(X, \text{xscores}) =
\begin{bmatrix}
cor(X_1, PC_1) & cor(X_1, PC_2) & \cdots & cor(X_1, PC_m) \
cor(X_2, PC_1) & cor(X_2, PC_2) & \cdots & cor(X_2, PC_m) \
\vdots & \vdots & \ddots & \vdots \
cor(X_p, PC_1) & cor(X_p, PC_2) & \cdots & cor(X_p, PC_m) \
\end{bmatrix}
]

4. 一个简单的数值模拟（用于验证）

假设我们有：

  • n = 5 个样本
  • p = 2 个原始变量：X1, X2
  • m = 2 个主成分得分：PC1, PC2

在 R 语言中：

# 假设的小数据
X <- matrix(c(1,2,3,4,5, 2,4,6,8,10), ncol=2)
colnames(X) <- c("X1", "X2")
xscores <- matrix(c(1,2,3,4,5, -1,1,-1,1,-1), ncol=2)
colnames(xscores) <- c("PC1", "PC2")
 
# 方法1：使用 cor() 函数直接得到 p x m 矩阵
result_matrix <- cor(X, xscores)
print(result_matrix)
# 输出将会是一个 2x2 的矩阵：
#           PC1        PC2
# X1  1.0000000  0.0000000
# X2  1.0000000  0.0000000
 
# 方法2：手动计算来验证，正如您所描述的
# 计算 X1 和 PC1 的相关系数 (对应结果矩阵的 [1,1] 位置)
cov_x1_pc1 <- cov(X[,1], xscores[,1])
std_x1 <- sd(X[,1])
std_pc1 <- sd(xscores[,1])
manual_r_11 <- cov_x1_pc1 / (std_x1 * std_pc1)
print(manual_r_11) # 会输出 1
 
# 您可以继续手动计算其他位置的值，会发现它们与方法1的 result_matrix 完全一致。

区别：

特点上一个代码rcc()
正则化可选 ridge（λI）ridge 或 shrinkage（Ledoit-Wolf 自适应收缩）
SVD 实现没有考虑 X/Y 行列数差异根据 p >= q 选择对 mat 或 t(mat) 做 SVD
输出典型变量得分 + 载荷典型相关系数、载荷、典型变量、解释方差比例、lambda 值、样本/变量名
缺失值NA 填 0，pairwiseshrinkage 模式下不允许 NA
统计检验没有也没有（MixOmics 只关注数值计算）
方差解释没有有 prop_expl_var（典型变量解释原始变量方差比例）

2️⃣ 为什么 mixOmics 中不做显著性检验

包数据类型假设背景是否做检验原因
vegan::cca()环境因子 X → 群落丰度 Y假设因果关系（环境因素影响微生物群落）有 permutation test可以检验约束效应是否显著，验证因果假设
mixOmics::rcc()高通量分子组数据（代谢物/转录本/蛋白）关注相关性/多变量关联无检验不假设因果关系，关注模式发现和潜变量相关性，p-value 意义不大

  1. rcc() 和 geigen+compu 核心算法一致：协方差矩阵 → 广义特征值/SVD → 典型变量 → 载荷。

  2. rcc() 更现代：

    • 支持 ridge/shrinkage 正则化
    • 自动选择 SVD 路径
    • 计算解释方差比例

  3. 不做统计检验的原因：

    • 高维分子数据通常只关注相关性模式
    • 没有明确因果假设
    • 检验 p-value 在高维情况下容易失真

  4. vegan 做检验是因为有环境→微生物的因果假设，可用 permutation test 验证。

下面把你前面提到的三类检验——整体模型检验、逐成分检验、特征层面检验——在“未通过（不显著）”时各自意味着什么，以及该怎么处理，一次说清楚。

1) 整体模型检验（整体典型相关是否显著）没通过

含义：

  • 在你当前的数据尺度/预处理/样本量下，转录组与代谢组之间没有显著的全局线性关联（至少不是强到超出随机噪声）。
  • 也可能是“有信号，但被噪声、批次、协变量或非线性关系掩盖”。

常见原因：

  • 样本量偏小、测序噪声大、批次效应/技术偏差未校正。
  • 变量数远大于样本数导致不稳定；缩放/归一化不当。
  • 真实关系是非线性或条件性（仅在某亚群或某条件成立）。
  • 配对不严（样本错配或时间不同步）。

下一步建议：

  • 强化预处理：批次校正（ComBat 等）、对数/方差稳定化、零膨胀/极端值处理。
  • 控制协变量（年龄、性别、批次、分组）或分层分析。
  • 提升稳健性：正则化/稀疏 CCA（sCCA）、核 CCA/非线性方法。
  • 聚合到通路/模块层面（GSVA、WGCNA 模块特征）再做 CCA，提高信噪比。
  • 增加样本量或在外部数据复现。

2) 逐成分检验（第 k 对典型变量是否显著）没通过

含义：

  • 第 对典型变量不带来新的稳定信息；在前面成分解释掉可检出的共享结构后，余下成分更像噪声或过拟合。
  • 常见现象：第一成分显著，后续成分不显著。

常见原因：

  • 设定的成分数过多；高维小样本导致后续成分不稳定。
  • 变量间共线性强，后续方向可解释方差极少。
  • 正则化不足/过强，导致方向要么噪声化，要么被抑制。

下一步建议：

  • 用置换/CV 决定截断的成分数（只保留显著成分）。
  • 采用稀疏 CCA（mixOmics 的 rgcca/spls/sPLS-CCA）减少噪声变量干扰。
  • 若关注解释性，宁可少而稳（1–2 个显著成分），并报告剩余成分不显著。
  • 若怀疑非线性，考虑核方法或先做非线性降维再 CCA。

3) 特征层面检验（单个基因/代谢物与典型变量的相关）没通过

含义：

  • 没有单个基因或代谢物与典型变量表现出稳定且可重复的强相关；信号可能弥散在很多弱效应上。
  • 另一种可能：多重检验校正后功效不足，个体特征达不到阈值。

常见原因：

  • 多重比较导致阈值严格，单个效应难以显著。
  • 特征不稳定（高相关、共线性）→ 单个变量 p 值摇摆。
  • 典型变量更多反映通路/模块级变化，而非个别特征。

下一步建议：

  • 从集合/通路层面解读：对按载荷排序的基因/代谢物做 KEGG/GO/Reactome 富集；即便单基因不显著，通路可能显著。
  • 做稳定性选择/自助法（bootstrap）评估载荷的重现性，汇报稳定特征而非一次性 p 值。
  • 使用稀疏 CCA（sCCA）以得到更可解释的稀疏载荷；或对高相关特征先做模块聚合。
  • 结合先验知识（代谢通路、调控网络）进行目标式验证，而不是纯数据驱动逐个特征显著性。

小结成一句话

  • 整体不显著：两组学之间缺乏稳健的全局线性共变结构（或被噪声/非线性/协变量遮蔽）。
  • 某成分不显著：该方向没有新信息，应减少成分或改用稀疏/非线性方法。
  • 特征不显著：信号更可能在通路/模块层面而非单个分子；用稀疏化、稳定性选择和富集分析提升可解释性。

如果你愿意，我可以把你现在的 CCA 结果跑一版完整报告模版（整体置换检验 + 成分显著性 + 特征稳定性 + 通路富集的图表清单），你只要替换数据矩阵就能直接产出可投报告。

CCA代码

注意如果要代码使用claude4.0 可以解释 1. 模型是否可靠 2. 解释比例 3. 相关性大小，这张图里面出现的就是最后关心的

library(mixOmics)
library(tidyverse)
 
# 计算变量稳定性函数
# =============================
bootstrap_stability <- function(X, Y, method = "ridge", lambda1 = 0.5, lambda2 = 0.05, nboot = 100) {
	nvar_X <- ncol(X)
	nvar_Y <- ncol(Y)
	loading_X_list <- matrix(0, nrow = nvar_X, ncol = nboot)
	loading_Y_list <- matrix(0, nrow = nvar_Y, ncol = nboot)
	
	for (b in 1:nboot) {
		idx <- sample(1:nrow(X), replace = TRUE)
		res_boot <- rcc(Y[idx, ], X[idx, ], method = method, lambda1 = lambda1, lambda2 = lambda2)
		loading_X_list[, b] <- res_boot$loadings$X[, 1] # 第一成分
		loading_Y_list[, b] <- res_boot$loadings$Y[, 1]
	}
	
	stability_X <- apply(loading_X_list, 1, function(x) cor(x, rep(mean(x), length(x))))
	stability_Y <- apply(loading_Y_list, 1, function(x) cor(x, rep(mean(x), length(x))))
	
	list(stability_X = stability_X, stability_Y = stability_Y)
}
# =============================
# 主函数
# =============================
cca_analysis <- function(X, Y, method = "ridge", lambda1 = 0.5, lambda2 = 0.05, nperm = 200, ncomp = 3, nboot = 100) {
	
	# 1. 原始 CCA
	res <- rcc(Y, X, method = method, lambda1 = lambda1, lambda2 = lambda2)
	# 动态确定可用成分数
	ncomp_real <- min(ncomp, ncol(res$variates$X), ncol(res$variates$Y))
	U <- res$variates$X[, 1:ncomp_real, drop = FALSE]
	V <- res$variates$Y[, 1:ncomp_real, drop = FALSE]
	
	
	# 2. 置换检验
	perm_cor <- replicate(nperm, {
		Yperm <- Y[sample(nrow(Y)), ]
		rcc(Yperm, X, method = method, lambda1 = lambda1, lambda2 = lambda2)$cor[1]
	})
	pval_global <- mean(perm_cor >= res$cor[1])
	
	# 各成分 p 值
	component_pvals <- map_dbl(1:length(res$cor), function(k) {
		perm_cork <- replicate(nperm, {
			Yperm <- Y[sample(nrow(Y)), ]
			rcc(Yperm, X, method = method, lambda1 = lambda1, lambda2 = lambda2)$cor[k]
		})
		mean(perm_cork >= res$cor[k])
	})
	
	variance_explained <- tibble(
		Component = 1:ncomp_real,
		VarExpl_X = res$prop_expl_var$X[1:ncomp_real],  # X 的解释率
		VarExpl_Y = res$prop_expl_var$Y[1:ncomp_real]   # Y 的解释率
	)
	# 4. 可视化
 
	
	# 6. 返回结果
	list(
		model = res,
		global_pval = pval_global,
		component_pvals = tibble(Component = 1:length(res$cor), CanonicalCor = res$cor, Pvalue = component_pvals),
		variance_explained = variance_explained
	)
}
# =============================
# 测试
# =============================
data(nutrimouse)
X <- nutrimouse$lipid
Y <- nutrimouse$gene
res_test <- cca_analysis(X, Y, method = "ridge", lambda1 = 1, lambda2 = 1, nperm = 200, ncomp = 3, nboot = 50)
res_test$global_pval
res_test$component_pvals
 
 
res_test$model$variates$X #典型变量
 
#等价
res_test$variance_explained
res_test$model$prop_expl_var
plotIndiv(res_test$model)
plotVar(res_test$model,style = "ggplot2")

载荷图代码 需要优化

###基础载荷图
 
loadings_X <- res_test$model$loadings$X[, 1:2]
 
loadings_Y <- res_test$model$loadings$Y[, 1:2]
 
# 合并数据
 
df_loadings <- rbind(
 
	data.frame(Feature = rownames(loadings_X), loadings_X, Dataset = "X"),
 
	data.frame(Feature = rownames(loadings_Y), loadings_Y, Dataset = "Y")
 
)
 
# 缩放因子，避免箭头过长
 
scale_factor <- 0.3
 
df_loadings$X1 <- df_loadings$X1 * scale_factor
 
df_loadings$X2 <- df_loadings$X2 * scale_factor
 
# 绘制 biplot（箭头 + 标签）
 
ggplot() +
 
	geom_segment(data = df_loadings, 
 
				 aes(x = 0, y = 0, xend = X1, yend = X2, color = Dataset),
 
				 arrow = arrow(length = unit(0.2, "cm")), alpha = 0.7) +
 
	geom_text(data = df_loadings, 
 
			  aes(x = X1, y = X2, label = Feature, color = Dataset),
 
			  check_overlap = TRUE, size = 3) +
 
	geom_vline(xintercept = 0, linetype = "dashed") +
 
	geom_hline(yintercept = 0, linetype = "dashed") +
 
	labs(x = "CCA Component 1", y = "CCA Component 2") +
 
	theme_minimal() +
 
	theme(legend.position = "bottom")
 
professional_biplot <- function(rcc_result, stability_results = NULL, 
 
								shrinkage_method = "stability", threshold = 0.1) {
 
	
 
	loadings_X <- rcc_result$loadings$X[, 1:2]
 
	loadings_Y <- rcc_result$loadings$Y[, 1:2]
 
	
 
	# 根据收缩方法处理载荷
 
	if (shrinkage_method == "stability" && !is.null(stability_results)) {
 
		# Type II收缩：基于稳定性
 
		loadings_X <- apply(loadings_X, 2, function(x) 
 
			shrink_by_stability(x, stability_results$stability_X, threshold))
 
		loadings_Y <- apply(loadings_Y, 2, function(x) 
 
			shrink_by_stability(x, stability_results$stability_Y, threshold))
 
	} else if (shrinkage_method == "magnitude") {
 
		# Type I收缩：基于载荷大小
 
		loadings_X <- apply(loadings_X, 2, function(x) x * importance_weights(x, threshold))
 
		loadings_Y <- apply(loadings_Y, 2, function(x) x * importance_weights(x, threshold))
 
	}
 
	
 
	# 创建数据框
 
	df_loadings <- rbind(
 
		data.frame(Feature = rownames(loadings_X), loadings_X, Dataset = "X"),
 
		data.frame(Feature = rownames(loadings_Y), loadings_Y, Dataset = "Y")
 
	)
 
	
 
	# 绘图
 
	ggplot() +
 
		geom_segment(data = df_loadings, 
 
					 aes(x = 0, y = 0, xend = X1, yend = X2, color = Dataset),
 
					 arrow = arrow(length = unit(0.2, "cm")), alpha = 0.7) +
 
		geom_text(data = df_loadings[abs(df_loadings$X1) > 0.05 | abs(df_loadings$X2) > 0.05,], 
 
				  aes(x = X1, y = X2, label = Feature, color = Dataset),
 
				  check_overlap = TRUE, size = 3) +
 
		geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.5) +
 
		geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
 
		labs(x = paste0("CCA Component 1 (", round(rcc_result$cor[1], 3), ")"),
 
			 y = paste0("CCA Component 2 (", round(rcc_result$cor[2], 3), ")"),
 
			 title = paste("Professional CCA Biplot -", shrinkage_method, "shrinkage")) +
 
		theme_minimal() +
 
		theme(legend.position = "bottom")
 
}
 
#专业载荷图
 
professional_biplot(res_test$model)

载荷图直接用这个

plotVar(res_test$model,style = "ggplot2",plot = F)

关系图谱

反向链接