主成分回归分析的EXCEL实现

2015-05-30林建华

宁波职业技术学院学报 2015年5期

关键词：多元回归分析主成分分析

林建华

摘要： EXCEL是一款功能非常强大的办公软件，文章利用其内置的公式和函数给出了主成分回归分析的完整算法和详细过程，得到的计算结果与专业统计软件给出的相同。因此，应用EXCEL实现主成分回归分析是可行的。

关键词：主成分分析；多元回归分析； EXCEL

中图分类号： O 245 文献标志码： A 文章编号： 1671-2153（2015）05-0097-05

0 引言

主成分回归是一种利用主成分作回归的方法，它既可以消除变量间的多重共线性，又能够保留全部的解释变量，回归方程的拟合效果也能够得到提高。主成分回归分析运算量较大，一般需要用专业统计软件来处理。这在一定程度上限制了不熟悉专业软件的使用者对主成分回归分析的实际应用。EXCEL是目前一款非常流行的实用办公软件，它不但具有功能强大的电子表格，而且还包含大量的函数来供数据分析和数据运算，以解决研究或实际工作中的问题。目前，虽然有应用EXCEL做主成分分析和多元线性回归的一些公开文献[1-3]，但尚未发现有利用EXCEL做主成分回归分析的具体算例及完整的实现过程。鉴于EXCEL使用的普遍性和广泛性，本文试图利用文献[4]提供的算例对应用EXCEL实现主成分回归分析的主要步骤和过程作详细的阐述。

1 基本原理

1.1 确定主成分

设有n个样品，每个样本观测p个指标，原始数据矩阵为x=（x1，x2，…，xp），将原始数据标准化，可得

X=X11 X12 … X1pX21 X22 … X2p… … … …Xn1 Xn2 … Xnp=（X1，X2，…，Xp）。（1）

计算标准化数据矩阵X的相关矩阵R=（rij）p×p，解特征方程|R-λI|=0，可得到R的特征根λ1≥λ2≥…λp≥0，在此基础上分别求出对应于特征值λi（i=1，2，…，p）的单位特征向量αi=（α1i，α2i，…，αpi）。最后，求得p个主成分Fi=α1iX1+α2iX2+…+αpiXp（i=1，2，…，p）。

1.2 进行回归分析

用标准化的原变量观测数据计算前m（mY=θ0+θ1F1+…+θmFm，（2）

由于主成分Fi（i=1，2，…，m）是标准化变量X1，X2，…，Xp表示的方程，所以再把m个主成分Fi代入式（2）得到：

Y=β0+β1X1+…+βpXp，（3）

其中β0=θ0，βi=θ1αi1+θ2αi2+…+θmαim（i=1，2，…，p）。

对回归方程进行统计检验，包括拟合优度检验（R2检验）、方程总体线性显著性检验（F检验）和变量的显著性检验（t检验）。

将式（3）中的X1，X2，…，Xp和Y还原为用原始变量x1，x2，…，xp和y来表示，即得到用原始变量表示的回归方程：

y=b0+b1x1+…+bpxp。（4）

值得注意的是，只有在原变量间存在较强的多重共线性时，采用主成分回归才会收到比较好的效果，而不是在任何情况下做主成分回归都有满意的效果，每种方法都有它的适用范围和局限性。

2 EXCEL实现过程

在EXCEL工作表Sheet1的单元格区域A1：E12，输入文献[4]所提供的原始数据，如表1所示。为了确保计算结果的精度，本文约定除在图表和有关方程式中只显示数字的几位小数外，中间运算过程中的数字尽可能保留多位小数。

2.1 原始数据的标准化处理

按式（5）对表1提供的原始数据x1，x2，x3和y进行中心标准化处理，即

Xij=（xij-j） Yij=（yij-j）。（5）

打开EXCEL工作表Sheet1，在单元格F2中输入公式：“=（B2-AVERAGE（B2：B12））/SQRT（DEVSQ（B2：B12））”，按回车键确认，单元格F2中就会显示数值-0.4774。然后，用鼠标按住F2右下角垂直下拖至F12单元格，便可以生成自变量x1的中心标准化数值X1。同理，可得到自变量x2和x3及因变量y的中心标准化数值X2，X3和Y，如表2所示。

2.2 计算样本相关系数矩阵

单击“数据”菜单中的“数据分析”，在弹出的“数据分析”窗口中选择“相关系数”，点击“确定”后，打开对话框，在“输入区域（I）：”输入“F2：H12”，在“输出区域（O）：”输入“F15：H17”，点击“确定”按钮，就会在输出区域给出相关系数矩阵的下三角的数值。相关系数矩阵系对称矩阵，上三角的数值与下三角的数值相等，可以采用复制、选择性粘贴、转置在另外单元格区域生成上三角的数值，之后，再把上三角的数值复制粘贴到相关系数矩阵的上三角的空白部分，计算结果如表3所示。

为方便以后的运算与使用，需要对EXCEL中所建立的相关系数矩阵A和单位矩阵I进行命名。选定相关系数矩阵所在的单元格区域G16：I18，单击编辑栏左端“名称”栏，输入自定义名称“A”后按回车键确认，以后在当前工作簿的所有工作表中，名称“A”就是指单元格区域G16：I18。同理，在单元格区域G23：I25建立3阶单位矩阵，并将其命名为“I”。

2.3 计算相关系数矩阵A的特征值及其单位特征向量

求解特征方程|A-λI|=0，可得到A的3个非负的特征值λ1、λ2、λ3，且λ1≥λ2≥λ3≥0。特征值λi（i=1，2，3）是主成分的方差。由Gersgorin圆盘定理[5]可知，λ1、λ2、λ3落在区间[0，2]内。用EXCEL求解λi（i=1，2，3）分两步，先确定每个特征值所在的小区间，然后利用“规划求解”命令求出每个特征值的更精确的近似值。

操作步骤如下：在单元格A1输入0，A2输入0.001（视求解精度而定），用鼠标按住A1：A2右下角垂直下拖至A20001单元格，在单元格B1输入公式：“=MDETERM（A-A1×I）”，用鼠标按住B1右下角垂直下拖至B20001单元格。观察B1：B20001单元格区域数值的正负符号的变化情况，每一次符号改变都说明这两个返回值之间有一个特征值。据此可得3个特征值λi（i=1，2，3）所在的小区间为[0.002，0.003]，[0.9981，0.9982]和[1.9991，0.9992]。单击“数据”菜单中的“规划求解”命令，打开“规划求解参数”对话框，在“设置目标：（T）”中输入“B1”，点击“目标值：（V）输入“0”，在“通过更改可变单元格：（B）中，输入或选择单元格区域“A1：A3”，单击“添加（A）”按钮，在“遵守约束：（U）栏中，分别添加条件：A1>=0.002，A1<=0.003，A2>=0.9981，A2<=0.9982，A3>=1.9991，A3<=1.9992，B2=0，B3=0，单击“求解（S）”按钮完成操作。此时单元格A1，A2，A3中的值分别显示为：0.0026908908，0.9981541760，1.9991549345，它们就是相关系数矩阵A的3个特征值，即λ1=1.999155，λ2=0.998154和λ3=0.002691。

对应于λi（i=1，2，3）的特征向量可由下面的齐次方程组解出：

（r11-λi）ψ1+r12ψ2+r13ψ3=0r11ψ1+（r12-λi）ψ2+r13ψ3=0r11ψ1+r12ψ2+（r13-λi）ψ3=0，（6）

将相关矩阵A的值和λi值代入方程组，并令ψ3=1，把方程组左边的最后一项移到等式的右边作为常数项，用前面的2个方程联立求解得到齐次方程的一组解，对解向量进行标准化可得到特征值λi对应的单位特征向量。

下面以λ1=1.999155为例，采用逆矩阵方法求解特征向量，详细过程如下（见表4）。

（1）输入系数矩阵和常数项矩阵。将λ1代入方程组，在单元格区域A1：B2输入系数矩阵，在区域C1：C2输入常数矩阵。

（2）求系数矩阵的逆矩阵。在输入系数矩阵和常数矩阵后，选取单元格区域A5：B6，在活动单元格输入公式：“=MINVERSE（A1：B2）”，按下Ctrl+Shift+Enter键，区域A5：B6就会显示系数矩阵的逆矩阵区域。

（3）求齐次方程的一组解。选取单元格区域D5：D6，在活动单元格输入公式：“=MMULT（A5：B6，C1：C2）”，按下Ctrl+Shift+Enter键，单元格D5和D6分别显示数字“0.9997”和“0.0616”，并在单元格D7输入“1”，这样就得到齐次方程的一组解。

（4）对解向量进行标准化。利用单位化公式αi=ψi/对解向量进行标准化，得到单位化特征向量αi=（α1i，α2i，α3i）。在单元格E5输入：“=D5/SQRT（SUMSQ（D5：D7））”，按下Ctrl+Shift+Enter键，用鼠标按住E5右下角垂直下拖至E7单元格，可以得到λ1的单位特征向量为α1=（0.7063，0.0435，0.7065）。依此类推，将λ2和λ3代入方程组，改变A1：B2主对角线单元格的值，求得其他两个特征值所对应的单位特征向量。即λ2的单位特征向量为α2=（-0.0357，0.999，-0.0258）；λ3的单位特征向量为α3=（-0.7070，-0.0070，0.7072）。值得注意的是，EXCEL中特征向量的正负号选择是一个比较重要的问题，具体如何确定可参考文献[6]，本文恕不赘述。

2.4 计算方差贡献率和累计贡献率

第k（k=1，2，3）个主成分的方差贡献率为

λk λj，前m个主成分的累计贡献率为λj λj。一般情况下，若m个主成分的累计贡献率超过85%，就认为m个主成分基本包含了原来指标的信息，即可以取前m个主成分进行分析和评价。方差贡献率的计算比较简单，如第1主成分的方差贡献率为1.9991/3=66.64%，前2个主成分的方差累计贡献率（1.9991+0.9981）/3=99.91%，因前2个主成分的方差累计贡献率为99.91%（>85%），故保留前2个主成分，而把第3个主成分剔除。

2.5 确定主成分表达式及计算其得分

根据λ1和λ2对应的单位特征向量，可得到以下2个主成分：

F1=0.7063X1+0.0435X2+0.7065X3F2=-0.0357X1+0.999X2-0.0258X3。（7）

将λ1和λ2对应的单位特征向量，输入单元格区域J3：J5和J9：J11。在单元格K2输入公式：“=MMULT（F2：H12，J3：J5）”，按下Ctrl+Shift+Enter键，用鼠标按住K2右下角垂直下拖至K12单元格，可得到第一主成分F1的得分。同理，可在单元格区域L2：L12得到第二主成分F2的得分，计算主成分F1和F2得分，如表5所示。

2.6 建立主成分回归模型

单击“数据”菜单中的“数据分析”命令，选中“数据分析”这一选项框，在分析工具（A）中选择“回归”，单击“确定”按钮。弹出“回归”对话框，在“Y值输入区域（Y）：输入“I2：I12”

在“X值输入区域（X）：输入“K2：L12”，再对其他一些选项进行适当选择，单击“确定”按钮完成操作（见图1）。输出的线性回归结果如图2所示。

对回归方程进行统计检验。从图2的输出结果可知，R2=0.9883，R2=0.9854，两者都接近于1。这表明自变量和因变量的线性相关性较强。F分布和t分布的临界值不用查表，可以用EXCEL的函数生成。在显著性水平α=0.05下，Fα（2，8）=FINV（0.05，2，8）=4.46，F=337.23>Fα（2，8），t（8）=TINV（0.05，8）=2.306，t1=25.486，t2=4.993>t（8），回归方程通过R2检验、F检验和t检验。

从图2提供的回归结果可知，回归系数估计值为θ1=0.690，θ2=0.191，常数项θ0≈0，得到回归方程为

Y=0.69F1+0.191F2-3.783E-17，（8）

将F1和F2代入式（8），可得：

Y=0.4804X1+0.2211X2+0.4826X3，（9）

将式（5）代入式（9），可得：

=0.4804（）+

0.2211（）+0.4826（），（10）

对式（10）进行整理，可得以原始变量表示的回归方程：

y=-9.1301+0.07278x1+0.60922x2+0.10626x3。（11）

这一结果与SAS[7]和Visual Basic[8]编程得到的计算结果完全相同，与SPSS[9]处理结果也近似相等。这就验证了上述利用EXCEL做主成分回归的方法是正确的，得到的结果是可靠的。

3 结束语

主成分回归计算过程较复杂，运算量较大。EXCEL是一款最为常用的数据处理软件，功能非常强大，完全可以满足主成分回归分析的计算需要，前面的计算例子已经验证了这一点。因此，对于一些不熟悉专业统计软件的使用者来说，利用EXCEL做主成分回归具有极大的便利性。