1 利用Excel2000进行主成分分析
第一步,录入数据,并对进行标准化。
【例】一组古生物腕足动物贝壳标本的两个变量:长度和宽度。
图1 原始数据和标准化数据及其均值、方差 (取自张超、杨秉庚《计量地理学基础》)
计算的详细过程如下:
⑴ 将原始数据绘成散点图(图2)。主持分分析原则上要求数据具有线性相关趋势——如果数据之间不相关(即正交),则没有必要进行主成分分析,因为主成分分析的目的就是用正交的变量代替原来非正交的变量;如果原始数据之间为非线性关系,则有必要对数据进行线性转换,否则效果不佳。从图2 可见,原始数据具有线性相关趋势,且测定系数R2=0.4979,相应地,相关系数R=0.7056。
⑵ 对数据进行标准化。标准化的数学公式为
1
*xijxijxjj
这里假定按列标准化,式中
n1nxijxij,ij(xijxj)2Var(xij)
ni1i1分别为第j列数据的均值和标准差,xij为第i行(即第i个样本)、第j列(即第j个变
*量)的数据,xij为相应于xij的标准化数据,n25为样本数目。
原始数据的散点图2520y = 0.7686x + 2.31742R = 0.4979宽度1510500510长度15图2 原始数据的散点图
2025
标准化数据的散点图y = 0.7056x + 2E-162R = 0.49793.0000002.000000宽度1.0000000.000000-2.00000-1.000000.0000001.0000002.0000003.000000-1.00000000-2.000000长度图3 标准化数据的散点图
对数据标准化的具体步骤如下:① 求出各列数据的均值,命令为average,语法为:
2
average(起始单元格:终止单元格)。如图1所示,在单元格B27中输入“=AVERAGE(B1:B26)”,确定或回车,即得第一列数据的均值x110.88;然后抓住单元格B27的右下角(光标的十字变细)右拖至C27,便可自动生成第二列数据的均值x210.68。
②求各列数据的方差。命令为varp,语法同均值。如图1所示,在单元格B28中输入“=VARP(B2:B26)”,确定或回车,可得第一列数据的方差Var(x1)19.4656,右拖至C28生成第二列数据的方差Var(x2)23.0976。
③ 求各列数据的标准差。将方差开方便得标准差。也可利用命令stdevp直接生成标准差,语法和操作方法同均值、方差,不赘述。
④ 标准化计算。如图1所示,在单元格D2中输入“=(B2-$B$27)/$B$29”,回车可得第一列第一个数据“3”的标准化数值-1.786045,然后按住单元格D2的右下角下拖至D26,便会生成第一列数据的全部标准化数值;按照单元格D2的右下角右拖至E2,就能生成第二列第一个数据“2”的标准化数据-1.806077,抓住单元格E2的右下角下拖至E26便会生成第二列数据的全部标准化数值。
⑤ 作标准化数据的散点图(图3)。可以看出,点列的总体趋势没有变换,两种数据的相关系数与标准化以前完全相同。但回归模型的截距近似为0,即有a0,斜率等于相关系数,即有bR。
⑶ 求标准化数据的相关系数矩阵或协方差矩阵。求相关系数矩阵的方法是:沿着“工具(T)”→“数据分析(D)”的路径打开“分析工具(A)”选项框(图4),确定,弹出“相关系数”对话框(图5),在“输入区域”的空白栏中输入标准化数据范围,并以单元格G1为输出区域,具体操作方法类似于回归分析。确定,即会在输出区域给出相关
图4 分析工具选项框
图5 相关系数对话框
3
系数矩阵的下三角即对角线部分,由于系对称矩阵,上三角的数值与下三角相等,故未给出(图6),可以通过“拷贝——转置——粘帖”的方式补充空白部分。
图6 标准化数据的相关系数和协方差
求协方差的方法是在“分析工具”选项框中选择“协方差”(图7),弹出“协方差” 选项框(图8),具体设置与“相关系数”类似,不赘述。结果见图6,可以看出,对于标准化数据而言,协方差矩阵与相关系数矩阵完全一样。因此,二者任取其一即可。
图7 在分析工具选项框中选择“协方差”
图8 协方差选项框
⑷ 计算特征根。我们已经得到相关系数矩阵为
0.70561C,
0.70561而二阶单位矩阵为
10I,
01 4
于是根据公式det(IC)0,我们有
100110.70560.7056110.70560.705610
按照行列式化为代数式的规则可得
(1)20.70562220.50210
根据一元二次方程的求根公式,当b24ac0时,我们有
bb24ac
2a据此解得11.7056,20.2944(对于本例,显然11R,21R)。这便是
相关系数矩阵的两个特征根。
⑸ 求标准正交向量。将1代入矩阵方程(IC)0,得到
0.70560.7056100.70560.70560 2在系数矩阵IC中,用第一行加第二行,化为
0.70560.705610 0020由此得12,令11,则有21,于是得基础解系
1,单位化为e110.7071单位化的公式为ei10.7071i1222(i1,2)。
完全类似,将2代入矩阵方程(IC)0,得到
0.70560.7056100.70560.70560 2用系数矩阵的第二行减去第一行,化为
0.70560.70561000
02于是得到12,取11,则有21,因此得基础解系为
10.7071 2,单位化为e210.7071这里e1、e2便是标准正交向量。
⑹ 求对角阵。首先建立标准正交矩阵P,即有
P[e10.70710.7071 e2]0.70710.7071 5
该矩阵的一个特殊性质便是PTP1,即矩阵的转置等于矩阵的逆。根据DPTCP,可知
0.70560.70710.70711.705600.70710.70711D0.70560.70710.70710
0.70710.707110.2934下面说明一下利用Excel进行矩阵乘法运算的方法。矩阵乘法的命令为mmult,语法是
mmult (矩阵1的单元格范围,矩阵2的单元格范围)。例如,用矩阵PT 与矩阵C相乘,首先选择一个输出区域如G1:H2,然后输入“=mmult(A1:B2,C1:D2)”,然后按下“Ctrl+Shift+Enter”键(图9),即可给出
1.206044 1.206044 0.20817 -0.20817
再用乘得的结果与P阵相乘,便得对角矩阵
1.705603 0
0 0.294397
如果希望一步到位也不难,选定输出区域如C3:D4,然后输入“=mmult(mmult(A1:B2,C1:D2),E1:F2)” (图10),同时按下“Ctrl+Shift+Enter”键,立即得到结果(图11)。显然,对角矩阵对角线的数值恰是相关系数矩阵的特征值。
图9 矩阵乘法示例
图10 矩阵连乘的命令与语法
至此,标准化的原始变量x与主成分之间z之间可以表作
x10.7056x11x2z11x20.70560z11.7056z2 0.2944z20 6
显然z1与z2之间正交。
图11 乘法结果:对角矩阵
⑺ 根据特征根计算累计方差贡献率。现已求得第一特征根为11.7056,第二特征根为
20.2944,二者之和刚好就是矩阵的维数,即有12m2,这里m=2为变量数
目(注意前面的n=25为样本数目)。比较图6或图10中给出的相关系数矩阵C与图11中给出的对角矩阵D可以看出,Tr.(C)=1+1=2,Tr.(D)=1.7056+0.2944=2,即有Tr.(C)= Tr.(D),可见将相关系数亦即协方差矩阵转换为对角矩阵以后,矩阵的迹(trace,即对角线元素之和)没有改变,这意味着将原始变量化为主成分以后,系统的信息量没有减少。现在问题是,如果我们只取一个主成分代表原来的两个变量,能反映原始变量的多少信息?这个问题可以借助相关系数矩阵的特征根来判断。利用Excel容易算出,第一特征根占特征根总和即矩阵维数的85.28%(见下表),即有
特征根 累计值 百分比 累计百分比 1.705603 1.705603 85.28% 85.28% 0.294397 2 14.72% 100.00%
也就是说:
1:1.7056,1/m1.7056/285.28% 2:0.2944,2/m0.2944/m14.72%
12:2,(12)/m2/2100%
这表明,如果仅取第一个主成分,可以反映原来数据85.28%的信息——换言之,舍弃第二个主成分,原来数据的信息仅仅损失14.72%,但分析变量的自由度却减少一个,整个分析将会显得更加简明。
⑻ 计算主成分载荷。根据公式jjej,容易算出
0.70710.9235 11.70560.9235
0.7071
20.29440.70710.38370.3837
0.7071
⑼ 计算公因子方差和方差贡献。根据上述计算结果可以比较公因子方差和方差贡献。再考虑全部的两个主成分的时候,对应于1和2的公因子方差分别为 V1j2ij0.923520.383721
7
V2ij0.92352(0.3837)21
2j对应于第一主成分z1和第二主成分z2的方差贡献分别为
ij0.923520.923521.7056 CV1iiCV2ij0.38372(0.3837)20.2944
可以看出(图12): 第一,方差贡献等于对应主成分的特征根,即有
CVjj 第二,公因子方差相等或彼此接近,即有
V1V2
第一,公因子方差之和等于方差贡献之和,即有
VCViijjm2
第一个规律是我们决定提取主成分数目的判据与之一,第二个规律是我们判断提取主成分数目是否合适的判据之一,第三个规律是我们判断提取主成分后是否损失信息的判据之一。去掉次要的主成分以后,上述规律理当仍然满足。这时如果第二个规律不满足,就意味着主成分的提取是不合适的。此外,上述规律也是我们检验计算结果是否正确的判据之一。
图12 公因子方差、方差贡献的计算结果及其与特征根的贡献
⑽ 计算主成分得分。根据主成分与原始变量的关系,应有
ZPTX
或者
XPZ
对于本例而言,式中
x1z1e11e120.70710.7071X,Z,Pe1e20.70710.7071
eexz212222TT这里e1e11e12,e2e21e22为前面计算的标准化特征向量。于是有
z10.70710.7071x1z0.70710.7071x
22化为代数形式便是
z10.7071x10.7071x2
z20.7071x10.7071x2
式中的x均为标准化数据。对ZPTX进行转置,可得
8
ZTXTP
图13 计算特征向量的公式及语法
图14 计算主成分得分
根据这个式子,利用Excel计算主成分得分的步骤如下:
9
① 将特征向量复制到标准化数据的附近;
② 选中一个与标准化数据占据范围一样大小的数值区域(如G2:H26);
③ 输入如下计算公式“=mmult(标准化数据的范围,特征向量的范围)”,在本例中就是“=MMULT(B2:C26,E2:F3)”(图13); ④ 同时按下“Ctrl+Shift+Enter”键。
⑤ 计算主成分得分的均值和方差,可以发现,均值为0(由于误差之故,约等于0),方差等于特征根。
⑥ 最后,可以对主成分得分进行标准化。已知主成分得分的均值为0,我们不按总体方差进行标准化,而按样本方差进行标准化。
图15 主成分得分的标准化结果
样本方差的计算公式为
1nVar(xj)(xijxj)2 n1i1相应地,标准差为
10
jVar(xj)1n(xijxj)2 n1i1标准化公式同前面给出的一样。结果见表15。注意,这里之所以按样本方差进行标准化,
主要目的是为了与SPSS的计算结果进行比较。
分别以z1、z2为坐标轴,将主成分得分(包括标准化的得分)点列标绘于坐标图中,可以发现,点列分布没有任何趋势:回归结果表明,回归系数和相关系数均为零,即有a0,b0,R0(图16,图17)。这从几何图形上显示:主成分之间是正交的,即有cos0(试将图16、图17与图2、图3对比)。
主成分得分的空间分布1.5000001.0000000.500000第二主成分得分y = -7E-17x - 2E-162R = 2E-320.000000-3.000000-2.000000-1.0000000.0000001.0000002.0000003.0000004.000000-0.500000-1.000000-1.500000第一主成分得分图16 主成分得分的相关系数为零
主成分得分的空间分布(标准化)2.52第二主成分得分1.510.50-3-2-1-0.50-1-1.5-2y = -2E-16x - 4E-172R = 3E-32123第一主成分得分图17 主成分得分的相关系数为零(标准化)
最后可以验证因子载荷即为(标准化)原始数据与主成分得分之间的相关系数,容易算出
11
(x1,z1)Correl(x1,z1)0.9235, (x2,z1)Correl(x2,z1)0.9235, (x1,z2)Correl(x1,z2)0.3837, (x2,z2)Correl(x2,z2)0.3837
图表标题x1-z143210-2-1-10-2-3图18 x1与z1的关系及其回归方程
z1 = 1.206x1 + 3E-17R = 0.85282z1得分线性 (z1得分)123
图表标题x2-z143210-2-1-10-2-3图19 x2与z1的关系及其回归方程
z1 = 1.206x2 - 3E-16R = 0.8528z1得分线性 (z1得分)1232
12
图表标题x1-z21.510.50-2-1-0.5-1-1.5图20 x1与z2的关系及其回归方程
z2 = 0.2082x1 - 1E-16R = 0.147220123z2得分线性 (z2得分)
图表标题x2-z21.510.50-2-1-0.5-1-1.5图21 x2与z2的关系及其回归方程
z2 = -0.2082x2 - 2E-16R = 0.147220123z2得分线性 (z2得分)
回归方程为
z11.206x1 z11.206x2 z20.2082x1 z20.2082x2
方程的系数恰是以下矩阵的元素
0.70561.206041.206040.70710.70711GPC0.70560.208170.20817
0.70710.70711 13