主成分分析|因子分析

在处理实际问题中，在很多情况下，不同指标之间可能存在一定的相关性，即变量反映样本的信息有一定重叠
由于指标较多再加上信息的重叠，势必增加了问题的复杂性
希望通过克服相关性、重叠性，用较少的变量代替原来较多的变量，而这些较少的变量可以反映原来众多变量的大部分信息，这实际上是一种降维的思想

主成分分析

主成分分析是一种通过降维技术将多个变量化为少数几个主成分的统计分析方法。这些主成分通常表示为原始变量的某种线性组合，能够反映原始变量的绝大部分信息
主成分分析在几何上是一个坐标转换，在二维空间中有明显的几何意义

示例

假设共有n个样品，每个样品有两个指标$x_1,x_2$，在坐标系中，这n个点沿$x_1$和$x_2$方向都有较大的分量，若仅考虑其中一个分量，包含在另一分量中的信息将会损失，因此直接舍弃某个分量不是降维的有效方法
若将该坐标系按逆时针方向旋转某个角度$\theta$编程新坐标系，变换公式为
$$\{\begin{array}{c}{Y}_1=X_1\cos \theta +X_2\sin \theta \\ Y_2=-X_1\sin \theta +X_2\cos \theta \end{array}$$

总体主成分的数学模型

设$X=(X_1,X_2,\dots ,X_P{)}^T$为P维随机向量
$$cov(X)=\sum =({\sigma }_{ij}{)}_{p\times p}=E[(X-E(x))(X-E(x){)}^T]$$
为其协方差矩阵，是非负定矩阵
构造的线性组合
$$Y_1=a_1^{T}X=a_{11}{X}_1+a_{12}{X}_2+\cdots +a_{1p}{X}_p$$
确定
$$a_1=(a_{11},a_{12},\dots ,a_{1p}{)}^T$$
$$Var(Y_1)=Cov(a_1^{T}X,a_1^{T}X)=a_1^T\sum a_1$$
使得$Var(Y_1)$达到最大

• 必须对$a_1$加以限制，否则$Var(Y_1)$无界。如$a_1^T{a}_1=1$
  如此确定的随机变量$Y_1=a_1^{T}X$称为X的第一主成分

如果第一主成分在$Y_1$方向上的分散性还不足以反映原变量的分散性
则构造
$$Y_2=a_2^{T}X=a_{21}{X}_1+a_{22}{X}_2+\cdots +a_{2p}{X}_p$$
限制$a_2^T{a}_2=1$
满足
$$Cov(Y_1,Y_2)=Cov(a_1^{T}X,a_2^{T}X)=a_1^T\sum a_2=0$$
$a_2$的模是1，$Y_1,Y_2$不相关
使得$Var(Y_2)$达到最大
如此确定的随机变量$Y_2$称为$X$的第二主成分
一般地$Y_1,Y_2,\dots ,Y_{k-1}$还不足以反映原变量的信息
则构造
$$Y_k=a_k^{T}X=a_{k1}^T{X}_1+a_{k2}^T{X}_2+\cdots +a_{kp}{X}_p$$
限制$a_k^T{a}_k=1$
$$Cov(Y_k,Y_i)=Cov(a_k^{T}X,a_i^{T}X)=a_k^T\sum a_i=0(i=1,2,\dots ,k-1)$$
使得$Var(Y_k)$达到最大
如此确定的随机变量$Y_k$称为X的第k主成分

按照上述方法，最多可以构造出p个方差大于0的主成分

总体主成分的求法

$\sum$是$X=(X_1,X_2,\dots ,X_p{)}^T$的协方差矩阵，其特征值按大小顺序排列为${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_p\ge 0$，相应的正交单位化特征向量为$e_1,e_2,\dots ,e_p$，则X的第k个主成分为
$$Y_k=e_k^{T}X=e_{k1}{X}_1+e_{k2}{X}_2+\cdots +e_{kp}{X}_p$$
其中
$$e_k=(e_{k1},e_{k2},\dots ,e_{kp}{)}^T$$
这时
$$Var(Y_k)=e_k^T\sum e_k={\lambda }_k{e}_k^T{e}_k={\lambda }_k,k=1,2,\dots ,p$$
$$Cov(Y_j,Y_k)=Cov(e_j^T,e_k^{T}X)=e_j\sum e_k={\lambda }_k{e}_j^T{e}_k=0,j\ne k$$

总体主成分的性质

主成分的协方差矩阵及总方差

记$Y=(Y_1,Y_2,\dots ,Y_P{)}^T$为p个主成分构成的随机向量，则$Y=P^{T}X$
其中$P=(e_1,e_2,\dots ,e_p)$
其协方差矩阵为
$$Cov(Y)=Cov(P^{T}X)=P^T\sum P=\Lambda =Diag({\lambda }_1,\dots ,{\lambda }_p)$$
各主成分的总方差为
$$\begin{array}{c}{\sum }_{i=1}^{P}Var(Y_1)={\sum }_{i=1}^P{\lambda }_i=tr\left(P^T\sum P\right)\\ =tr\left(\sum P{P}^T\right)=tr\left(\sum \right)={\sum }_{i=1}^{P}Var(X_i)\end{array}$$
即主成分分析把p这个原始变量$X_1,X_2,\dots ,X_p$的总方差分解成p个不相关变量$Y_1,Y_2,\dots ,Y_p$的方差和，且使得
$$Var(Y_k={\lambda }_k(k=1,2,\dots ,p)$$

主成分的贡献率和累计贡献率

$$\frac{{\lambda }_k}{{\sum }_{i=1}^p{\lambda }_i}=\frac{Var(Y_k)}{{\sum }_{i=1}^{p}Var(X_i)}$$
描述了第k个主成分提取的$X_1,X_2,\dots ,X_p$的总信息的份额，我们称此为第k个主成分$Y_k$的贡献率
前m个主成分的贡献率之和
$$\frac{{\sum }_{i=1}^m{\lambda }_i}{{\sum }_{i=1}^p{\lambda }_i}$$
称为$Y_1,Y_2,\dots ,Y_m$的累计贡献率，它表明前m个主成分综合的$X_1,X_2,\dots ,X_p$信息的能力

例1

解
算得$\sum$的特征值及单位正交特征向量
$$\begin{array}{c}{\lambda }_1=5.8284,{\lambda }_2=2.0000,{\lambda }_3=0.1716\\ e_1=(-0.3827,0.9239,0.0000{)}^T,Y_1=-0.3827X_1+0.9239X_2\\ e_2=(0.0000,0.0000,1.0000{)}^T,Y_2=X_3\\ e_3=(0.9239,0.3827,0.0000{)}^T,Y_3=0.9239X_1+0.3827X_2\end{array}$$
第一主成分的贡献率为
$$\frac{5.83}{5.83+2.00+0.17}=\frac{5.83}{1+5+2}=73\%$$
前两个主成分的累计贡献率为
$$\frac{5.83+2.00}{5.83+2.00+0.17}=\frac{5.83+2.00}{1+5+2}=98\%$$
因此，用前两个主成分代替原来的三个变量，其信息的损失仅为2%，是很小的

标准化变量的主成分

在实际问题中，不同的变量有不同的量纲，由于不同的量纲会引起各变量取值的分散性程度差异较大
各变量取值的大小也会影响其方差的大小
这时，变量的总方差主要受方差较大的变量控制，若进行主成分分析，则优先照顾了方差较大的变量
会给主成分的解释带来困难，有时还会造成不合理的结果
通常将原变量进行标准化再做主成分分析

设
$${\mu }_k=E(X_k),{\sigma }_{kk}=Var(X_k),k=1,2,\dots ,p$$
标准化变量
$$X_k^{\ast }=\frac{X_k-{\mu }_k}{\sqrt{{\sigma }_{kk}}},k=1,2,\dots ,p$$
这时一切对于$1\le k\le p$均有$Var(X_k^{\ast })=1$
令
$$X^{\ast }=(X_1^{\ast },X_2^{\ast },\dots ,X_p^{\ast }{)}^T$$
则其协方差矩阵$\rho =({\rho }_{ij}{)}_{p\times p}=Cov(X^{\ast })$
为X的相关系数矩阵
$${\rho }_{ij}=E(X_i^{\ast },X_j^{\ast })=\frac{Cov(X_i,X_j)}{\sqrt{{\sigma }_{ii}{\sigma }_{jj}}}$$
X的第i个主成分为
$$Y_k^{\ast }=(e_k^{\ast }{)}^T{X}^{\ast }=e_{k1}^{\ast }{X}_1^{\ast }+\cdots +e_{kp}^{\ast }{X}_p^{\ast }$$
并且
$$\sum _{i=1}^{p}Var(Y_i^{\ast })=\sum _{i=1}^p{\lambda }_i^{\ast }=\sum _{i=1}^{p}Var(X_i^{\ast })=p$$
其中${\lambda }_1^{\ast }\ge {\lambda }_2^{\ast }\ge \cdots \ge {\lambda }_p^{\ast }$为$\rho$的特征值
$e_i^{\ast }=(e_{i1}^{\ast },e_{i2}^{\ast },\dots ,e_{ip}^{\ast })$为相应于${\lambda }_i^{\ast }$的单位正交特征向量
第i个主成分贡献率为$\frac{{\lambda }_i^{\ast }}{P}$
前m个主成分累计贡献为${\sum }_{i=1}^m\frac{{\lambda }_i^{\ast }}{P}$

例2

解
算得$\sum$的特征值及单位正交特征向量
$${\lambda }_1=100.16,{\lambda }_2=0.84$$
$$e_1=(0.040,0.999{)}^T,Y_1=0.040X_1+0.999X_2$$
$$e_2=(0.999,-0.040{)}^T,Y_2=0.999X_1-0.040X_2$$
此时，第一主成分的贡献率为$\frac{100.16}{2}=99.2\%$
$$e_1^{\ast }=(0.707,0.707{)}^T$$
$$Y_1^{\ast }=0.707X_1^{\ast }+0.707X_2^{\ast }=0.707(X_1-{\mu }_1)+0.0707(X_2-{\mu }_2)$$
$$e_2^{\ast }=(0.707,-0.707{)}^T$$
$$Y_2^{\ast }=0.707X_1^{\ast }-0.707X_2^{\ast }=0.707(X_1-{\mu }_1)-0.0707(X_2-{\mu }_2)$$
此时，第一主成分的贡献率为
$$\frac{1.4}{2}=70\%$$

样本主成分

实际问题中$\sum$或$\rho$未知，需用样本估计
设
$$x_i=(x_{i1},x_{i2},\dots ,x_{ip}{)}^T,i=1,2,\dots ,n$$
容量为n的样本，用样本协方差矩阵及样本相关系数矩阵
$$S=\frac{1}{n-1}\sum _{k=1}^n(x_k-\bar{x})(x_k-\bar{x}{)}^T$$
$$R=(r_{ij}{)}_{p\times p}=\left(\frac{S_{ij}}{\sqrt{S_{ii}}\sqrt{S_{jj}}}\right)=(S_{ij}{)}_{p\times p}$$
估计$\sum ,\rho$
其中
$$\bar{x}=({\bar{x}}_1,{\bar{x}}_2,\dots ,{\bar{x}}_p{)}^T,{\bar{x}}_j=\frac{1}{n}\sum _{i=1}^n{x}_{ij},j=1,2,\dots ,p$$
$$S_{ij}=\frac{1}{n-1}\sum _{k=1}^n(x_{ki}-{\bar{x}}_i)(x_{kj}-{\bar{x}}_j),i,j=1,2,\dots ,p$$
S为样本协方差矩阵，它的特征值为$$\begin{gathered} {\hat{\lambda }}_1\ge \\ hat{lambda}_2\ge \cdots \ge {\hat{\lambda }}_p\ge 0 \end{gathered}$$，相应的正交单位化特征向量为${\hat{e}}_1,{\hat{e}}_2,\dots ,{\hat{e}}_p$
样本主成分
$$Y_k=e_k^{T}X=e_{k1}{X}_1+e_{k2}{X}_2+\cdots +e_{kp}{X}_p,k=1,2,\dots ,p$$
n个观测值$X_i=(x_{i1},x_{i2},\dots ,x_{ip}{)}^T(k=1,2,\dots ,n)$代入第k个样本主成分，得$Y_i$得n个观测值$Y_{1i},Y_{2i},\dots ,Y_{ni}$称为第k个样本主成分的得分
第k个主成分的贡献率为
$$\frac{{\hat{\lambda }}_k}{{\sum }_{j=1}^p{\hat{\lambda }}_k}$$
前m个主成分的累计贡献率为
$$\frac{{\sum }_{j=1}^m{\hat{\lambda }}_k}{{\sum }_{j=1}^p{\hat{\lambda }}_k}$$
通常从样本相关矩阵R出发进行主成分分析，从$\rho$出发求得的主样本成分称标准化样本主成分

例3

解
从样本协方差矩阵出发进行主成分分析，算得
$${\hat{\lambda }}_1=100.004139{\hat{e}}_1=(0.559147,0.421287,0.714046{)}^T$$
$${\hat{\lambda }}_2=25.324480{\hat{e}}_2=(0.827674,-0.333483,-0.451382{)}^T$$
$${\hat{\lambda }}_3=1.568048{\hat{e}}_3=(0.047960,0.843390,-0.535157{)}^T$$
$$Y_1=0.559147x_1+0.421287x_2+0.714046x_3$$
$$Y_2=0.827674x_1-0.333483x_2-0.451382x_3$$
各样本主成分的贡献率为
$$\frac{100.004}{136.896}=80.36\%\frac{25.324}{136.896}=18.52\%\frac{1.568}{136.896}=1.15\%$$
前两个主成分的贡献率达到了98.86%，实际运用中可以只取前两个主成分
当一个学生的$y_1$较大时，可以推断他较高或较胖或又高又胖；反之，当一个学生比较魁梧时，所对应的$y_1$也比较大，所以第一成分是反映学生身材是否魁梧的综合指标，我们一般称之为大小因子
同理，第二主成分反映学生体型特征的综合指标，我们一般称之为瘦高还是胖矮因子

主成分分析中注意的问题

1. 一般要求累计贡献流程达到85%，常见的情况2~3个主成分
2. 一般要将各原始变量作标准化处理，也就是相关系数矩阵，通常从相关系数矩阵出发进行主成分分析

主成分分析过程

1. 将原始数据标准化
2. 建立相关系数矩阵
3. 计算相关系数阵的特征值${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_p\ge 0$，对应的单位正交特征向量$a_1,a_2,\dots ,a_p$
4. 由累计贡献率确定主成分个数m，并写出主成分
   $$z_i=a_{i}x,i=1,2,\dots ,m$$

主成分分析的应用

综合评价

主成分分析能从选定的指标体系中归纳出大部分信息，可以根据主成分提供的信息进行综合评价
利用主成分进行综合评价时，对主成分进行加权综合，权数根据其方差贡献率确定
$$z={\rho }_1{z}_1+{\rho }_2{z}_2+\cdots +{\rho }_m{z}_m$$
本质上，综合评价函数是对原始指标的线性综合，从计算主成分到加权，经过两次线性运算后得到综合评价函数

主成分回归

将各主成分作为新自变量代替原来自变量作回归分析。用主成分分析筛选变量，可用较少的计算量获得选择最佳变量子集合的效果，并且由于各主成分两两不相关，不存在由于多重共线性带来的影响

主成分聚类

主要用于样品的聚类，先构建主成分，再进行聚类
保留多少主成分取决于保留部分的累计方差在方差总和中所占百分比，它包括着前几个主成分概括信息的多少。
在实践中，粗略规定一个百分比便可以决定保留几个主成分，如果多保留一个主成分，但累计方差增加无几，便不再保留
设${\lambda }_1,{\lambda }_2,\dots ,{\lambda }_p$为主成分的特征值，则前K个方差累计贡献率为
$$\frac{{\lambda }_1+{\lambda }_2+\cdots +{\lambda }_k}{{\sum }_{{\lambda }_i}}$$
一般累计方差贡献率大于85%时不在增加新的主成分

因子分析法

因子分析法概述

与主成分分析法都基于统计分析法
主成分分析是通过坐标变换提取主成分，将一组具有相关性的变量变换为一组独立的变量，将主成分表示为原始观察变量的线性组合
因子分析法是要构造因子模型，将原始观察变量分解为因子的线性组合，因此因子分析法是主成分分析法的发展

因子分析法的模型

设有n个样品，每个样品观察p个指标，这p个指标有较强的相关性
为了消除由于观察量纲的差异以及数量级不同所造成的影响，对样本观察数据进行标准化处理，使标准化后的变量均值为0，方差为1
用x表示标准化后的变量向量，$F_1,F_2,\dots ,F_m$表示标准化的公共因子
因子模型的表达式
$$x_1=a_{11}{F}_1+a_{12}{F}_2+\cdots +a_{1m}{F}_m+e_1$$
$$x_2=a_{21}{F}_1+a_{22}{F}_2+\cdots +a_{2m}{F}_m+e_2$$
$$\dots \dots$$
$$x_p=a_{p1}{F}_1+a_{p2}{F}_2+\cdots +a_{pm}{F}_m+e_p$$

因子分析法的假设

$X=(X_1,X_2,\dots ,X_p{)}^{\prime }$是可观测随机向量，其均值向量$E(X)=0$，协方差矩阵$cov(X)=\sum$，且协方差矩阵$\sum$和相关矩阵$R$相等
$F=(F_1,F_2,\dots ,F_p{)}^{\prime }(m<p)$是不可观测随机向量，且均值向量$E(F)=0$，协方差矩阵$cov(X)=I$，即向量$F$的各分量是相关独立的
$\epsilon =({\epsilon }_1,{\epsilon }_2,\dots ,{\epsilon }_p{)}^{\prime }$与$F$相互独立，且$E(\epsilon )=0$，其协方差矩阵为对角阵
其矩阵形式为
$$X=AF+e$$
其中$A=(a_{ij})$，$a_{ij}$为因子载荷。
数学上可以证明，因子载荷$a_{ij}$就是第i变量与第j因子的相关系数，反映了第i变量在第j因子上的重要性
F称为X的公共因子或潜因子，矩阵A称为因子载荷矩阵，e称为X的特殊因子

变量的共同度(公共方差)

反应全部公共因子对原有变量$x_i$总方差的解释说明比例
原有变量$x_i$的共同度为因子载荷矩阵$A$的第i行元素的平方和
即
$$h_i^2=\sum _{j=1}^m{a}_{ij}^2$$
$h_i^2$越接近1，原有变量在标准化下的方差为1，说明公共因子解释原有变量$x_i$的信息越多

公共因子$F_j$的方差贡献
公共因子$F_j$的方差贡献定义为因子载荷矩阵$A$中的第$j$列元素的平方和
即
$$g_j^2=\sum _{i=1}^p{a}_{ij}^2$$
反映了公共因子$F_j$对原有变量总方差的解释能力，其值越高说明因子的重要程度越高

因子载荷的求解

1. 主成分法
2. 主轴因子法
3. 极大似然法

因子旋转

不管用何种方法确定的初始因子载荷矩阵，他们都不是唯一的
设$F_1,F_2,\dots ,F_m$是初始公共因子，则可以建立他们如下的线性组合得到一组新的公共因子$F_1^{\prime },F_2^{\prime },\dots ,F_m^{\prime }$
$$F_1^{\prime }=d_{11}{F}_1+d_{12}{F}_2+\cdots +d_{1m}{F}_m$$
$$F_2^{\prime }=d_{21}{F}_1+d_{22}{F}_2+\cdots +d_{2m}{F}_m$$
$$\dots \dots$$
$$F_m^{\prime }=d_{m1}{F}_1+d_{m2}{F}_2+\cdots +d_{mm}{F}_m$$
使得$F_1^{\prime },F_2^{\prime },\dots ,F_m^{\prime }$彼此相互独立，同时也可以很好地解释原始变量之间的相互关系
因子旋转分为正交旋转和斜交旋转
正交旋转由初始载荷矩阵右乘一正交阵得到，经过正交旋转而得到的新的公共因子任然保持相互独立的性质
斜交旋转放弃了因子之间彼此独立这个限制，因而可能达到更为简洁的形式，实际意义也更容易解释

因子得分

当因子模型建立之后，往往需要反过来考察每一个样品的性质及样品之间的相互关系
因子得分就是公共因子$F_1,F_2,\dots ,F_m$在每一个样品上的得分

因子得分的计算

用回归的思想求出线性组合系数的估计值
即：建立如下公共因子为因变量，原始变量为自变量的回归方程
$$F_j={\beta }_{j1}{X}_1+{\beta }_{j2}{X}_2+\cdots +{\beta }_{jp}{X}_p,j=1,2,\dots ,m$$
此处因为原始变量与公共因子变量均为标准化变量，所以回归模型中不存在常数项
在最小二乘意义下，可以得到的F估计值
$$\hat{F}=A^{\prime }{R}^{-1}X$$
A，因子载荷矩阵
R，原始变量的相关矩阵
X，原始变量向量

因子分析的步骤

1. 根据研究问题选取原始变量
2. 对原始变量进行标准化，并求其相关矩阵，分析变量之间的相关性
3. 求解初始公共因子及其因子载荷矩阵
4. 因子旋转
5. 计算因子得分
6. 根据因子得分进行进一步分析

Matlab

[pc, score] = princomp(x)
[pc, score, latent] = princomp(x)
[pc, latent, explain] = pcacova(covx)

pc，主成分系数
score，原x矩阵在主成分空间的表示
latent，x协方差矩阵的特征值
explain，方差贡献率