多元因子分析（MFA）及其在R中实现
（Abdi et al, 2013）
MFA已被广泛用于多种领域，例如：
调查分析，对象为被调查的个体，变量为调查问题，问题可按主题划分为变量集；
感官分析，对象是食品，第一组变量包括味觉变量（甜，苦等），第二组包括化学变量（pH值、葡萄糖速率等）；
生态学分析，对象为观测地域，第一组变量描述了环境特性，第二组描述物种数量；
多组学分析，对象是样本，不同的变量集可以为转录组、蛋白组、代谢组等测量数据；
时间序列，可以将各变量集定义为每个变量在不同时间的测量值；
本篇简介R包FactoMineR的MFA方法，以及结合factoextra包实现可视化。分析中将展示如何揭示重要的变量，这些变量对解释数据集的方差起最大的作用。
数据集
FactoMineR包的内置数据集wine，记录了21种葡萄酒的信息。
library(FactoMineR)#数据集，详情 ?wine
data(wine)
head(wine[1:6])
该数据集中，每一行代表一种葡萄酒，共21行。
每一列代表了葡萄酒的属性信息，共31列。其中前2列是分类变量，label（Saumur、Bourgueil或Chinon）代表了葡萄酒的产地，soil（Reference、Env1、Env2或Env4）代表了植物的土壤环境；第3-31列是定量变量，代表了葡萄酒各类品质特征的打分。
整体而言，对于数据集中所有代表葡萄酒属性的变量，可将它们组合为6组变量集，包含1组分类变量和5组定量变量：
第1-2列，共2列，包括葡萄酒产地等属性，该变量集代表葡萄酒“来源属性”（origin）；
第3-7列，共5列，包括摇动前的气味强度、香气质量、果味等属性，该变量集代表葡萄酒“摇动前的气味属性”（odor before shaking）；
第8-10列，共3列，包括视觉上的区分度等属性，该变量集代表葡萄酒“视觉属性”（visual）；
第11-20列，共10列，包括摇动后的气味强度、香气质量、果味等属性，该变量集代表葡萄酒“摇动后的气味属性”（odor after shaking）；
第21-29列，共9列，包括酸度等属性，该变量集代表葡萄酒“味觉属性”（taste）；
第30-31列，共2列，包括质量等属性，该变量集代表葡萄酒“质量属性”（quality）。
接下来期望根据这些组合属性，对21种葡萄酒的品质进行综合评价，以及寻找它们之间的共同特征或差异，即可通过MFA来完成。
R包FactoMineR的MFA
FactoMineR中，MFA通过MFA()函数执行。
#执行 MFA，详情 ?MFA#其中，摇动前后的气味、视觉和味觉作为 active groups，来源和质量作为 supplementary groups
res.mfa <- MFA(wine, group = c(2, 5, 3, 10, 9, 2), type = c('n', 's', 's', 's', 's', 's'),
name.group = c('origin', 'odor_before_shaking', 'visual', 'odor_after_shaking', 'taste', 'quality'),
num.group.sup = c(1, 6), graph = TRUE)
summary(res.mfa)
其中，group指定列的集合作为变量集（数据集的列需要提前按组排列好）。
type指定变量标准化方法，与各变量集一一对应。MFA的目标是整合描述相同观察对象的不同变量集，由于各变量的量纲不一致，为了平衡每组变量的影响（或者说为了使各变量集具有可比性），需对各变量集执行标准化处理。对于定量变量集（将用于PCA），通常通过将每个变量除以其标准差进行标准化（即Z-scores）；对于分类变量集（将用于MCA），每组通过将其所有元素除以称为其第一奇异值（该值是标准差矩阵等效项）的量进行标准化。type参数中，n即用于这种对分类变量的标准化，s用于定量变量的标准化，详情参阅帮助文档。
name.group为变量集指定名称。
结果概要显示，前两轴承载了68.87%的总方差，具有较高的代表性。
接下来对于对象（葡萄酒）和变量（属性）间的关系，以及变量对排序空间的贡献等，可主要通过排序图评估。
graph = TRUE参数，可在计算过程中自动输出图形。不过我们有更便捷的可视化方法，见下文。
同样地，对于MFA主要信息的提取，除了可直接在结果中提取外，也和排序图的展示一样使用其它工具辅助完成，见下文。
#提取主要信息，例如res.mfa
res.mfa$eig #各轴特征值
R包factoextra的可视化方法
接下来展示使用factoextra包来帮助解释和可视化MFA结果。
尽管下述方法并非全部来自factoextra包。
特征值/方差
首先查看MFA的特征值概况。
#借助 factoextra 包更好地提取数据和可视化
library(factoextra)
#提取每个维度（轴）的特征值/承载的方差
eig.val <- get_eigenvalue(res.mfa)
head(eig.val)
#可视化各轴特征值/方差
fviz_screeplot(res.mfa)
eigenvalue为各MFA轴的特征值，即各轴承载的方差；variance.percent为各轴特征值与所有轴特征值总和的比例，即各轴的方差贡献率；cumulative.variance.percent为累积贡献率。
变量集整体间的关系及其对排序空间的贡献
MFA结果中，可首先观测各变量集整体之间的关系。
#变量集结果，用于反映变量集整体对 MFA 排序空间的贡献
group <- get_mfa_var(res.mfa, 'group')
group
#各变量集作为整体，计算了它们与各 MFA 轴的相关性
head(group$correlation)
#排序坐标也可用于反映各变量集与各 MFA 轴的相关程度，前两轴为例展示
head(group$coord)
fviz_mfa_var(res.mfa, 'group', axes = 1:2)
对应MFA命令行参数，摇动前后的气味、视觉和味觉作为active groups，图中红色显示；来源和质量作为supplementary groups，图中绿色显示。
根据变量集在排序图中的位置，可以评估各变量集与MFA轴之间的相关性及相对贡献。第一轴上4个active groups的坐标几乎相同，意味着它们对第一维的贡献相似；对于第二轴，摇动前后的气味具有相对较大的坐标，表明它们对第二维的贡献最大。
#直接显示了各变量集对各排序轴的贡献度head(group$contrib)
#可视化各变量集对前两轴的贡献度
fviz_contrib(res.mfa, 'group', axes = 1)
fviz_contrib(res.mfa, 'group', axes = 2)
#变量集之间的相关性矩阵
res.mfa$group$RV
#相关图
library(corrplot)
corrplot(res.mfa$group$RV, method = 'number', number.cex = 0.8, diag = FALSE, tl.cex = 0.8)
corrplot(res.mfa$group$RV, add = TRUE, type = 'upper', method = 'pie', diag = FALSE, tl.pos = 'n', cl.pos = 'n')
#这种相关性可通过置换检验确定重要性
#例如摇动前的气味（原始数据集的 3-7 行）和摇动后的气味（原始数据集的 11-20 行）之间的相关性的置换检验
#记得执行 Z-scores 标准化变量集
coeffRV(scale(wine[3:7]), scale(wine[11:20]))$p.value
各变量间的关系及其对排序空间的贡献
上述查看了变量集整体间的关系，接下来可以查看各变量间的关系。
#各具体变量对排序空间的贡献，数据结构同上述“变量集整体”quanti.var <- get_mfa_var(res.mfa, 'quanti.var')
quanti.var
#例如同样以排序图可视化各变量与 MFA 轴的关系，前两轴为例展示
head(quanti.var$coord) #排序坐标
fviz_mfa_var(res.mfa, 'quanti.var', axes = 1:2, palette = 'jco', col.var.sup = 'violet', repel = TRUE)
图中颜色代表了变量所属的变量集，定量变量展示为向量。相关圈（correlation circle）代表所有定量变量范数标准化向量。
两个变量之间的夹角为锐角，则二者存在正相关；钝角为负相关；趋于正交则相关性很低。变量与哪个轴的夹角越小，表明其越作用于哪个轴；变量在某轴的投影长度越长，则表明其对该轴的贡献越大。
#直接显示了各变量对各排序轴的贡献度head(quanti.var$contrib)
#展示 top20 变量对前两轴的贡献度
fviz_contrib(res.mfa, choice = 'quanti.var', axes = 1, top = 20, palette = 'jco')
fviz_contrib(res.mfa, choice = 'quanti.var', axes = 2, top = 20, palette = 'jco')
红色虚线表示平均贡献值，便于观测哪些变量大于平均贡献。
#可选使用渐变色，在排序图中按变量贡献对 active groups 变量上色#这里方便观测数据，将所有变量以点展示
fviz_mfa_var(res.mfa, 'quanti.var', axes = 1:2, col.var = 'contrib', geom = c('point', 'text'),
gradient.cols = c('#00AFBB', '#E7B800', '#FC4E07'), col.var.sup = 'violet', repel = TRUE)
select_val <- dimdesc(res.mfa, axes = 1:2, proba = 0.001)
summary(select_val)
#保留最具典型的一些变量的排序图
varsig <- res.mfa$quanti.var$cor
varsig <- subset(varsig, rownames(varsig) %in% unique(c(rownames(select_val$Dim.1$quanti), rownames(select_val$Dim.2$quanti))))
plot(varsig[ ,1:2], asp = 1, type = 'n', xlim = c(-1, 1), ylim = c(-1, 1))
abline(h = 0, lty = 3)
abline(v = 0, lty = 3)
symbols(0, 0, circles = 1, inches = FALSE, add = TRUE)
arrows(0, 0, varsig[,1], varsig[,2], length = 0.08, angle = 20)
for (v in 1:nrow(varsig)) {
if (abs(varsig[v,1]) > abs(varsig[v,2])) {
if (varsig[v,1] >= 0) pos <- 4 else pos <- 2
} else {
if (varsig[v,2] >= 0) pos <- 3 else pos <- 1
}
text(varsig[v,1], varsig[v,2], labels=rownames(varsig)[v], pos = pos)
}
对象间以及对象和变量间的关系
观测对象（葡萄酒）在排序空间中的位置及相互关系，以及对象和变量的关系。
#查看排序对象，数据结构类上上述ind <- get_mfa_ind(res.mfa)
ind
#例如查看对象排序坐标
head(ind$coord)
#展示对象在前两轴中的排序图
#如果不想在图中展示分类变量，可使用参数 invisible = 'quali.var'
fviz_mfa_ind(res.mfa, axes = 1:2, palette = c('#00AFBB', '#E7B800', '#FC4E07'),
habillage = 'Label', addEllipses = TRUE, ellipse.type = 'confidence', repel = TRUE)
分类变量将默认和对象一起展示在该图中。分类变量显示为黑色，Env1、Env2、Env3代表变量集中的土壤类别，Saumur、Bourgueuil和Chinon代表酒的产地。在该图中，对象（葡萄酒）按其产地着色显示。
若两种葡萄酒的品质较为相似，则它们在排序图中的位置相互接近。
可结合上文中变量集或变量对各排序轴的贡献，评价葡萄酒特征。例如，葡萄酒T1和T2在第二轴上对应较高的数值，而葡萄酒摇动前后的气味两组变量集对第二轴具有较高的贡献，那么可知T1和T2应当属于更为“芳香浓郁”类型。
以及评估分类变量和对象以及定量变量间的关系。例如，分类变量Env4在第二轴上具有较高坐标值，靠近葡萄酒T1和T2，暗示改土壤环境中的葡萄更具香甜品质。
#其它可选可视化方案，例如分别按两组分类变量特征，在图中对代表葡萄酒的点标注颜色fviz_ellipses(res.mfa, axes = 1:2, c('Label', 'Soil'), repel = TRUE)
在排序图中展示MFA的全局PCA与各组变量集PCA中对象得分的关系。
#MFA 图中对象位置与各组变量集 PCA 中的对象位置，前两轴为例展示fviz_mfa_ind(res.mfa, axes = 1:2, partial = 'all')
#可选择部分对象展示
fviz_mfa_ind(res.mfa, axes = 1:2, partial = c('1DAM', '1VAU', '2ING'))
图中黑点代表MFA排序空间内的对象坐标；与对象点相连的不同颜色的线则代表了各变量集各自的PCA中对象的坐标，它们的形心即为MFA中对象坐标。
偏轴图
偏轴（partial axes）代表了所有变量集各自PCA（或MCA）所得特征向量（排序轴）投影到MFA排序空间中的情况，图中圆圈的半径（等于1）代表标准化偏轴最大的长度。
#偏轴图，前两轴为例展示
fviz_mfa_axes(res.mfa, axes = 1:2)
通过该图可评估哪些变量集的特征对MFA排序空间具有更高的贡献，这些变量集的前两轴与MFA的前两轴是否高度相关，MFA的各维度主要受哪种类型的属性集支配等。
参考资料