如何在R语言中进行最大协方差分析（MCA）并绘制系数图？

我之前也碰到过这个问题——R生态里确实没有专门封装MCA的现成包，但其实MCA的核心逻辑就是交叉协方差矩阵的奇异值分解（SVD），完全可以手动实现，步骤清晰且不复杂。下面我会一步步带你完成分析，包括数据预处理、MCA计算和结果可视化。

核心原理回顾 #

MCA的核心步骤你已经提到了：

• 对两个时间序列数据集（记为X和Y，均为n×k的矩阵，n是时间点数量，k是变量数）做去均值处理
• 计算交叉协方差矩阵 C = t(X) %*% Y / (n-1)
• 对C执行奇异值分解：svd(C)，得到的左奇异向量U对应X的模态系数，右奇异向量V对应Y的模态系数，奇异值S代表各耦合模态的协方差强度

具体实现步骤（附代码） #

1. 模拟/准备数据 #

先模拟两个对齐的时间序列作为示例，假设我们有100个时间点，X和Y都是单变量序列：

# 设置随机种子保证可复现
set.seed(123)
n <- 100  # 时间点数量
# 模拟X序列：带趋势的噪声
X <- ts(seq(0, 5, length.out = n) + rnorm(n, 0, 0.5))
# 模拟Y序列：与X相关的滞后序列加噪声
library(data.table) # 用到shift函数，没装的话先运行install.packages("data.table")
Y <- ts(shift(X, n = 3) + rnorm(n, 0, 0.3))

2. 数据预处理（去均值） #

MCA需要处理距平数据，所以先对两个序列去均值：

X_anom <- scale(X, center = TRUE, scale = FALSE)
Y_anom <- scale(Y, center = TRUE, scale = FALSE)

3. 计算交叉协方差矩阵并执行SVD #

# 计算交叉协方差矩阵
C <- t(X_anom) %*% Y_anom / (n - 1)
# 执行奇异值分解
mca_svd <- svd(C)
# 提取核心结果
U <- mca_svd$u  # X的模态系数
V <- mca_svd$v  # Y的模态系数
S <- mca_svd$d  # 奇异值（协方差强度）

4. 计算MCA时间系数（投影后的时间序列） #

时间系数代表每个模态在时间上的变化，是原序列投影到模态向量上的结果：

PCx <- X_anom %*% U  # X序列的MCA时间系数
PCy <- Y_anom %*% V  # Y序列的MCA时间系数

绘制MCA系数图 #

1. 模态时间系数耦合关系图 #

这是最常用的MCA可视化方式，用来展示两个序列第一模态的耦合程度：

# 用ggplot2绘制（没装的话先运行install.packages("ggplot2")）
library(ggplot2)
df <- data.frame(PCx = PCx[,1], PCy = PCy[,1], time = 1:n)

ggplot(df, aes(x = PCx, y = PCy)) +
  geom_point(color = "#2c3e50", size = 2) +
  geom_smooth(method = "lm", color = "#e74c3c", se = FALSE, lwd = 1.2) +
  labs(x = "X的MCA时间系数（第一模态）", 
       y = "Y的MCA时间系数（第一模态）",
       title = "MCA第一模态时间系数耦合关系") +
  theme_minimal()

2. 原序列与MCA时间系数对比图 #

可以直观看到模态系数和原序列的对应关系：

# 用基础绘图分栏展示
par(mfrow = c(2,1))
plot(X, main = "原X序列 vs MCA时间系数", col = "#3498db", lwd = 1.5)
lines(PCx[,1], col = "#e74c3c", lwd = 2)
legend("topleft", legend = c("原序列", "MCA时间系数"), 
       col = c("#3498db", "#e74c3c"), lwd = 2, bty = "n")

plot(Y, main = "原Y序列 vs MCA时间系数", col = "#3498db", lwd = 1.5)
lines(PCy[,1], col = "#e74c3c", lwd = 2)
legend("topleft", legend = c("原序列", "MCA时间系数"), 
       col = c("#3498db", "#e74c3c"), lwd = 2, bty = "n")
par(mfrow = c(1,1))

注意事项 #

• 如果你的时间序列是多变量的（比如X是n×p，Y是n×q），代码逻辑完全一致，只是U会是p×r，V是q×r（r是模态数，取min(p,q)）
• 可以通过S/sum(S)计算各模态的协方差贡献占比，筛选出主要模态进行分析
• 去均值是关键步骤，否则协方差矩阵会包含均值的干扰，导致结果偏差

内容的提问来源于stack exchange，提问作者Yang Yang