第六章 主成分分析与因子分析 6.1 主成分分析 6.2 因子分析
6.1 主成分分析 6.1.1 主成分分析的概念与步骤 6.1.2 使用INSIGHT模块作主成分分析 6.1.1 主成分分析的概念与步骤 6.1.2 使用INSIGHT模块作主成分分析 6.1.3 使用“分析家”作主成分分析 6.1.4 使用PRINCOMP过程进行主成分分析
6.1.1 主成分分析的概念与步骤 1. 主成分分析基本思想 6.1.1 主成分分析的概念与步骤 1. 主成分分析基本思想 主成分分析是数学上对数据降维的一种方法。其基本思想是设法将原来众多的具有一定相关性的指标(比如p个指标),重新组合成一组新的互不相关的综合指标来代替原来指标。通常数学上的处理就是将原来p个指标作线性组合,作为新的综合指标。但是这种线性组合,如果不加限制,则可以有很多,应该如何去选取呢?
在所有的线性组合中所选取的F1应该是方差最大的,故称F1为第一主成分。如果第一主成分不足以代表原来p个指标的信息,再考虑选取F2即选第二个线性组合。为了有效地反映原有信息,F1已有的信息就不需要再出现在F2中,用数学语言表达就是要求Cov(F1,F2)=0。称F2为第二主成分,依此类推可以构造出第三、第四、…、第p个主成分。
2. 主成分分析的数学模型 设有n个样品(多元观测值),每个样品观测p项指标(变量):X1,X2,…,Xp,得到原始数据资料阵: 其中Xi = (x1i,x2i,…,xni)',i = 1,2,…,p。
Fi = a1iX1 + ai2X2 +…+apiXp i = 1,2,…,p 用数据矩阵X的p个列向量(即p个指标向量)X1,X2,…,Xp作线性组合,得综合指标向量: 简写成: Fi = a1iX1 + ai2X2 +…+apiXp i = 1,2,…,p
为了加以限制,对组合系数ai' = (a1i,a2i,…,api)作如下要求: 1) Fi与Fj(ij, i, j = 1, …, p)互不相关,即Cov(Fi,Fj) = ai'ai = 0,其中Σ是X的协方差阵。 2) F1是X1,X2,…,Xp的一切线性组合(系数满足上述要求)中方差最大的,即 ,其中c = (c1,c2,…,cp)' F2是与F1不相关的X1,X2,…,Xp一切线性组合中方差最大的,…,Fp是与F1,F2,…,Fp-1都不相关的X1,X2,…,Xp的一切线性组合中方差最大的。
满足上述要求的综合指标向量F1,F2,…,Fp就是主成分,这p个主成分从原始指标所提供的信息总量中所提取的信息量依次递减,每一个主成分所提取的信息量用方差来度量,主成分方差的贡献就等于原指标相关系数矩阵相应的特征值i,每一个主成分的组合系数 ai' = (a1i,a2i,…,api) 就是相应特征值i所对应的单位特征向量ti。方差的贡献率为 ,i越大,说明相应的主成分反映综合信息的能力越强。
3. 主成分分析的步骤 (1) 计算协方差矩阵 计算样品数据的协方差矩阵:Σ = (sij)pp,其中 i,j = 1,2,…,p (2) 求出Σ的特征值及相应的特征向量 求出协方差矩阵Σ的特征值12…p>0及相应的正交化单位特征向量: 则X的第i个主成分为Fi = ai'X i = 1,2,…,p。
(3) 选择主成分 在已确定的全部p个主成分中合理选择m个来实现最终的评价分析。一般用方差贡献率 解释主成分Fi所反映的信息量的大小,m的确定以累计贡献率 达到足够大(一般在85%以上)为原则。
(4) 计算主成分得分 计算n个样品在m个主成分上的得分: i = 1,2,…,m (5) 标准化 实际应用时,指标的量纲往往不同,所以在主成分计算之前应先消除量纲的影响。消除数据的量纲有很多方法,常用方法是将原始数据标准化,即做如下数据变换: 其中 , ,j = 1,2,…,p。标准化后的数据阵记为X*,其中每个列向量(标准化变量)的均值为0,标准差为1,数据无量纲。
Fj = a1jX1* + a2jX2* +...+ apjXp* j = 1,2,…,m 标准化后变量的协方差矩阵(Covariance Matrix)Σ = (sij)pp,即原变量的相关系数矩阵(Correlation Matrix)R= (rij)pp: i,j = 1,2,…,p 此时n个样品在m个主成分上的得分应为: Fj = a1jX1* + a2jX2* +...+ apjXp* j = 1,2,…,m
6.1.2 使用INSIGHT模块作主成分分析 【例6-1】全国沿海10个省市经济指标的主成分分析 表6-1 全国沿海10个省市经济综合指标 假设表6-1中数据已经存放在数据集Mylib.jjzb中,试对各地区的经济发展水平进行主成分分析。 地区 GDPx1 人均GDPx2 工业增加值x3 第三产业增加值x4 固定资产投资x5 基本建设投资x6 社会消费品零售总额x7 海关出口总额x8 地方财政收入x9 辽宁 5458.2 13000 1376.2 2258.4 1315.9 529 123.7 399.7 山东 10550 11643 3502.5 3851 2288.7 1070.7 3181.9 211.1 610.2 河北 6076.6 9047 1406.7 2092.6 1161.6 597.1 1968.3 45.9 302.3 天津 2022.6 22068 822.8 960 703.7 361.9 941.4 115.7 171.8 江苏 10636.3 14397 3536.3 3967.2 2320 1141.3 3215.8 384.7 643.7 上海 5408.8 40627 2196.2 2755.8 1970.2 779.3 2035.2 320.5 709 浙江 7670 16570 2356.5 3065 2296.6 1180.6 2877.5 294.2 566.9 福建 4682 13510 1047.1 1859 964.5 397.9 1663.3 173.7 272.9 广东 11769.7 15030 4224.6 4793.6 3022.9 1275.5 5013.6 1843.7 1201.6 广西 2455.4 5062 367 995.7 542.2 352.7 1025.5 15.1 186.7
1. 使用INSIGHT模块做主成分分析的步骤 1) 在INSIGHT模块中打开数据集Mylib.jjzb;选择菜单“Analyze”“Multivariate(Y X)(多元分析)”,打开“Multivariate(Y X)”对话框; 2) 将做主成分分析的变量x1~x9选为Y变量,将变量diqu选为Label变量,如图所示。
图6-1 多元分析对话框 3) 单击“Method”按钮,在打开的对话框中可以选择计算协方差矩阵的特征值或是计算相关系数矩阵的特征值。系统默认计算相关系数矩阵的特征值和特征向量,单击“OK”按钮返回。
4) 单击“Output”按钮,在打开的对话框(图左)中包括“Descriptive Statistics”选项、“Bivariate Plots”选项以及各种多元分析的选项。选中“Principal Component Analysis”复选框,单击下面的“Principal Component Options”按钮,打开“Principal Component Options”对话框,选中“Eigenvectors”复选框,取消“Correlations(Structure)”复选框,如图右所示。
2. 主成分的结果分析 输出的数字分析结果有4个部分:简单统计量、相关系数矩阵、相关系数矩阵的特征值以及相关系数矩阵的特征向量。
3) 图6-5给出相关系数矩阵的特征值(Eigenvalue)、上下特征值之差(Difference)、各主成分的方差贡献率(Proportion)以及累积贡献率(Cumulative)。 相关系数矩阵的特征值即各主成分的方差,可以看出,第一主成分的方差贡献率为80.11%,前两个主成分的累积贡献率已达92.33%,因此,只需用前面2个主成分就可以概括这组数据。
4) 图6-6给出相关系数矩阵的两个最大特征值的特征向量,据此可以写出第一和第二主成分得分: PCR1 = 0.35x1* + 0.04x2* + 0.36x3* + 0.37x4* + 0.37x5* + 0.35x6* + 0.36x7* + 0.30x8* + 0.36x9* PCR2 = -0.21x1* + 0.94x2* – 0.01x3* – 0.05x4* + 0.10x5* – 0.02x6* – 0.14x7* + 0.05x8* + 0.18x9* 对于第一主成分而言,除了x2(人均GDP)外,各变量所占比重均在0.3左右以上,因此第一主成分(Prin1)主要由x1、x3~x9八个变量解释;而第二主成分则主要由x2这一个变量解释。
5) 选择菜单“Edit(编辑)”“Observations(观测)”“Label in Plots”,在弹出的对话框中选中所有diqu变量值,单击“OK”按钮返回,显示结果中的散点图上出现地区名; 图中看出,上海在第二主成分PCR2的得分远远高于其他省市,而在第一主成分PCR1的得分则处于中间。广东、江苏、山东和浙江则在第1主成分的得分上位于前列。
6) 回到INSIGHT的数据窗口,可以看到前两个主成分的得分情况(如图6-8左)。 单击数据窗口左上角的箭头,在弹出的菜单中选择“Sort(排序)”选项,在打开的对话框中选定排序变量PCR1,并单击“Asc/Des”按钮将其设为降序(Des),如图6-8所示。
单击“OK”按钮返回,得到按第一主成分排序的结果如图6-9左所示。同样方法可以得到按第二主成分排序的结果如图6-9右所示。 从第一主成分排序情况来看,沿海19省市经济发展状况综合排名前5位的省市依次为:广东、江苏、山东、浙江、上海;从第二主成分排序情况来看,人均GDP排名前5位的省市依次是:上海、天津、浙江、广东、福建。
6.1.3 使用“分析家”作主成分分析 【例6-2】某企业为了了解其客户的信用程度,评价客户的信用等级,采用信用评估常用的5C方法,5C的目的是说明顾客违约的可能性。 1) 品格x1,指客户的信誉。 2) 能力x2,指客户的偿还能力。 3) 资本x3,指客户的财务势力和财务状况。 4) 附带的担保品x4。 5) 环境条件x5,指客户的外部因素。
通过专家打分,得到10个客户5项指标的得分如表6-3所示。 表6-2 10个客户5项指标的得分 假设表6-2中数据已经存放在数据集Mylib.xydj中,试对各客户的信用等级进行评估。 客户编号ID x1 x2 x3 x4 x5 1 76.5 81.5 76 75.8 71.7 6 85 79.2 80.3 84.4 2 70.6 73 67.6 68.1 78.5 7 94 87.5 89.5 92 3 90.7 87.3 91 80 8 84.6 66.9 68.8 64.8 66.4 4 77.5 73.6 70.9 69.8 74.8 9 57.7 60.4 57.4 60.8 65 5 85.6 68.5 70 62.2 10 69.2 64.9 68.9
1. 使用“分析家”做主成分分析的步骤 1) 在“分析家”中打开数据集Mylib.xydj; 2) 选择菜单“Statistics(统计)”“Multivariate(多元分析)”“Principal Components(主成分分析)”,打开“Principal Components”对话框; 3) 在对话框中输入主成分分析的变量,如图所示。
4) 单击“Statistics(统计)”按钮,打开“Principal Components:Statistics”对话框; 在“# of components:”右边的框中指定主成分的个数4,如图右。单击“OK”返回;
5) 单击“Save Data”按钮,打开“Principal Components:Save Data”对话框,在该对话框中可选择存储数据。 选中“Create and save scores data”,如图6-11所示。单击“OK”返回;
6) 单击“Plots”按钮,打开“Principal Components:Plots”对话框,可以设置图形输出。 ● 在“Scree Plot (碎石图)”选项卡中(图左),选中“Create scree plot(建立碎石图)”复选框。 ● 在“Component Plot (成分图)”选项卡中(图右),选中“Create component Plot(建立成分图)”复选框。
2. 主成分的结果分析 输出的数字分析结果包括4个部分:简单统计量、相关系数矩阵、相关系数矩阵的特征值以及相关系数矩阵的特征向量。 1) 图6-13给出变量的简单统计量,图中显示5项指标中品格、能力和附带担保品是最为重要的,其标准差高出其他变量。
2) 图6-14给出各变量之间的相关系数矩阵。可以看出,能力与资本、附带担保品有着较强的相关性,表明客户的偿还能力与其财务实力、财务状况和抵押资产有着重要的关系。
3) 图6-15给出相关系数矩阵的特征值(Eigenvalues)、上下特征值之差(Difference)、各主成分的方差贡献率(proportion)以及累积贡献率(Cumulative)。 相关系数矩阵的特征值即各主成分的方差,可以看出,第一主成分的方差贡献率为84.22%,第二主成分的方差贡献率为7.67%,第三主成分的方差贡献率为5.95%。说明第一主成分已经具有足够多的方差贡献率,可以很好地概括这组数据。
在“分析家”左边的管理窗口中双击“Scree plot”项,打开的“Scree plot”对话框显示前4个特征值的“碎石图”,很直观地看到第一主成分远远大于其它特征值,说明第一主成分已经代表了绝大部分信息。
4) 图6-16给出相关系数矩阵的特征向量,由最大特征值所对应的特征向量可以写出第一主成分的表达式。 Prin1 = 0.4135x1* + 0.4729x2* + 0.4656x3* + 0.4547x4* + 0.4265x5* 利用特征向量各分量的值可以对主成分进行解释,对于第一主成分而言,各变量所占比重大致相等,且均为正数,说明第一主成份是对所有指标的一个综合测度,作为综合的信用等级指标,可以用来排序。
5) 在“分析家”窗口中,双击左边项目管理中的“Scores Table”项,打开“Scores Table”对话框; 选择菜单“File”“Save as By SAS Name”,将其保存为数据表Scores;然后,在VIEWTABLE中打开该表;选择菜单“Data”“Sort”,按主成分Prin1排序,结果如表6-3所示。 表6-3 客户的信用等级 在正确评估了顾客的信用等级后,就能正确制定出对其的信用期、收账政策等,这对于加强应收账款的管理大有帮助。 客户编号 1 2 3 4 5 6 7 8 9 10 第一主成分得分 3.17 -9.01 25.09 -4.36 -6.41 13.62 35.88 -10.34 -33.80 -13.83 名次
6.1.4 使用PRINCOMP过程进行主成分分析 1. PRINCOMP过程的功能简介 由特征向量得出相应的主成分,用少数几个主成分代替原始变量,并计算主成分得分。 2) 主成分的个数可以由用户自己确定,主成分的名字可以用户自己规定,主成分得分是否标准化可由用户规定。
3) 输入数据集可以是原始数据集、相关阵、协方差阵等。输入为原始数据时,还可以规定从协方差阵出发还是从相关阵出发进行分析,由协方差阵出发时方差大的变量在分析中起到更大的作用。 4) 该过程还可生成两个输出数据集:一个包含原始数据及主成分得分,它可作为主成分回归和聚类分析的输入数据集;另一个包含有关统计量,类型为TYPE = CORR或COV的输出集,它也可作为其他过程的输入SAS集。
2. PRINCOMP过程的格式 PRINCOMP过程的常用格式如下: PROC PRINCOMP <选项列表>; VAR 变量列表; [WEIGHT 变量列表;] [FREQ 变量列表;] [PARTIAL 变量列表;] [BY 变量列表;] RUN;
1) PROC PRINCOMP语句用来规定输入输出和一些运行选项,其选项及功能见表6-4。 其中: 1) PROC PRINCOMP语句用来规定输入输出和一些运行选项,其选项及功能见表6-4。 表6-4 PROC PRINCOMP语句的选项 2) VAR语句指定用于主成分分析的变量,变量必须为数值型(区间型)变量。缺省使用DATA = 输入数据集中所有数值型变量进行主成分分析。 DATA = 输入数据集,可以是原始数据集,也可以是TYPE = CORR,COV的数据集; OUT = 输出包含原始数据和主成分得分的数据集; OUTSTAT = 统计量输出数据集; COVARIANCE | COV 要求从协方差阵出发计算主成分,缺省为从相关阵出发计算。 N = 要计算的主成分个数,缺省时全部计算。 STANDARD | STD 要求在OUT = 的数据集中把主成分得分标准化为单位方差。缺省时主成分得分的方差为相应特征值。 PREFIX = 主成分名字的前缀,缺省时为PRIN1、PRIN2…。
3. 应用实例 【例6-3】对全国30个省市自治区经济发展基本情况的八项指标作主成分分析,原始数据如表6-5。 表6-5 全国30个省市自治区经济发展基本情况 省份 GDPx1 居民消费水平x2 固定资产投资x3 职工平均工资x4 货物周转量x5 居民消费价格指数x6 商品零售价格指数x7 工业总产值x8 北京 1394.89 2505 519.01 8144 373.9 117.3 112.6 843.43 天津 920.11 2720 345.46 6501 342.8 115.2 110.6 582.51 河北 2849.52 1258 704.87 4839 2033.3 115.8 1234.85 山西 1092.48 1250 290.9 4721 717.3 116.9 115.6 697.25 内蒙 832.88 1387 250.23 4134 781.7 117.5 116.8 419.39 辽宁 2793.37 2397 387.99 4911 1371.1 116.1 114 1840.55 吉林 1129.2 1872 320.45 4430 497.4 114.2 762.47 黑龙江 2014.53 2334 435.73 4145 824.8 114.3 1240.37 上海 2462.57 5343 996.48 9279 207.4 118.7 113 1642.95 江苏 5155.25 1926 1434.95 5943 1025.5 2026.64 浙江 3524.79 2249 1006.39 6619 754.4 116.6 113.5 916.59
省份 GDPx1 居民消费水平x2 固定资产投资x3 职工平均工资x4 货物周转量x5 居民消费价格指数x6 商品零售价格指数x7 工业总产值x8 安徽 2003.58 1254 474 4609 908.3 114.8 112.7 824.14 福建 2160.52 2320 553.97 5857 609.3 115.2 114.4 433.67 江西 1205.11 1182 282.84 4211 411.7 116.9 115.9 571.84 山东 5002.34 1527 1229.55 5145 1196.6 117.6 114.2 2207.69 河南 3002.74 1034 670.35 4344 1574.4 116.5 114.9 1367.92 湖北 2391.42 571.68 4685 849 120 116.6 1220.72 湖南 2195.7 1408 422.61 4797 1011.8 119 115.5 843.83 广东 5381.72 2699 1639.83 8250 656.5 114 111.6 1396.35 广西 1606.15 1314 382.59 5105 556 118.4 116.4 554.97 海南 364.17 1814 198.35 5340 232.1 113.5 111.3 64.33 四川 3534 1261 822.54 4645 902.3 118.5 117 1431.81 贵州 630.07 942 150.84 4475 301.1 121.4 117.2 324.72 云南 1206.68 334 5149 310.4 121.3 118.1 716.65 西藏 55.98 1110 17.87 7382 4.2 117.3 5.57 陕西 1000.03 1208 300.27 4396 500.9 600.98 甘肃 553.35 1007 114.81 5493 507 119.8 468.79 青海 165.31 1445 47.76 5753 61.6 118 116.3 105.8 宁夏 169.75 1355 61.98 5079 121.8 117.1 115.3 新疆 834.57 1469 376.95 5348 339 119.7 116.7 428.76
假定上述数据已经存放在数据集Mylib.jjfz中。 (2) 执行主成分分析的PRINCOMP过程 (1) 数据集 假定上述数据已经存放在数据集Mylib.jjfz中。 (2) 执行主成分分析的PRINCOMP过程 对数据集jjfz执行主成分分析的PRINCOMP过程代码如下: proc princomp data = Mylib.jjfz n = 4 out = w1 outstat = w2; var x1-x8; proc print data = w1; run;
(3) 结果分析 在各变量之间的相关系数矩阵中可以看出,有较强相关性的变量依次为: GDP(x1)与固定资产投资(x3)之间的相关系数为0.9506; GDP(x1)与工业总产值(x8)之间的相关系数为0.8737; 固定资产投资(x3)与工业总产值(x8)之间的相关系数为0.7919; 居民消费价格指数(x6)与商品零售价格指数(x7)之间的相关系数为0.7628; 货物周转量(x5)与工业总产值(x8)之间的相关系数为0.6586,等等。
图6-18给出相关系数矩阵的特征值、上下特征值之差、各主成分对方差的贡献率以及累积的贡献率。 相关系数矩阵的特征值即各主成分的方差,可以看出,第一主成分对方差的贡献率为46.94%,第二主成分对方差的贡献率为27.46%,第三主成分对方差的贡献率为15.19%,之后的主成分的贡献率为0.05。前三个主成分的累积贡献率为89.58%,因此,对第四主成分以后的主成分完全可以忽略不计,用前三个主成分就可以很好地概括这组数据。
图6-19给出相关系数矩阵前4大特征值对应的特征向量,由此可以写出前三个主成分的表达式: 图6-19 原始变量对于各个主成分的因子载荷量 图6-19给出相关系数矩阵前4大特征值对应的特征向量,由此可以写出前三个主成分的表达式: Prin1 = 0.46x1* + 0.31x2* + 0.47x3* + 0.24x4* + 0.25x5* – 0.26x6* – 0.32x7* + 0.42x8* Prin2 = 0.26x1* – 0.40x2* + 0.11x3* – 0.49x4* + 0.50x5* + 0.17x6* + 0.40x7* + 0.29x8* Prin3 = 0.11x1* + 0.25x2* + 0.19x3* + 0.33x4* – 0.25x5* + 0.72x6* + 0.40x7* + 0.19x8*
可见,第一主成分中x3、x1、x8的系数最大;第二主成分中x5、x7具有较大的正系数,x4、x2则具有较大的负系数;第三主成分中x6的系数最大,远远超过其他指标的影响。因此,可以把第一主成分看成是由固定资产投资(x3)、GDP(x1)、工业总产值(x8)所刻画的反映经济发展水平的综合指标;把第二主成分看成是由货物周转量(x5)、职工平均工资(x4)、居民消费水平(x2)、商品零售价格指数(x7)所刻画的与人民生活水平有关的综合指标;把第三主成分单独看成是居民消费价格指数(x6)的影响指标。 最后输出的是数据集w1,其中包含前4个主成分Prin1~Prin4的得分。
按第一主成分和第二主成分的得分作图,又称为载荷图,代码如下: (4) 主成分的散点图 按第一主成分和第二主成分的得分作图,又称为载荷图,代码如下: proc plot data=w1 vpct=80; plot prin1*prin2 $ diqu='*'/ haxis=-3.5 to 3 by 0.5 HREF=-2,0,2 vaxis=-3 to 4.5 by 1.5 VREF=-2,0,2; run; 显示如图6-20。
广东、江苏、上海、山东的第一主成分取值较高,说明这些省市的经济发展水平较高,其次是浙江、辽宁、河北、河南、北京、天津等。 由于在第二主成分中职工平均工资与居民消费水平具有负的载荷量,因此处于右半图中的河北、河南、山东等地的职工平均工资与居民消费水平较低,商品零售价格指数较高;而左半图中上海、天津、海南、北京等地的职工平均工资与居民消费水平较高,商品零售价格指数较低。
6.2 因子分析 6.2.1 因子分析的概念与步骤 6.2.2 使用INSIGHT模块作因子分析 6.2 因子分析 6.2.1 因子分析的概念与步骤 6.2.2 使用INSIGHT模块作因子分析 6.2.3 使用FACTOR过程进行因子分析
6.2.1 因子分析的概念与步骤 1. 因子分析模型 设p维可观测的随机向量X = (X1,...,Xp)'(假定Xi为标准化变量,即E(Xi) = 0,Var(Xi) = 1,i = 1,2,…,p)表示为
或 X = AF + ε 上式称为因子模型,其中F1、F2、…、Fm称为公共因子,简称因子,是不可观测的变量;待估的系数阵A称为因子载荷阵,aij(i = 1,2,…,p;j = 1,2,…,m)称为第i个变量在第j个因子上的载荷(简称为因子载荷); ε称为特殊因子,是不能被前m个公共因子包含的部分。并且满足:cov(F,ε) = 0,即F,ε不相关; D(F) = Im,即F1、F2、…、Fm互不相关,方差为1;D(ε) = diag(12,22,…,p2),即ε1、ε2、…、εp互不相关,方差不一定相等,εi~N(0,i2)。 因子分析的目的就是通过模型X = AF + ε以F代替X,由于m < p,从而达到降维的愿望。
2. 因子分析模型中的几个统计特征 (1) 因子载荷aij的统计意义 由Xi = ai1F1 +…+ aimFm + εi,两边同乘以Fj,再求数学期望: E(XiFj)=ai1E(F1Fj)+…+aijE(FjFj)+…+aimE(FmFj)+E(εiFj) 从而有 rij = E(XiFj) = aij 即载荷矩阵中第i行,第j列的元素aij是第i个变量与第j个公共因子的相关系数,反映了第i个变量与第j个公共因子的相关程度。|aij| 1,绝对值越大,相关程度越高。在这种意义上公共因子解释了观测变量间的相关性。
(2) 变量共同度的统计意义 因子载荷矩阵第i行的元素平方和: 称为变量Xi的共同度(i = 1,2,…,p)。 对Xi = ai1F1 +…+ aimFm + εi两边求方差: 显然,若因子方差hi2大,剩余方差i2必小。而hi2大就表明Xi对公因子的共同依赖程度大。设Var(Xi) = 1,即所有的公共因子和特殊因子对变量Xi的贡献为1。如果hi2非常靠近1,则i2非常小,此时因子分析的效果好,从原变量空间到公共因子空间的转化性质好。可见hi2反映了变量Xi对公共因子F的依赖程度,故称hi2为变量Xi的共同度。
(3) 公共因子Fj方差贡献的统计意义 因子载荷矩阵A中各列元素的平方和: 称为公共因子Fj对X的贡献,是衡量Fj相对重要性的指标,qj2越大表明Fj对X的贡献越大。
3. 因子载荷矩阵的估计方法 给定p个相关变量X1,...,Xp的观测数据阵X,由X = AF + ε易推出 ∑ = AA' + D 其中∑ = D(X)为X的协方差阵,A = (aij)为p m的因子载荷阵,D = diag(12,22,…,p2)为p阶对角阵。 由p个相关变量的观测数据可得到协差阵的估计,记为S。为了建立因子模型,首先要估计因子载荷aij和特殊方差i2。常用的参数估计方法有以下三种:主成分法,主因子法和极大似然法。
(1) 主成分法 设样品协方差阵S的特征值为λ1≥λ2≥…≥λp≥0,u1,u2,…,up,为对应的标准化特征向量,当最后p–m个特征值较小时,S可近似地分解为: 其中 为pm阵, ,即得因子模型的一个解。载荷阵A中的第j列和X的第j个主成分的系数相差一个倍数(j = 1,…,m),故这个解称为主成分解。
(2) 主因子法 主因子方法是对主成分方法的修正,设R = AA' + D,则R* = R – D = AA'称为约相关矩阵,若已知特殊因子方差的初始估计 ,也就是已知变量共同度的估计: 则R*对角线上的元素是,而不是1。即:
计算R. 的特征值和特征向量,取前m个正特征值λ1. ≥λ2. ≥…≥λp. > 0,相应的特征向量为u1. ,u2. ,…,up R* = AA' 其中 ,令 (i = 1,…,p), 则A和D为因子模型的一个解,这个解称为主因子解。
在实际中特殊因子方差(或变量共同度)是未知的。以上得到的解是近似解。为了得到近似程度更好的解,常常采用迭代主因子法。即利用上面得到的 D* = diag( ) 作为特殊因子方差的初始估计,重复上述步骤,直到解稳定为止。 变量共同度hi2常用的初始估计有以下几种方法: ● 取第i个变量与其他所有变量的多重相关系数的平方; ● 取第i个变量与其他变量相关系数绝对值的最大值; ● 取1,它等价于主成分解。
(3) 极大似然法 假定公共因子F和特殊因子ε服从正态分布,那么可得到因子载荷阵和特殊因子方差的极大似然估计,设p维观测向量X(1),...,X(n)为来自正态总体Np(μ,∑)的随机样品,则样品似然函数为μ,∑的函数L(μ,∑)。 设∑= AA' + D,取μ = ,则似然函数为A,D的函数:(A,D),求A,D使达最大。为保证得到唯一解,可附加计算上方便的唯一性条件:A'D-1A = 对角阵,用迭代方法可求得极大似然估计A和D。
4. 因子旋转(正交变换) 所谓因子旋转就是将因子载荷矩阵A右乘一个正交矩阵T后得到一个新的矩阵A*。它并不影响变量Xi的共同度hi2,却会改变因子的方差贡献qj2。因子旋转通过改变坐标轴,能够重新分配各个因子解释原始变量方差的比例,使因子更易于理解。
设p维可观测向量X满足因子模型:X = AF +ε。T为正交阵,则因子模型可写为 X = ATT'F +ε = A*F* +ε 其中A* = AT,F* = T'F。 易知,∑ = AA' + D = A*A*' + D(其中A* = AT)。这说明,若A,D是一个因子解,任给正交阵T,A* = AT,D也是因子解。在这个意义下,因子解是不惟一的。 由于因子载荷阵是不惟一的,所以可对因子载荷阵进行旋转。目的是使因子载荷阵的结构简化,使载荷矩阵每列或行的元素平方值向0和1两极分化,这样的因子便于解释和命名。
有三种主要的正交旋转法:四次方最大法、方差最大法和等量最大法。这些旋转方法的目标是一致的,只是策略不同。 如果两种旋转模型导出不同的解释,这两种解释不能认为是矛盾的。倒不如说是看待相同事物的两种不同方法,是在公因子空间中的两个不同点。只取决于惟一的一种你认为是正确旋转的任何结论都是不成立的。 在统计意义上所有旋转都是一样的,即不能说一些旋转比另一些旋转好。因此,在不同的旋转方法之间进行的选择必须根据非统计观点,通常选择最容易解释的旋转模型。
Fji = j1xi1 + j2xi2 +…+ jpxip (j = 1,2,…,k) 5. 因子得分 计算因子得分的途径是用原有变量来描述因子,第j个因子在第i个样品上的值可表示为: Fji = j1xi1 + j2xi2 +…+ jpxip (j = 1,2,…,k) 式中,xi1,xi2,…,xip分别是第1,2,…,p个原有变量在第i个样品上的取值,j1,j2,…,jp分别是第j个因子和第1,2,…,k个原有变量间的因子值系数。可见,它是原有变量线性组合的结果(与因子分析的数学模型正好相反),因子得分可看作各变量值的加权(j1,j2,…,jp)总和,权数的大小表示了变量对因子的重要程度。
于是有: Fj = j1X1+j2X2+…+jpXp (j = 1,2,…,k) 上式称为因子得分函数。由于因子个数k小于原有变量个数p,故式中方程的个数少于变量的个数。因此,对因子值系数通常采用最小二乘意义下的回归法进行估计。可将上式看作是因子变量Fj对p个原有变量的线性回归方程(其中常数项为0)。可以证明,式中回归系数的最小二乘估计满足:Bj = Aj'R-1,其中 Bj = (j1,j2,…,jp),Aj' = (a1j,a2j,…,apj)为第1,2,…,p个变量在第j个因子上的因子载荷,R-1为原有变量的相关系数矩阵的逆矩阵。 由上式计算出因子变量Fj的因子值系数,再利用因子得分函数可算出第j个因子在各个样品上的因子得分。
6.2.2 使用INSIGHT模块作因子分析 【例6-4】今有20个盐泉,盐泉的水化学特征系数值见表6-6。试对盐泉水的化学分析数据作因子分析。 表6-6 盐泉水化学特征系数的数据 利用因子分析法,可揭示观察数据中7个指标之间的相互关系,寻找潜在的影响因子,并用这些潜在因子对原指标之间的相关关系进行解释。假定表6-6的数据已经存入数据集mylib.yq中。 序号 x1(矿化度) x2(Br·103/Cl) x3(K·103/盐) x4(K·103/Cl) x5(Na/K) x6(Mg·102/Cl) x7(Na/Cl) 1 11.835 0.480 14.360 25.210 25.21 0.810 0.98 2 45.596 0.526 13.850 24.040 26.01 0.910 0.96 … 19 304.092 0.283 0.789 1.357 438.36 0.193 1.01 20 202.446 0.042 0.741 1.266 309.77 0.290 0.99
1. 使用INSIGHT模块做因子分析的步骤 在INSIGHT模块中打开数据集Mylib.yq。 (1) 求相关系数阵及其特征值 选择菜单“Analyze”“Multivariate(Y X)(多元分析)”,打开“Multivariate(Y X)”对话框。将变量x1~x7选为Y变量,如图所示。
单击“Output”按钮,选中“Principal Component Analysis(主成分分析)”复选框,如图所示。 单击下面的“Principal Component Options(主成分选项)”按钮,打开“Principal Component Options”对话框,确认“Correlations(Structure)(相关(结构))”复选框被选中(默认状态),单击“OK”按钮返回;
两次单击“OK”按钮,得到因子分析结果。输出的数字分析结果包括5个部分:简单统计量、相关系数矩阵、相关系数矩阵的特征值以及默认的两个因子载荷阵等。 其中相关系数阵及其特征值等如图6-22所示。 结果显示,前三个特征值的方差贡献率依次为:0.6063、0.1788、0.1315。
(2) 建立因子载荷阵 由于前三个特征值的累积贡献率已达91.66%,故取前三个特征值建立因子载荷阵。 选择菜单“Tables”“Principal Components”,在弹出的“Principal Component Analysis”对话框中选择“3”个因子,及“Correlations(Structure)”选项,单击“OK” ,得到因子载荷阵如图所示。 由于第1、2公因子的载荷中有一些数值在0.5附近的中等载荷,其意义含糊不清,故考虑作因子旋转。
(3) 因子旋转 重新回到INSIGHT的数据窗口,选择菜单“Analyze”“Multivariate(Y X)”,打开“Multivariate(Y X)”对话框,将变量x1~x7选为Y变量。 首先,单击“Method”按钮,在打开的对话框中单击“Rotation Options”按钮,打开“Rotation Options”对话框,选择旋转方式为“Quartimax(最大四分位法)”,并修改“Components”的值为3,如图6-25所示。
然后,单击“Output”按钮,在打开的对话框中单击“Principal Component Analysis”复选框下面的“Principal Component Options”按钮,打开“Principal Component Options”对话框。 选中“Component Rotation”复选框(图左),单击“Rotation Options”按钮,打开“Rotation Options”对话框,增加选中“Output Component Scores”复选框和“Communality Estimates”复选框如图右所示。
结果包括正交旋转矩阵(Orthogonal Rotation Matrix)、旋转后的因子载荷阵(Rotation Correlations (Structure))(图左),以及各变量的共同度(图右)。 在数据集窗口还可以看到旋转前后的因子得分。
2. 因子分析的结果分析 从旋转后的因子模型(即因子载荷阵)中可以看出,相对于旋转前的因子模型,第一个公因子在x1、x5上的载荷增加,而在x2、x6、x7三个指标上的载荷明显减少。公因子1的载荷有正有负,正载荷主要是x5和x1,它们是钠盐形成的显示;负载荷主要是x3和x4,它们表示了钾盐形成的必要物质来源。 第二个公因子在x6(Mg·102/Cl)、x7(Na/Cl)两个指标上的载荷明显增加,这说明第二公因子是钾盐形成的条件的显示。 第三个公因子中起主要作用的是x2(Br·103/Cl),它是钾盐或钾矿化的一个环境标志。
回到INSIGHT数据窗口,用鼠标单击左上角的三角箭头,在弹出的菜单中选择“Extract”,打开“Extract”对话框,按下“Ctrl”键,用鼠标选定ID、RT1、RT2和RT3,如图6-28左所示,单击“OK”按钮,得到只包含编号及旋转后因子得分的数据子集如图6-29右。
利用数据窗口的排序功能,依次按三种公因子排序结果如图6-30所示。
图6-31是根据样品的因子得分,取RT1和RT2两个因子轴作因子得分图。 可见20个盐泉除第3号和7号外可分为三类:第一类为第14~20号盐泉,它们以第一因子轴上得分高,F2上得分绝对值低为特征;第二类为第8~13号盐泉,它们以F1上得分绝对值小,F2上得分为较大的负值为特征;第三类为第1~6号盐泉,它们以F1上得分为较大负值为特征。这三类表示三种不同的盐泉。
6.2.3 使用FACTOR过程进行因子分析 1. FACTOR过程简介 PROC FACTOR DATA = <数据集> <选项>; VAR <原始变量>; [PRIORS <共性值列表>;] [PARTIAL <变量列表>;] [FREQ <变量>;] [WEIGHT <变量>;] [BY <变量列表>;] RUN;
(1) PROC FACTOR语句 PROC FACTOR语句标志FACTOR过程的开始,同时还可通过设置其他语句定义数据集、指定具体分析方法和过程等。可设置的选项及其功能见表6-7。 通常只需要VAR语句作为PROC FACTOR语句的附加选项,其余均可省略。 (2) VAR语句 VAR语句用来指定需要分析的数值变量。如果该句省略,那么在其他语句中未做特殊规定的所有数值变量都将被分析。
(3) PARTIAL语句 如果想将因子分析建立在偏相关阵或协差阵的基础上,可用PARTIAL语句,以便程序将PARTIAL语句列出的变量的效果从整体分析中划分出来。 (4) PRIOR语句 PRIOR语句为每一个变量指定一个从0.0到1.0之间的初始共性方差估计值。第一个数值对应于VAR语句中的第一个变量,第二个数值对应第二个变量,依次类推。给出的数值个数必须与变量个数相等。 可以用“PROC FACTOR”语句中的“PRIORS =”选项指定各种各样的共性方差估计方法。
2. PROC SCORE得分过程 FACTOR过程的输出结果包括特征值情况、因子载荷、公因子解释比例,等等。为了计算因子得分,一般在PROC FACTOR语句中加一个SCORE选项和“OUTSTAT = 输出数据集”选项,然后用如下的得分过程计算公因子得分。 PROC SCORE DATA = <原始数据集> SCORE = <FACTOR过程的输出数据集> OUT = <得分输出数据集>; VAR <用来计算得分的原始变量集合>; RUN;
3. 实例分析 【例6-5】2004年31个省市自治区经济发展基本情况的八项指标,原始数据如表6-8所示。 表6-8 31个省市自治区经济发展基本情况 假定上述数据存放在数据集Mylib.jjfz中,试对经济发展基本情况的八项指标作因子分析。 地区 GDP x1 工业生产总值x2 固定资产投资x3 居民消费水平x4 货物周转量x5 居民消费价格指数x6 商品零售价格指数x7 职工平均工资x8 北京 4283.31 1290.16 2528.21 1354.23 537.7 100.952 99.249 29674 天津 2931.88 1436.73 1245.66 806.11 11223 102.255 100.83 21754 河北 8768.79 4086.43 3218.76 2619.18 4029.2 104.25 103.17 12925 山西 3042.41 1568.47 1443.88 1147.27 1417.7 104.14 103.06 12943 … 新疆 2200.15 745 1147.15 663.52 727.8 102.712 100.75 14484
结果给出8个变量的简单统计量,相关阵(略),相关阵的特征值、累计贡献(如图所示)。 (1) 主成分解 PROC FACTOR DATA= Mylib.jjfz SIMPLE CORR; Var x1 – x8; TITLE '8个经济指标的分析'; TITLE2 '主成分解'; RUN; 结果给出8个变量的简单统计量,相关阵(略),相关阵的特征值、累计贡献(如图所示)。
前两个主成分解释了84.60%的方差,按照缺省的选择因子个数的准则MINEIGEN,取大于1的特征值,所以取两个因子。 因子模型(factor pattern,或称因子载荷阵)为最重要的结果之一,如图所示。 它们是用公因子表示原始变量的回归系数。第一公因子在所有8个变量上都有正的载荷,可见这个因子反应了经济发展规模的影响,但载荷有大有小。第二公因子在居民消费价格指数和零售商品价格指数上有大的正载荷,反映了价格指标的影响。
结果还给出了公因子解释能力的估计(图6-34): 图6-34 各变量的共同度 Variance Explained by Each Facor给出了公因子对原始变量的解释能力(方差贡献)的量度,Final Communality Estimates:Total是两个公因子对原始变量的解释能力的总和。 最后一行给出每个原始变量的共同度,由于变量x5被两个因子解释的信息不够多,于是考虑选取3个公因子。
在FACTOR语句中增加选项n = 3,重做分析: (2) 选择公因子 在FACTOR语句中增加选项n = 3,重做分析: PROC FACTOR DATA= Mylib.jjfz n = 3; Var x1 – x8; RUN; 主要结果如图,三个公因子对原始变量的解释能力总和为7.421649,较前有所增加。 而且每一个变量的共同度都比较大,所以可以认为三个公因子可以很好地解释原始变量中的信息。但是由因子载荷系数看出,三个因子的含义不易解释,于是考虑作因子旋转。
方差最大旋转的旋转阵和旋转后的因子载荷阵如图6-36所示。 (3) 因子旋转 为了得到更好的因子解释,在上面的PROC FACTOR语句中再加上一个ROTATE = PROMAX旋转选项,这样将在得到主成分解后再进行方差最大正交旋转(VARIMAX),并加了一个REORDER选项使输出时把原始变量受相同因子影响的放在一起: PROC FACTOR DATA= Mylib.jjfz n = 3 ROTATE = VARIMAX REORDER; Var x1 – x8; RUN; 方差最大旋转的旋转阵和旋转后的因子载荷阵如图6-36所示。
从结果看得到的因子比正交旋转前有了较大改进。第一因子在GDP、工业生产总值、固定资产投资、居民消费水平这些指标上的载荷更加明确,因此第一因子可以被称为总量因子;第二因子在商品零售价格指数、居民消费价格指数和职工平均工资指标上的载荷也更加明确,第二因子可以被称为消费因子;尤其是第三因子在货物周转量指标上的载荷得到确立。
(4) 因子得分 为了产生因子得分,需要在FACTOR过程中使用SCORE选项和OUTSTAT=选项输出得分系数数据集。比如,为了计算方差最大正交旋转的因子得分,并做因子得分散点图,可以用如下程序: PROC FACTOR DATA=Mylib.jjfz n=3 ROTATE=VARIMAX REORDER SCORE OUT=OUT; RUN; PROC PRINT DATA=Out; VAR diqu factor1 factor2 factor3; PROC PLOT DATA = Out; PLOT factor2*factor1 $diqu = '*'/href=0 vref=0;
用回归法得到的因子得分系数如图所示。由此可以写出三个因子得分函数: F1 = 0.2655x1 + 0.2579x2 + 0.2808x3 + 0.2570x4 – 0.0161x7 – 0.0096x6 + 0.0530x8 – 0.1580x5 F2 = –0.01660x1 –0.0227x2 –0.0724x3 +0.0279x4 + 0.4154x7 + 0.3944x6 – 0.3292x8 + 0.1409x5 F3 = – 0.0693x1 – 0.0445x2 – 0.1561x3 – 0.0664x4 + 0.1784x7 + 0.1217x6 – 0.0055x8 + 1.1378x5
将31个地区的观测值代入以上因子得分函数,即得各地区的因子得分,本例由FACTOR过程自动给出计算结果,如图6左所示。根据前两个因子得分所作的因子得分图如图右所示。
第4象限中的北京、上海、广东、江苏等地区经济发展总量高,消费价格低,工资高;而第2象限中的云南、贵州、宁夏、海南等地区的经济发展总量低,消费价格高,工资低;第1象限中的河南、湖北、湖南等地区的经济发展总量高,消费价格高,工资低;第3象限中的西藏、天津、甘肃、青海等地区的经济发展总量低,消费价格低,工资高。