背景

在用统计分析方法研究多变量的课题时,变量个数太多就会增加课题的复杂性。人们自然希望变量个数较少而得到的信息较多

在很多情形,变量之间是有一定的相关关系的,当两个变量之间有一定相关关系时,可以解释为这两个变量反映此课题的信息有一定的重叠

在许多领域的研究与应用中,通常需要对含有多个变量的数据进行观测,收集大量数据后进行分析寻找规律。多变量大数据集无疑会为研究和应用提供丰富的信息,但是也在一定程度上增加了数据采集的工作量。更重要的是在很多情形下,许多变量之间可能存在相关性,从而增加了问题分析的复杂性。如果分别对每个指标进行分析,分析往往是孤立的,不能完全利用数据中的信息,因此盲目减少指标会损失很多有用的信息,从而产生错误的结论。

因此需要找到一种合理的方法,在减少需要分析的指标同时,尽量减少原指标包含信息的损失,以达到对所收集数据进行全面分析的目的。由于各变量之间存在一定的相关关系,因此可以考虑将关系紧密的变量变成尽可能少的新变量,使这些新变量是两两不相关的,那么就可以用较少的综合指标分别代表存在于各个变量中的各类信息。主成分分析与因子分析就属于这类降维算法。

概念

PCA(Principal Component Analysis),即主成分分析方法,是一种使用最广泛的数据降维算法。PCA的主要思想是将n维特征映射到k维上,这k维是全新的正交特征也被称为主成分,是在原有n维特征的基础上重新构造出来的k维特征。PCA的工作就是从原始的空间中顺序地找一组相互正交的坐标轴,新的坐标轴的选择与数据本身是密切相关的。其中,第一个新坐标轴选择是原始数据中方差最大的方向,第二个新坐标轴选取是与第一个坐标轴正交的平面中使得方差最大的,第三个轴是与第1,2个轴正交的平面中方差最大的。依次类推,可以得到n个这样的坐标轴。通过这种方式获得的新的坐标轴,我们发现,大部分方差都包含在前面k个坐标轴中,后面的坐标轴所含的方差几乎为0。于是,我们可以忽略余下的坐标轴,只保留前面k个含有绝大部分方差的坐标轴。事实上,这相当于只保留包含绝大部分方差的维度特征,而忽略包含方差几乎为0的特征维度,实现对数据特征的降维处理。

讲人话

主成分分析法是运用“降维”思想,把多个指标变换成少数综合指标的多元统计方法,这里的综合指标就是主成分。每个主成分都是原始变量的线性组合,彼此相互独立,并保留了原始变量绝大部分信息。其本质是通过原始变量的相关性,寻求相关变量的综合替代对象,并且保证了转化过程中的信息损失最小。

背景和概念看不懂?没关系!我来用人话讲解一遍。

PCA主成分析可以用来解决什么问题:

  1. 存在很多个个体、很多个变量,你想通过变量来对个体进行区分,找出个体之间的差别;
  2. 变量太多太多了,不可能为了区分每两个个体之间都需要用上所有的变量,变量与变量之间反应的信息也许还会有所重叠。

举个栗子

图片来源:恶霸小猴子

关于一个人的描述可以用很多中方式,但是很多描述其实是有所重叠的,我们可以把一个人的所有形容词归为几大类(降维),这样的话我我们对一个人进行描述的话仅需要使用这几个大类就可以了。

数学上的表达

图片来源:同济小旭学长

对于二维空间中点的描述需要两个坐标(x,y),如何对坐标轴进行变换,是的对这些点的位置的描述仅需要一维数据(降维)?

PCA主成分析的目标:只保留一个轴的时候(二维降到一维),信息保留最多。

怎么样最好:找到数据分布最分散的方向(方差最大)作为主成分(坐标轴)

首先:去中心化(把坐标原点放在数据中心),也就是数据的中心化与标准化,可以看这篇文章:数据中心化与标准化

然后:找坐标系(找到方差最大的方向),问题就是怎么找到方差最大的方向?(这一步是通过算法实现的,没看懂,需要数学线性代数的知识,我的都还给老师了,不过对于非计算机专业的,只是想使用PCA画图,这一步完全可以跳过)

找到一个能够反应尽量多的点信息的坐标,同时也要使得在新坐标下,点与点之间的差异性可以得到最大的保留。

找坐标的过程
找坐标的过程

主成分与原始变量之间的关系:

  • 主成分是原始变量的线性组合;
  • 主成分的数量相对于原始数量更少;
  • 主成分保留了原始变量的大部分信息。

PCA算法步骤总结

在进行之前最好先检验下数据之间的相关性:

首先进行KMO和Bartlett的检验,判断是否可以进行主成分分析。 对于KMO值:0.8上非常合适做主成分分析,0.7-0.8之间一般适合,0.6-0.7之间不太适合,0.5-0.6之间表示差,0.5下表示极不适合,对于 Bartlett的检验(p < 0.05,严格来说p < 0.01),若显著性小于0.05或0.01,拒绝原假设,则说明可以做主成分分析,若不拒绝原假设,则说明这些变量可能独立提供一些信息,不适合做主成分分析。

from factor_analyzer.factor_analyzer import calculate_bartlett_sphericity

# 巴特利特P值小于0.01,KMO值大于0.6;说明此数据适合做因子分析。
chi_square_value, p_value = calculate_bartlett_sphericity(data_)  # 计算巴特利特P值
chi_square_value, p_value

from factor_analyzer.factor_analyzer import calculate_kmo

kmo_all, kmo_model = calculate_kmo(data_)  # KMO检验
kmo_model

目标:把n维的数据降为k维。

步骤:

1、输入数据矩阵$X \in R^{m \times n}$,每一行是一个样本$x^{(i)}$;

2、必须地,去均值(中心化),使$X$每一列均值为0;

3、可选地,特征归一化(标准化);

X = (data - data.mean()) / data.std(ddof=0)  # 使用ddof=1样本标准偏差或ddof=0总体标准偏差作为参数来控制。

4、计算协方差矩阵$\sum \sum^{m}_{i=1} (x^{(i)})^{T}x^{(i)} = \frac{1}{m} X^T X$;

SIGMA = (data_.T @ data_) / data_.shape[0]  # 计算协方差矩阵(代表原矩阵中各列之间的相关性),@用来计算矩阵之间的乘法
# SIGMA  # 各列之间的相关性

5、计算奇异值分解:

U, S, V = np.linalg.svd(SIGMA)
# U

6、投影矩阵:

P = U[:, :2]  # 这里取U的前两列就是前两个投影的方向
# P

7、降维后的数据矩阵$Z = XP$;

Z = X @ P
Z

输出:降维后的矩阵$Z$,投射矩阵$P$。

泰坦尼克号数据集举例分析

讲解:颜色红色代表生还,蓝色代表死亡。可以看到红色(生还)多集中在横轴PC0<0区域,观察降维后的数据,符合横轴PC0<0的对应原7维中的Age、SibSp、Parch、Fare,解释Age值越大(大部分中年人)、SibSp(兄弟姐妹越多)、Fare(票价越贵)生还率越高;蓝色点(死亡率)多集中在横轴PC0>0区域,原7维中的Pclass、Sex经过降维后PC0>0,说明Pclasss(船舱等级越高,3代表三等舱,是最便宜的)、Sex(男性为1,女性为0)越为男性,死亡率越高。

R语言实现代码

> install.packages("devtools")  # 安装包
> library(devtools)  # 加载上一步安装的依赖包
> install_github("vqv/ggbiplot") # 这个包是基于上个包的,因此要先安装第一个包并加载后才能安装这个包
> library(ggbiplot) # 加载上一步安装的依赖包
> setwd("C:/Users/myxc/Desktop/") # 设置R的工作文件夹
> getwd()  # 查看上一步设置是否成功(非必要步骤,一般都是可以成功的,只是检查一下)
[1] "C:/Users/myxc/Desktop"
> df = read.csv("C:/Users/myxc/Desktop/demo.csv", header= T,row.names=1) # 读取数据
> df  # 查看数据
   site_01 site_02 site_03 site_04
0  5.70225 2.31307 8.07974 3.34235
1  1.42797 6.34185 5.56293 1.21440
2  4.73692 6.45518 1.52564 5.61749
3  8.67942 0.85052 3.48302 2.37014
4  9.66726 0.67532 3.76219 8.67886
5  0.98044 4.69239 3.19871 8.72567
6  5.94397 8.58706 4.36327 4.18807
7  0.80913 3.48337 1.61548 4.77624
8  0.81183 2.53953 7.32616 5.59256
9  3.43219 8.91296 7.78453 2.23462
10 4.85237 5.88927 3.69173 6.79233
11 1.80877 7.10859 5.34309 9.34172
12 5.63572 9.56429 6.24929 1.41587
13 4.47662 5.21469 9.82324 8.81561
14 1.41973 3.60362 2.80987 7.70326
15 7.18351 2.48648 1.65911 9.97133
16 4.70689 6.98316 7.48182 5.75968
17 0.44505 1.74083 3.13251 0.21471
> group=c(rep("OH",9),rep("OL",9))  # 设置分组,在PCA图上,同一个分组的数据会被圈到一个圈圈里(csv文件一共18条数据,前9条数据为“OH”分组,后9条数据为“OL”分组)
> result <- prcomp(df, scale = TRUE)  # 将数据导入PCA算法进行加工
> result  # 查看经过算法加工后可以制图的数据
Standard deviations (1, .., p=4):
[1] 1.2124401 1.0014478 0.9296445 0.8141575

Rotation (n x k) = (4 x 4):
               PC1       PC2        PC3         PC4
site_01 -0.3218732 0.8189381 -0.3556268 -0.31506772
site_02  0.6170596 0.1264450  0.3479028 -0.69441544
site_03  0.5484323 0.5179058  0.1253226  0.64442987
site_04 -0.4635269 0.2124281  0.8583632  0.05683011
> ggscreeplot(result)  # 绘制碎石图
> ggbiplot(result, obs.scale = 1, var.scale = 1, groups =group, ellipse = TRUE, circle = TRUE, var.axes = F)  # 绘制PCA图
> ggbiplot(result, obs.scale = 1, var.scale = 1, groups =group, ellipse = TRUE,circle = TRUE, var.axes = F) + scale_color_brewer(palette = "Set1") + theme(legend.direction = 'horizontal', legend.position = 'top')  # 绘制可以自定义美化的PCA图
PCA图
PCA图

碎石图横轴前两项的数据36.8%、25.1%分别对应PCA图的X轴和Y轴。

PCoA

PCoA(Principal Co-ordinates Analysis)分析即主坐标分析,可呈现研究数据相似性或差异性的可视化坐标,是一种非约束性的数据降维分析方法,可用来研究样本群落组成的相似性或相异性。它与PCA类似,通过一系列的特征值和特征向量进行排序后,选择主要排在前几位的特征值,找到距离矩阵中最主要的坐标,结果是数据矩阵的一个旋转,它没有改变样本点之间的相互位置关系,只是改变了坐标系统。

两者之间的区别:PCA是基于样本的相似系数矩阵(如欧式距离)来寻找主成分,而PCoA是基于距离矩阵(欧式距离以外的其他距离)来寻找主坐标。

过程详解

1、综合评价:

首先进行KMO和Bartlett的检验,判断是否可以进行主成分分析。 对于KMO值:0.8上非常合适做主成分分析,0.7-0.8之间一般适合,0.6-0.7之间不太适合,0.5-0.6之间表示差,0.5下表示极不适合,对于 Bartlett的检验(p < 0.05,严格来说p < 0.01),若显著性小于0.05或0.01,拒绝原假设,则说明可以做主成分分析,若不拒绝原假设,则说明这些变量可能独立提供一些信息,不适合做主成分分析;

适合性检验(KMO检验):

值的范围因子分析适合情况
>0.9非常适合
0.8~0.9很适合
0.7~0.8适合
0.6~0.7勉强适合
0.5~0.6不太适合
<0.5不适合

2、通过分析方差解释表格和碎石图,确定主成分的数量方差解释表格主要是看主成分对于变量解释的贡献率(可以理解为究竟需要多少主成分才能把变量表达为100%),如果太低(如低于60%)则需要调整主成分数据,碎石图的作用是根据特征值下降的坡度来确认需要选择的主成分个数,这两者结合可用于确认或调整主成分个数;

  • 保留的主成分使得方差贡献率达到85%以上;
  • 保留的主成分的方差(特征值)大于1;
  • 碎石图绘制了关于各主成分及其特征值的图形,我们只需要保留图形中变化最大之处以上的主成分即可。

3、通过分析主成分载荷系数与热力图,可以分析到每个主成分中隐变量的重要性,如研究【多金属矿体】中25种有用元素的分布规律,其中各元素视为指标,假设前文确定得到5个主成分,主成分1中,SO、SO2、Na2S、HS、H2S主成分载荷系数较大,因此可将主成分1确定为硫化物成分,以此类推,也可结合具体业务进行各主成分的隐变量分析;

4、基于主成分载荷图通过将多主成分降维成双主成分或者三主成分,通过象限图的方式呈现主成分的空间分布。如果提取2个主成分时,无法呈现三维载荷主成分散点图,如果提取1个主成分时,无法显示主成分象限图;

5、通过分析成分矩阵,得出主成分成分公式与权重;

6、输出主成分分析法综合得分。

注意事项

  1. 主成分要求变量之间的共线性或相关关系比较强,否则不能通过 KMO 检验和 Bartlett 球形检验;
  2. 主成分分析倾向于降维,从而达到简化系统结构,抓住问题实质的目的。(可侧重于输出结果 2、输出结果 3、输出结果 8);
  3. 主成分分析时通常需要综合自己的专业知识,以及软件结果进行综合判断,即使是特征根值小于 1,也一样可以提取主成分;
  4. KMO 值为 null 不存在可能导致的原因为:

    • 样本量过少容易导致相关系数过高,一般希望分析样本量大于 5 倍分析项个数;
    • 各个分析项之间的相关关系过高或过低。

参考文献

[1] Scientific Platform Serving for Statistics Professional 2021. SPSSPRO. (Version 1.0.11)[Online Application Software]. Retrieved from https://www.spsspro.com.

[2] 何晓群.多元统计分析.北京:中国人民大学出版社,2012.

[3] 王 伟,赵 明.主成分分析法在航材分类指标体系构建中的应用[J].舰船电子工程,2019,39 (1): 118-120.

[4] 丁敬国,郭锦华. 基于主成分分析协同随机森林算法的热连轧带钢宽度预测[J]. 东北大学学报(自然科学版)2021,42(9):1268-1274,1289.

[5] 曹晓俊. 对我国上市银行经营业绩的分析——基于主成分分析、因子分析和聚类分析的方法 [J]. 宿州学院学报, 2016, 31(7): 5.