主成分分析和因子分析 PC, FA 及运行结果
欢迎投稿(荐稿)计量经济圈，计量相关都行
邮箱：econometrics666@sina.cn
编辑: @计量经济圈(ID: econometrics666)
征募令：只要你能够整理二篇文献发到计量经济圈公众号，你就可以直接进入我们的大社群，无需经过其他程序。
主成分分析
在许多领域的研究与应用中，往往需要对反映事物的多个变量进行大量的观测，收集大量数据以便进行分析寻找规律。多变量大样本无疑会为研究和应用提供了丰富的信息，但也在一定程度上增加了数据采集的工作量，更重要的是在多数情况下，许多变量之间可能存在相关性，从而增加了问题分析的复杂性，同时对分析带来不便。如果分别对每个指标进行分析，分析往往是孤立的，而不是综合的。盲目减少指标会损失很多信息，容易产生错误的结论。
因此需要找到一个合理的方法，在减少需要分析的指标同时，尽量减少原指标包含信息的损失，以达到对所收集数据进行全面分析的目的。由于各变量间存在一定的相关关系，因此有可能用较少的综合指标分别综合存在于各变量中的各类信息。主成分分析与因子分析就属于这类降维的方法。
主成分分析是设法将原来众多具有一定相关性（比如P个指标），重新组合成一组新的互相无关的综合指标来代替原来的指标。
主成分分析，是考察多个变量间相关性一种多元统计方法，研究如何通过少数几个主成分来揭示多个变量间的内部结构，即从原始变量中导出少数几个主成分，使它们尽可能多地保留原始变量的信息，且彼此间互不相关.通常数学上的处理就是将原来P个指标作线性组合，作为新的综合指标。
最经典的做法就是用F1（选取的第一个线性组合，即第一个综合指标）的方差来表达，即Var(F1)越大，表示F1包含的信息越多。因此在所有的线性组合中选取的F1应该是方差最大的，故称F1为第一主成分。如果第一主成分不足以代表原来P个指标的信息，再考虑选取F2即选第二个线性组合，为了有效地反映原来信息，F1已有的信息就不需要再出现在F2中，用数学语言表达就是要求Cov(F1, F2)=0，则称F2为第二主成分，依此类推可以构造出第三、第四，……，第P个主成分。
2. 问题描述
下表1是某些学生的语文、数学、物理、化学成绩统计：
首先，假设这些科目成绩不相关，也就是说某一科目考多少分与其他科目没有关系。那么一眼就能看出来，数学、物理、化学这三门课的成绩构成了这组数据的主成分（很显然，数学作为第一主成分，因为数学成绩拉的最开）。为什么一眼能看出来？因为坐标轴选对了！下面再看一组学生的数学、物理、化学、语文、历史、英语成绩统计，见表2，还能不能一眼看出来：
数据太多了，以至于看起来有些凌乱！也就是说，无法直接看出这组数据的主成分，因为在坐标系下这组数据分布的很散乱。究其原因，是因为无法拨开遮住肉眼的迷雾~如果把这些数据在相应的空间中表示出来，也许你就能换一个观察角度找出主成分。如下图1所示：
但是，对于更高维的数据，能想象其分布吗？就算能描述分布，如何精确地找到这些主成分的轴？如何衡量你提取的主成分到底占了整个数据的多少信息？所以，我们就要用到主成分分析的处理方法。
3. 数据降维
为了说明什么是数据的主成分，先从数据降维说起。数据降维是怎么回事儿？假设三维空间中有一系列点，这些点分布在一个过原点的斜面上，如果你用自然坐标系x,y,z这三个轴来表示这组数据的话，需要使用三个维度，而事实上，这些点的分布仅仅是在一个二维的平面上，那么，问题出在哪里？
如果你再仔细想想，能不能把x,y,z坐标系旋转一下，使数据所在平面与x,y平面重合？这就对了！如果把旋转后的坐标系记为x’,y’,z’，那么这组数据的表示只用x’和y’两个维度表示即可！当然了，如果想恢复原来的表示方式，那就得把这两个坐标之间的变换矩阵存下来。这样就能把数据维度降下来了！
但是，我们要看到这个过程的本质，如果把这些数据按行或者按列排成一个矩阵，那么这个矩阵的秩就是2！这些数据之间是有相关性的，这些数据构成的过原点的向量的最大线性无关组包含2个向量，这就是为什么一开始就假设平面过原点的原因！
那么如果平面不过原点呢？这就是数据中心化的缘故！将坐标原点平移到数据中心，这样原本不相关的数据在这个新坐标系中就有相关性了！有趣的是，三点一定共面，也就是说三维空间中任意三点中心化后都是线性相关的，一般来讲n维空间中的n个点一定能在一个n-1维子空间中分析！
上一段文字中，认为把数据降维后并没有丢弃任何东西，因为这些数据在平面以外的第三个维度的分量都为0。现在，假设这些数据在z’轴有一个很小的抖动，那么我们仍然用上述的二维表示这些数据，理由是我们可以认为这两个轴的信息是数据的主成分，而这些信息对于我们的分析已经足够了，z’轴上的抖动很有可能是噪声，也就是说本来这组数据是有相关性的，噪声的引入，导致了数据不完全相关，但是，这些数据在z’轴上的分布与原点构成的夹角非常小，也就是说在z’轴上有很大的相关性，综合这些考虑，就可以认为数据在x’,y’ 轴上的投影构成了数据的主成分！
课堂上老师谈到的特征选择的问题，其实就是要剔除的特征主要是和类标签无关的特征。而这里的特征很多是和类标签有关的，但里面存在噪声或者冗余。在这种情况下，需要一种特征降维的方法来减少特征数，减少噪音和冗余，减少过度拟合的可能性。
PCA的思想是将n维特征映射到k维上（k<n），这k维是全新的正交特征。这k维特征称为主成分，是重新构造出来的k维特征，而不是简单地从n维特征中去除其余n-k维特征。
二、PCA实例
现在假设有一组数据如下：
行代表了样例，列代表特征，这里有10个样例，每个样例两个特征。可以这样认为，有10篇文档，x是10篇文档中“learn”出现的TF-IDF，y是10篇文档中“study”出现的TF-IDF。
第一步，分别求x和y的平均值，然后对于所有的样例，都减去对应的均值。这里x的均值是1.81，y的均值是1.91，那么一个样例减去均值后即为（0.69,0.49），得到
第二步，求特征协方差矩阵，如果数据是3维，那么协方差矩阵是
这里只有x和y，求解得
对角线上分别是x和y的方差，非对角线上是协方差。协方差是衡量两个变量同时变化的变化程度。协方差大于0表示x和y若一个增，另一个也增；小于0表示一个增，一个减。如果ｘ和ｙ是统计独立的，那么二者之间的协方差就是０；但是协方差是０，并不能说明ｘ和ｙ是独立的。协方差绝对值越大，两者对彼此的影响越大，反之越小。协方差是没有单位的量，因此，如果同样的两个变量所采用的量纲发生变化，它们的协方差也会产生树枝上的变化。
第三步，求协方差的特征值和特征向量，得到
上面是两个特征值，下面是对应的特征向量，特征值0.0490833989对应特征向量为，这里的特征向量都归一化为单位向量。
第四步，将特征值按照从大到小的顺序排序，选择其中最大的k个，然后将其对应的k个特征向量分别作为列向量组成特征向量矩阵。
这里特征值只有两个，我们选择其中最大的那个，这里是1.28402771，对应的特征向量是(-0.677873399, -0.735178656)T。
第五步，将样本点投影到选取的特征向量上。假设样例数为m，特征数为n，减去均值后的样本矩阵为DataAdjust(m*n)，协方差矩阵是n*n，选取的k个特征向量组成的矩阵为EigenVectors(n*k)。那么投影后的数据FinalData为
FinalData(10*1) = DataAdjust(10*2矩阵) x 特征向量(-0.677873399, -0.735178656)T
得到的结果是：
这样，就将原始样例的n维特征变成了k维，这k维就是原始特征在k维上的投影。
上面的数据可以认为是learn和study特征融合为一个新的特征叫做LS特征，该特征基本上代表了这两个特征。上述过程如下图2描述：
正号表示预处理后的样本点，斜着的两条线就分别是正交的特征向量（由于协方差矩阵是对称的，因此其特征向量正交），最后一步的矩阵乘法就是将原始样本点分别往特征向量对应的轴上做投影。
整个PCA过程貌似及其简单，就是求协方差的特征值和特征向量，然后做数据转换。但是有没有觉得很神奇，为什么求协方差的特征向量就是最理想的k维向量？其背后隐藏的意义是什么？整个PCA的意义是什么？
三、PCA推导
先看下面这幅图：
在第一部分中，我们举了一个学生成绩的例子，里面的数据点是六维的，即每个观测值是6维空间中的一个点。我们希望将6维空间用低维空间表示。
先假定只有二维，即只有两个变量，它们由横坐标和纵坐标所代表；因此每个观测值都有相应于这两个坐标轴的两个坐标值；如果这些数据形成一个椭圆形状的点阵，那么这个椭圆有一个长轴和一个短轴。在短轴方向上，数据变化很少；在极端的情况，短轴如果退化成一点，那只有在长轴的方向才能够解释这些点的变化了；这样，由二维到一维的降维就自然完成了。
上图中，u1就是主成分方向，然后在二维空间中取和u1方向正交的方向，就是u2的方向。则n个数据在u1轴的离散程度最大（方差最大），数据在u1上的投影代表了原始数据的绝大部分信息，即使不考虑u2，信息损失也不多。而且，u1、u2不相关。只考虑u1时，二维降为一维。
椭圆的长短轴相差得越大，降维也越有道理。
1. 最大方差理论
在信号处理中认为信号具有较大的方差，噪声有较小的方差，信噪比就是信号与噪声的方差比，越大越好。如前面的图，样本在u1上的投影方差较大，在u2上的投影方差较小，那么可认为u2上的投影是由噪声引起的。
因此我们认为，最好的k维特征是将n维样本点转换为k维后，每一维上的样本方差都很大。
比如我们将下图中的5个点投影到某一维上，这里用一条过原点的直线表示（数据已经中心化）：
假设我们选择两条不同的直线做投影，那么左右两条中哪个好呢？根据我们之前的方差最大化理论，左边的好，因为投影后的样本点之间方差最大（也可以说是投影的绝对值之和最大）。
计算投影的方法见下图5：
图中，红色点表示样例，蓝色点表示在u上的投影，u是直线的斜率也是直线的方向向量，而且是单位向量。蓝色点是在u上的投影点，离原点的距离是<x,u>（即xTu或者uTx）。
2. 最小二乘法
我们使用最小二乘法来确定各个主轴（主成分）的方向。
对给定的一组数据（下面的阐述中，向量一般均指列向量）：
其数据中心位于:
数据中心化（将坐标原点移到样本点的中心点）：
中心化后的数据在第一主轴u1方向上分布散的最开，也就是说在u1方向上的投影的绝对值之和最大（也可以说方差最大），计算投影的方法上面已经阐述，就是将x与u1做内积，由于只需要求u1的方向，所以设u1也是单位向量。
在这里，也就是最大化下式：
由矩阵代数相关知识可知，可以对绝对值符号项进行平方处理，比较方便。所以进而就是最大化下式：
两个向量做内积，可以转化成矩阵乘法：
所以目标函数可以表示为：
括号里面就是矩阵乘法表示向量内积，由于列向量转置以后是行向量，行向量乘以列向量得到一个数，一个数的转置还是其本身，所以又可以将目标函数化为：
去括号：
又由于u1和i无关，可以拿到求和符外面，上式化简为：
学过矩阵代数的同学可能已经发现了，上式括号里面求和后的结果，就相当于一个大矩阵乘以自身的转置，其中，这个大矩阵的形式如下：
X矩阵的第i列就是xi
于是有：
所以目标函数最终化为：
其中的
我们假设
所以
下面我们求解最大值、取得最大值时u1的方向这两个问题。
先解决第一个问题，对于向量x的二范数平方为:
同样，目标函数也可以表示成映射后的向量的二范数平方：
把二次型化成一个范数的形式，由于u1取单位向量，最大化目标函数的基本问题也就转化为：对一个矩阵，它对一个向量做变换，变换前后的向量的模长伸缩尺度如何才能最大？我们有矩阵代数中的定理知，向量经矩阵映射前后的向量长度之比的最大值就是这个矩阵的最大奇异值，即：
式中，
针对本问题来说，
再解决第二个问题，对一般情况，设对称阵
的n个特征值分别为：
相应的单位特征向量为：
任取一个向量x，用特征向量构成的空间中的这组基表示为：
则：
所以：
针对第二个问题，我们取上式中的
证明完毕。
主成分所占整个信息的百分比可用下式计算：
式中分母为
有些研究工作表明，所选的主轴总长度占所有主轴长度之和的大约85% 即可，其实，这只是一个大体的说法，具体选多少个，要看实际情况而定。
3.意义
PCA将n个特征降维到k个，可以用来进行数据压缩，例如100维的向量最后可以用10维来表示，那么压缩率为90%。同样图像处理领域的KL变换使用PCA做图像压缩，人脸检测和匹配。
这里一大部分整理自：http://blog.csdn.net/zhongkelee/article/details/44064401
主成分分析Stata操作
读取数据：
主成分分析：
变量特征向量：
Screenplot:
预测前三个主成分：
预测出来的前三个主成分结果：
前三个主成分与原始变量的相关性：
因此，我们可以说第一个主成分变量与原始的8个变量都是相关的，而第二个主成分变量则只反映其中的3个变量，至于第三个主成分变量，它只反映1个变量headroom因素（我们把绝对值>0.5设为显著）。
因子分析
因子分析是研究如何以最少的信息丢失，将众多原始变量浓缩成少数几个因子变量，以及如何使因子变量具有较强的可解释性的一种多元统计分析方法。我们以下为例：
为了了解学生的学习能力，观测了n个学生p个科目的成绩，用X1,...,Xp表示p个科目（例如代数，几何，语文，英语......）。我们对这些资料进行归纳分析，得出全部科目X所共有的因子有m(m<p)个，如数学推导因子，记忆因子，计算因子等，分别记为F1,...Fm，即
X(i)=a(i1)F1+a(i2)F2+...+a(im)Fm+ε(i) (i=1,...,p)
用这m个不可观测的互不相关的公共因子F1...Fm和一个特殊因子ε(i)来描述原始可测的相关变量（科目）X1...Xp，并解释分析学生的学习能力。它们的系数a(i1),...a(im)称为因子载荷。这就是一个因子分析模型，即达到了降维又可以用于分类。
综上所述，因子分析是寻找潜在的起支配作用的因子模型的方法。因子分析是根据相关性大小把变量分组，使得同组内的变量之间相关性较高，但不同的组的变量相关性较低，每组变量代表一个基本结构，这个基本结构称为公共因子。对于所研究的问题就可试图用最少个数的不可测的所谓公共因子的线性函数与特殊因子之和来描述原来观测的每一分量。
通过因子分析得来的新变量是对每个原始变量进行内部剖析。因子分析不是对原始变量的重新组合，而是对原始变量进行分解，分解为公共因子和特殊因子两部分。具体地说，就是要找出某个问题中可直接测量的具有一定相关性的诸指标，如何受少数几个在专业中有意义、又不可直接测量到、且相对独立的因子支配的规律，从而可用各指标的测定来间接确定各因子的状态。因子分析只能解释部分变异，主成分分析能解释所有变异。
因子分析法是从研究变量内部相关的依赖关系出发，把一些具有错综复杂关系的变量归结为少数几个综合因子的一种多变量统计分析方法。它的基本思想是将观测变量进行分类，将相关性较高，即联系比较紧密的分在同一类中，而不同类变量之间的相关性则较低，那么每一类变量实际上就代表了一个基本结构，即公共因子。对于所研究的问题就是试图用最少个数的不可测的所谓公共因子的线性函数与特殊因子之和来描述原来观测的每一分量。
因子分析模型描述如下：
⑴X = (x1，x2，…，xp）￠是可观测随机向量，均值向量E(X)=0，协方差阵Cov(X)=∑，且协方差阵∑与相关矩阵R相等（只要将变量标准化即可实现）。
⑵F = (F1，F2，…，Fm）￠ （m<p）是不可测的向量，其均值向量E(F)=0，协方差矩阵Cov(F) =I，即向量的各分量是相互独立的。
⑶e = (e1，e2，…，ep）￠与F相互独立，且E(e)=0， e的协方差阵∑是对角阵，即各分量e之间是相互独立的，则模型：
x1 = a11F1+ a12F2 +…+a1mFm + e1
x2 = a21F1+a22F2 +…+a2mFm + e2
………
xp = ap1F1+ ap2F2 +…+apmFm + ep
称为因子分析模型，由于该模型是针对变量进行的，各因子又是正交的，所以也称为R型正交因子模型。
其矩阵形式为：x =AF + e
这里，
⑴m ￡ p；
⑵Cov(F，e)=0，即F和e是不相关的；
⑶D(F) = Im ，即F1，F2，…，Fm不相关且方差均为1；
D(e)=，即e1，e2，…，ep不相关，且方差不同。
我们把F称为X的公共因子或潜因子，矩阵A称为因子载荷矩阵，e 称为X的特殊因子。A = (aij），aij为因子载荷。数学上可以证明，因子载荷aij就是第i变量与第j因子的相关系数，反映了第i变量在第j因子上的重要性。
模型中F1，F2，…，Fm叫做主因子或公共因子，它们是在各个原观测变量的表达式中都共同出现的因子，是相互独立的不可观测的理论变量。公共因子的含义，必须结合具体问题的实际意义而定。e1，e2，…，ep叫做特殊因子，是向量x的分量xi(i=1，2，…，p）所特有的因子，各特殊因子之间以及特殊因子与所有公共因子之间都是相互独立的。
模型中载荷矩阵A中的元素（aij）是为因子载荷。因子载荷aij是xi与Fj的协方差，也是xi与Fj的相关系数，它表示xi依赖Fj的程度。可将aij看作第i个变量在第j公共因子上的权，aij的绝对值越大（|aij|￡1），表明xi与Fj的相依程度越大，或称公共因子Fj对于xi的载荷量越大。为了得到因子分析结果的经济解释，因子载荷矩阵A中有两个统计量十分重要，即变量共同度和公共因子的方差贡献。
因子载荷矩阵A中第i行元素之平方和记为hi2，称为变量xi的共同度。
它是全部公共因子对xi的方差所做出的贡献，反映了全部公共因子对变量xi的影响。hi2大表明x的第i个分量xi对于F的每一分量F1，F2，…，Fm的共同依赖程度大。将因子载荷矩阵A的第j列（ j =1，2，…，m）的各元素的平方和记为gj2，称为公共因子Fj对x的方差贡献。
gj2就表示第j个公共因子Fj对于x的每一分量xi(i= 1，2，…，p）所提供方差的总和，它是衡量公共因子相对重要性的指标。gj2越大，表明公共因子Fj对x的贡献越大，或者说对x的影响和作用就越大。如果将因子载荷矩阵A的所有gj2 ( j =1，2，…，m）都计算出来，使其按照大小排序，就可以依此提炼出最有影响力的公共因子。
因子旋转
建立因子分析模型的目的不仅是找出主因子，更重要的是知道每个主因子的意义，以便对实际问题进行分析。如果求出主因子解后，各个主因子的典型代表变量不很突出，还需要进行因子旋转，通过适当的旋转得到比较满意的主因子。
旋转的方法有很多，正交旋转（orthogonal rotation）和斜交旋转（oblique rotation）是因子旋转的两类方法。最常用的方法是最大方差正交旋转法（Varimax）。进行因子旋转，就是要使因子载荷矩阵中因子载荷的平方值向0和1两个方向分化，使大的载荷更大，小的载荷更小。因子旋转过程中，如果因子对应轴相互正交，则称为正交旋转；如果因子对应轴相互间不是正交的，则称为斜交旋转。常用的斜交旋转方法有Promax法等。
因子得分
因子分析模型建立后，还有一个重要的作用是应用因子分析模型去评价每个样品在整个模型中的地位，即进行综合评价。例如地区经济发展的因子分析模型建立后，我们希望知道每个地区经济发展的情况，把区域经济划分归类，哪些地区发展较快，哪些中等发达，哪些较慢等。这时需要将公共因子用变量的线性组合来表示，也即由地区经济的各项指标值来估计它的因子得分。
设公共因子F由变量x表示的线性组合为：
Fj = uj1 xj1+ uj2 xj2+…+ujpxjp j=1，2，…，m
该式称为因子得分函数，由它来计算每个样品的公共因子得分。若取m=2，则将每个样品的p个变量代入上式即可算出每个样品的因子得分F1和F2，并将其在平面上做因子得分散点图，进而对样品进行分类或对原始数据进行更深入的研究。
但因子得分函数中方程的个数m小于变量的个数p，所以并不能精确计算出因子得分，只能对因子得分进行估计。估计因子得分的方法较多，常用的有回归估计法，Bartlett估计法，Thomson估计法。
因子分析的核心问题有两个：一是如何构造因子变量；二是如何对因子变量进行命名解释。因此，因子分析的基本步骤和解决思路就是围绕这两个核心问题展开的。
（i）因子分析常常有以下四个基本步骤：
⑴确认待分析的原变量是否适合作因子分析。
⑵构造因子变量。
⑶利用旋转方法使因子变量更具有可解释性。
⑷计算因子变量得分。
（ii）因子分析的计算过程：
⑴将原始数据标准化，以消除变量间在数量级和量纲上的不同。
⑵求标准化数据的相关矩阵；
⑶求相关矩阵的特征值和特征向量；
⑷计算方差贡献率与累积方差贡献率；
⑸确定因子：
设F1，F2，…， Fp为p个因子，其中前m个因子包含的数据信息总量（即其累积贡献率）不低于80%时，可取前m个因子来反映原评价指标；
⑹因子旋转：
若所得的m个因子无法确定或其实际意义不是很明显，这时需将因子进行旋转以获得较为明显的实际含义。
⑺用原指标的线性组合来求各因子得分：
采用回归估计法，Bartlett估计法或Thomson估计法计算因子得分。
⑻综合得分
以各因子的方差贡献率为权，由各因子的线性组合得到综合评价指标函数。
F = (w1F1+w2F2+…+wmFm)/(w1+w2+…+wm )
此处wi为旋转前或旋转后因子的方差贡献率。
⑼得分排序：利用综合得分可以得到得分名次。
因子分析Stata操作
读取数据：
数据展示：
因子分析：
因子负荷：
Uniqueness：代表测量误差和未能被因子所解释的那部分，因此当它的值大于0.6，我们一般认为，这些因子未能很好地解释各个原始变量。这个示例里面，恰恰都出现了Uniqueness>0.6的情形，所以这里的三个因子并不能很好地解释原始变量。
主成分分析与因子分析不同点
1、因子分析中是把变量表示成各因子的线性组合，而主成分分析中则是把主成分表示成各变量的线性组合。
2、主成分分析的重点在于解释各变量的总方差，而因子分析则把重点放在解释各变量之间的协方差。
3、主成分分析中不需要有假设（assumptions），因子分析则需要一些假设。因子分析的假设包括：各个共同因子之间不相关，特殊因子（specific factor）之间也不相关，共同因子和特殊因子之间也不相关。
4、主成分分析中，当给定的协方差矩阵或者相关矩阵的特征值是唯一的时候，主成分一般是独特的；而因子分析中因子不是独特的，可以旋转得到不同的因子。
5、在因子分析中，因子个数需要分析者指定（spss根据一定的条件自动设定，只要是特征值大于1的因子进入分析），而指定的因子数量不同而结果不同。在主成分分析中，成分的数量是一定的，一般有几个变量就有几个主成分。和主成分分析相比，由于因子分析可以使用旋转技术帮助解释因子，在解释方面更加有优势。
大致说来，当需要寻找潜在的因子，并对这些因子进行解释的时候，更加倾向于使用因子分析，并且借助旋转技术帮助更好解释。而如果想把现有的变量变成少数几个新的变量（新的变量几乎带有原来所有变量的信息）来进入后续的分析，则可以使用主成分分析。当然，这种情况也可以使用因子得分做到。所以这种区分不是绝对的。
在算法上，主成分分析和因子分析很类似，不过在因子分析中所采用的协方差矩阵的对角元素不再是变量的方差，而是和变量对应的共同度（变量方差中被各因子所解释的部分）。
Suggested readings
动态面板回归和软件操作，单位根和协整检验（Dynamic Panel Data）
政策效应评估的四种方法（Policy evaluation）
倾向匹配分析深度（Propsensity matching analysis）
写在后面：各位圈友，咱们的计量经济圈社群里面资料和计量咨询都很多，希望大家能够积极加入咱们这个大家庭（戳这里）。我们已经邀请社群里的圈友建立了微信群与圈圈对话，进去之后一定要看“群公告”，不然接收不了群息。
帮点击一下下面的小广告，谢谢支持！