第二节自变量重新组合与主成分回归 上节讲的岭估计是一种有偏估计,它牺牲了无偏性的性质,但使估计的均方误差大为缩小, 从而使估计稳定而接近真实值。它主要用于自变量高度相关(复共线)的场合,也可以有于剔除 变量。这一节要介绍的主成分估计是由 w.F. Massy于1965年提出来的。它在上述诸方面近似 于岭估计,但处理资料的方法不一样。 、主成分回归的概念 先说去掉变量的原则。在回归模型 Y=lBo+XB+E, E(a)=0, Var(a)=ol (2.2.1) 中若有某一个变量,比如说X1,它的变差很小,即 (X1-X)2 则我们有理由把它并入常数项,实际就是把它剔除了 再谈主成分的思想。它包含两步:一是将自变量重新线性组合,二是在新的组合变量中去 掉那些变差小的,而留下主成分。设模型(221)的设计矩阵资料已中心化,即IX=0。模型 相关矩阵XX的特征根为A1≥2≥…≥Am,它们是下面特征方程的m个解 LXX-dm=0 (22.3) 特征根λ;对应的特征向量p满足下面的线性方程 LXX-A ImIp=0 (2.24) 选取标准化的向量P,使P1P1=1。如果这m个特征根都是单根,对应的m个特征向量一定 是两两正交的。若特征根有重根,我们也可以使特征向量正交化。我们利用这m个特征向量(每 个特征向量是m维)对原来变量作线性组合: 2,=pnX1+p2X2+…+pmXm=和,j=1…m 这样算完成了主成分回归的第一步:自变量重新组合 Z,=1,…;m也是变量,按其变差大小排序,去掉变差太小的Z,留下变差大的前r个 Z,便得到Y关于Z,产=1,…r的回归方程 Y=la+z+8 (2.2.6) 其回归系数就是所谓主成分估计。这就完成了主成分回归的第二步 至于这第二步具体如何计算,与原有相关矩阵的特征根λ1…,λm关系如何,要看下面的 分析
11 第二节 自变量重新组合与主成分回归 上节讲的岭估计是一种有偏估计，它牺牲了无偏性的性质，但使估计的均方误差大为缩小， 从而使估计稳定而接近真实值。它主要用于自变量高度相关(复共线)的场合，也可以有于剔除 变量。这一节要介绍的主成分估计是由 W. F. Massy 于 1965 年提出来的。它在上述诸方面近似 于岭估计，但处理资料的方法不一样。 一、主成分回归的概念 先说去掉变量的原则。在回归模型 n Y X E Var I 2 0 = 1 + + , ( ) = 0, ( ) = （2.2.1） 中若有某一个变量，比如说 X1,它的变差很小，即 = − n i X i Xi 1 2 ( 1 ) 0 （2.2.2） 则我们有理由把它并入常数项，实际就是把它剔除了。 再谈主成分的思想。它包含两步：一是将自变量重新线性组合，二是在新的组合变量中去 掉那些变差小的，而留下主成分。设模型(2.2.1)的设计矩阵资料已中心化，即 1X = 0 。模型 相关矩阵 X′X 的特征根为λ1≥λ2≥…≥λm,它们是下面特征方程的 m 个解 | X X − I m |= 0 （2.2.3） 特征根λi 对应的特征向量 pi 满足下面的线性方程 | X X − i I m | pi = 0 （2.2.4） 选取标准化的向量 pi，使 1 ' pi pi = 。如果这 m 个特征根都是单根，对应的 m 个特征向量一定 是两两正交的。若特征根有重根，我们也可以使特征向量正交化。我们利用这 m 个特征向量(每 个特征向量是 m 维)对原来变量作线性组合： Z j = p j1X1 + p j2X2 ++ p j m X m = Xp j , j = 1, ,m （2.2.5） 这样算完成了主成分回归的第一步：自变量重新组合。 Zj, j=1，…, m 也是变量，按其变差大小排序，去掉变差太小的 Zi, 留下变差大的前 r 个 Zj, 便得到 Y 关于 Zj, j=1，…, r 的回归方程 = + + 1 1 0 n rr Y Z （2.2.6） 其回归系数就是所谓主成分估计。这就完成了主成分回归的第二步。 至于这第二步具体如何计算，与原有相关矩阵的特征根λ1,…,λm 关系如何，要看下面的 分析
、主成分的确定 由于X已中心化,故B==∑y。对于新的组合变量(225,乙共有m个为一组 n 每组共进行了n次观测:Z1,产=1,…m,P=1,…,n。因为 1Z=1(如)=(1Xp,=0 故知Z,=1,…m也已中心化。因而 2n=0 (2.2.8) 于是新变量Z户=1,…m的变差度量为 ∑(z-2)2=∑2=Z乙,=pX,j=1 我们的任务是要将上面m个变差排序,去掉较小者。为此先叙述并证明下面的引理: 引理221设A为m阶对称阵,其特征根为1≥M2≥…≥Am,对应m个线性无关的特 征向量为p,p,…pm不妨设这些特征向量已经正交标准化。则对任意m维向量b,有 (1) max bAb= PrAp,=1: (2)mx、bAb=p4p+1=4° 证明(1)因为p,…pm是R中的标准正交基,故对任意b∈Rm,b可表为p1;…pm的线性组 合。设P=(p…m)则P为正交阵,P'P=E,P'AP=A,而b=pC,且b′b=C’P PC=C’C,于是 max Ab= max CPAPC= max C Ccel =mx(C21+…+Cnn)≤mx(λ1(C2+…+C2)=1 等号可在令C1=1,C2=…=Cm=0时达到,即在b=p方向达到 (2)因为b′p=0,=1,…k故b′P=(b′p;…,b'p,b'pk+1…b'pm)=(0,…;0,b pk+1…b′pm)而b′P=CP′P=C′E=(C1…CkC+1…Cm),故Cl=…=C=0。于是 bAb=mX、 CPAPC=mx(C4A+…+Cnn) ≤maxk+1(Ck+1+…+Cm)=入+1
12 二、主成分的确定 由于 X 已中心化，故 = = = n i Yi n Y 1 0 1 ˆ 。对于新的组合变量(2.2.5)，Zi 共有 m 个为一组， 每组共进行了 n 次观测：Zij, j=1，…,m, i=1,…,n。因为 1Z j = 1(Xp j ) = (1X ) p j = 0 （2.2.7） 故知 Zj, j=1，…,m 也已中心化。因而 = = = n i j Zij n Z 1 0 1 （2.2.8） 于是新变量 Zj, j=1，…,m 的变差度量为 Z Z Z Z Z p X Xp j m n i j j j j n i i j j i j ( ) , 1, , 1 1 − 2 = 2 = = = = = （2.2.9） 我们的任务是要将上面 m 个变差排序，去掉较小者。为此先叙述并证明下面的引理： 引理 2.2.1 设 A 为 m 阶对称阵，其特征根为λ1≥λ2≥…≥λm, 对应 m 个线性无关的特 征向量为 p1, p2,…,pm,不妨设这些特征向量已经正交标准化。则对任意 m 维向量 b，有 （1） 1 1 1 1 max = = = b Ab p Ap b b ； （2） 1 1 1 0, 1, , 1 max + + + = = = = k k = k b p i k b b b Ab p Ap i 。 证明 (1) 因为 p1,…,pm 是 R m 中的标准正交基，故对任意 b∈R m , b 可表为 p1,…,pm 的线性组 合。设 P= (p1,…,pm), 则 P 为正交阵，P′P=E, P′AP=Λ, 而 1 1 = m m mm b pC ，且 b′b=C′P′ PC=C′C, 于是 1 2 2 1 1 1 2 1 2 1 1 1 1 1 1 max( ) max( ( )) max max max = + + + + = = = = = = = = m C C m m C C m b b C C C C C C C C Ab C P APC C C 等号可在令 C1=1, C2=…=Cm=0 时达到，即在 b=p1 方向达到。 （2） 因为 b′pi=0, i=1,…,k, 故 b′P=(b′p1,…,b′pk, b′pk+1…b′pm)= (0,…,0, b′ pk+1…b′pm),而 b′P=C′P′P=C′E= (C1…Ck Ck+1…Cm)，故 C1=…=Ck=0。于是 1 2 2 1 1 1 2 1 2 1 1 0 1 0, 1, 1 max ( ) max max max( ) 1 + + + = + + = = = = = = = + + = = = + + k k m k C C k k m m C C C C C C b p i k b b C C b Ab C P APC C C i k
等号可在令Ck+1=1,其余C均为0时达到 证毕 这个引理结论(1)的几何含意是,二次型b′Ab在球面b′b=1的最大值为A的最大特征 根1,且在b=p1方向达到。结论(2)的几何含意是,在球面b′b=1上,除去p1方向,在剩余 各方向的线性组合中,二次型b′Ab在b=p方向达到最大值λ2,再除去p,p方向,在剩余 各方向的线性组合中,二次型b′Ab在b=p3方向达到最大值λ3等等。如图2221所示。 图2221 现在回到正题,要给(22.9)表达的新的组合变量的变差排序。因为X′X的特征根为1 A2≥…≥m,PP1=1,故由引理 mxp(XP=41,当P,=P (22.10) 我们称Xp1=Z1=p1+…+pmYm为原自变量X1,…Ym的第一主成分,它的变差为λ1。类似地 可得第二主成分Z2,其变差为A2,…,最后一个主成分为Zm,变差为λm。要去掉变差小的 就是去掉最小或较小特征根所对应的组合变量。因为XX非负定,其特征根非负,较小的特 征根是接近于0的。不妨设去掉了m-r个组合变量,留下了r个。于是模型(2.26成为 r=lao +z(ra(r)+E (2.2.11) 这里Z()=(Z1,…,Z),a(r)=(a1…,a)。又记A(n)=diag(1,2,…A),则典则形式(2.2 的最小二乘解为 r)=anzar 而典则形式(226)的解为 (ar:0) 由于z=XP,故a=P′B(Za=XP),于是回到原模型(22.1),令P(n)为P的前r列,得 B(r)=Pa= P(r)a(r) 我们称之为B的主成分估计 实际计算过程可以归纳如下:
13 等号可在令 Ck+1=1，其余 Ci 均为 0 时达到。 证毕 这个引理结论(1)的几何含意是，二次型 b′Ab 在球面 b′b=1 的最大值为 A 的最大特征 根λ1，且在 b=p1 方向达到。结论(2)的几何含意是，在球面 b′b=1 上，除去 p1 方向，在剩余 各方向的线性组合中，二次型 b′Ab 在 b=p2 方向达到最大值λ2, 再除去 p1, p2 方向，在剩余 各方向的线性组合中，二次型 b′Ab 在 b=p3 方向达到最大值λ3,等等。如图 2.2.2.1 所示。 图 2.2.2.1 现在回到正题，要给(2.2.9)表达的新的组合变量的变差排序。因为 X′X 的特征根为λ1 ≥λ2≥…≥λm, p j p j = 1 ，故由引理， 1 1 max ( ) = = j j p p p X X p j j ，当 p j = p1 （2.2.10） 我们称 Xp1=Z1=p11X1+…+p1mXm 为原自变量 X1,…,Xm 的第一主成分，它的变差为λ1。类似地， 可得第二主成分 Z2，其变差为λ2，…，最后一个主成分为 Zm，变差为λm。要去掉变差小的， 就是去掉最小或较小特征根所对应的组合变量。因为 X′X 非负定，其特征根非负，较小的特 征根是接近于 0 的。不妨设去掉了 m-r 个组合变量，留下了 r 个。于是模型(2.2.6)成为 =1 + ( )( ) + 0 Y Z r r （2.2.11） 这里 = nr Z(r) (Z1，…,Zr), ( ) ( , , ) 1 1 r r r = 。又记Λ(r)=diag(λ1,λ2，…,λr),则典则形式(2.2.11) 的最小二乘解为 r r Z r Y ' ( ) 1 ( ) ( ) ˆ − = （2.2.12） 而典则形式(2.2.6)的解为 ˆ ( ˆ 0) ( ) = r 由于 Z=XP, 故α=P′β(Zα=Xβ), 于是回到原模型(2.2.1)，令 P(r)为 P 的前 r 列，得 ( ) ˆ ( ) ˆ( ) ˆ r = P = P r r 我们称之为β的主成分估计。 实际计算过程可以归纳如下： P1 P2 P3 λ1 λ2 λ3
原模型Y0=X0)B0+g; 中心化Y=F)-y0,X=X0)-X(0 中心化模型Y=XB+ε; 计算相关阵X′X 及其特征根A1≥2…≥km∑= 特征向量p1,p2…,pm 取主成分r个,使(11+…+M)X≥75%;Pv=(P,p2;…p);Z=XPo; 典则化模型Y=1a+Za()+E:a()=AZ1F,M(r)=dng(1…,λ) 原中心标准化模型主成分估计B(r)=Pa 需要指出的是,舍掉的那些近似为0的特征根以及相应的主成分,正好反映了原来自变 量的复共线关系。因为若λ≈0,则≈0,这就是p1X1+…+ Pink=0,是一个复共线关系。体 会一下岭回归与主成分回归,都在处理复共线关系。岭回归是沿相关阵主对角加一常数,从而 使最小特征根变大,主成分回归则干脆去掉这些小特征根。一补一泻,治的是同一疾病 算例222法国有关进口总额的经济分析 我们这里选用 Malinvand于1966年作出的研究法国有关进口总额经济分析的实例,这个 实例曾被本书的全书参考书目中文献[8]选用。选一个老实例的原因是计算比较复杂,我们 用自编的程序实现了主成分回归的计算,与原例资料吻合,从而验证我们程序的正确性。有了 正确的程序,代入别的资料计算不是再简单不过的事情吗? 该例资料有11组,分别为法国自1949年至1959年11个年份的。资料打印在下面程序执 行过程中。考虑的因变量Y资料列是法国进口总额。自变量有三个,法国国内生产总值X1存 储量Ⅺ2,总消费量X。现在要作模型拟合,先考虑主成分回归。 程序行印出相关阵的三个特征根经自动按大小排序为1=1.9992,2=0.9982,y=0.0026 选择保留特征根个数为2。于是程序计算出典则化模型 y=Za -+ 的估计:α1=0.6899,α≥=0.1920。它告诉我们,对于典则化模型资料,回归方程是 =0.6899Z1+0.1920Z2。这个结果在文献[8]里没有印出。接着程序回到中心标准化模型 F=XB+ε 计算出它的主成分估计:B1=0.4806,B2=0.2215B3=0.4825。这个资料是根据公式 B=Pa计算出来的,它与文献[8]所列的结果吻合。它告诉我们,对于中心标准化资 料,回归方程是Y=0.4806X1+0.2215Z2+0.4825X3 接下来,程序对典则化模型γ=zα十ε进行一般最小二乘回归统计分析计算。先打印出典 则化模型资料,再打印出计算分析结果。注意它的自变量只有两个。代表的主成分是
14 原模型 Y (0)=X (0)β(0)+ε； 中心化 (0) (0) (0) (0) Y = Y −Y , X = X − X 中心化模型 Y=Xβ+ε； 计算相关阵 X′X； 及其特征根 λ1≥λ2≥…≥λm, * 1 = = m i i 特征向量 p1, p2,…, pm； 取主成分 r 个，使(λ1+…+λr)/λ*≥75%；P(r)= (p1, p2,…,pr)；Z(r)=XP(r)； 典则化模型 = + + Y 1 0 Z(r) (r) ； ˆ( ) , ( ) ( , , ) ( ) 1 1 (r) r diag r r = Z Y r = − ； 原中心标准化模型主成分估计 ( ) ( ) ( ) ˆ ˆ Pr r r = 。 需要指出的是，舍掉的那些近似为 0 的特征根以及相应的主成分，正好反映了原来自变 量的复共线关系。因为若λj≈0，则 Xpj≈0,这就是 pj1X1+…+pjmXm=0, 是一个复共线关系。体 会一下岭回归与主成分回归，都在处理复共线关系。岭回归是沿相关阵主对角加一常数，从而 使最小特征根变大，主成分回归则干脆去掉这些小特征根。一补一泻，治的是同一疾病。 算例 2.2.2 法国有关进口总额的经济分析 我们这里选用 Malinvand 于 1966 年作出的研究法国有关进口总额经济分析的实例，这个 实例曾被本书的全书参考书目中文献［8］选用。选一个老实例的原因是计算比较复杂，我们 用自编的程序实现了主成分回归的计算，与原例资料吻合，从而验证我们程序的正确性。有了 正确的程序，代入别的资料计算不是再简单不过的事情吗? 该例资料有 11 组，分别为法国自 1949 年至 1959 年 11 个年份的。资料打印在下面程序执 行过程中。考虑的因变量 Y 资料列是法国进口总额。自变量有三个，法国国内生产总值 X1,存 储量 X2，总消费量 X3。现在要作模型拟合，先考虑主成分回归。 程序行印出相关阵的三个特征根:经自动按大小排序为λ1=1.9992,λ2=0.9982,λ3= 0.0026。 选择保留特征根个数为 2。于是程序计算出典则化模型 Y=Zα＋ε 的估计：α1=0.6899，α2=0.1920。它告诉我们，对于典则化模型资料，回归方程是 Y=0.6899Z1+0.1920Z2。这个结果在文献［8］里没有印出。接着程序回到中心标准化模型 Y=Xβ＋ε 计算出它的主成分估计： 0.4825 ˆ 0.2215, ˆ 0.4806, ˆ 1 = 2 = 3 = 。这个资料是根据公式 ( ) ( ) ˆ ˆ = Pr r 计算出来的，它与文献［8］所列的结果吻合。它告诉我们，对于中心标准化资 料，回归方程是 Y = 0.4806X1＋0.2215Z2＋0.4825X3。 接下来，程序对典则化模型 Y=Zα＋ε进行一般最小二乘回归统计分析计算。先打印出典 则化模型资料，再打印出计算分析结果。注意它的自变量只有两个。代表的主成分是