真正的好东西偏最小二乘回归多元线性回归分析+典型相关分析+主成分分析.docx

1、真正的好东西偏最小二乘回归多元线性回归分析+典型相关分析+主成分分析偏最小二乘回归是一种新型的 多元统计数据分析方法 ，它与 1983 年由伍德和阿巴诺等人首次提出。 近十年来， 它在理论、 方法和应用方面都得到了迅速的发展。密西根大学的弗耐尔教授称偏最小二乘回归为第二代回归分析方法。偏最小二乘回归方法在统计应用中的重要性主要的有以下几个方面：（1）偏最小二乘回归是一种 多因变量对多自变量 的回归建模方法。（2）偏最小二乘回归可以较好地解决许多以往用普通多元回归无法解决的问题。在普通多元线形回归的应用中， 我们常受到许多限制。 最典型的问题就是自变量之间的多重相关性 。如果采用普通的最小二乘方

2、法， 这种变量多重相关性就会严重危害参数估计， 扩大模型误差， 并破坏模型的稳定性。 变量多重相关问题十分复杂， 长期以来在理论和方法上都未给出满意的答案， 这一直困扰着从事实际系统分析的工作人员。 在偏最小二乘回归中开辟了一种有效的技术途径， 它利用对系统中的数据信息进行分解和筛选的方式， 提取对因变量的解释性最强的综合变量， 辨识系统中的信息与噪声， 从而更好地克服变量多重相关性在系统建模中的不良作用。（3）偏最小二乘回归之所以被称为第二代回归方法，还由于它可以实现多种数据分析方法的 综合应用 。偏最小二乘回归 =多元线性回归分析 +典型相关分析 +主成分分析由于偏最小二乘回归在建模的同时

3、实现了数据结构的简化， 因此，可以在二维平面图上对多维数据的特性进行观察， 这使得偏最小二乘回归分析的图形功能十分强大。 在一次偏最小二乘回归分析计算后， 不但可以得到多因变量对多自变量的回归模型， 而且可以在平面图上直接观察两组变量之间的相关关系， 以及观察样本点间的相似性结构。 这种高维数据多个层面的可视见性， 可以使数据系统的分析内容更加丰富， 同时又可以对所建立的回归模型给予许多更详细深入的实际解释。一、 偏最小二乘回归的建模策略 原理 方法11.1 建模原理设有 q 个因变量 y1,., 和 p 自变量 x1 ,., xp 。为了研究因变量和自变量yq的统计关系,我们观测了 n 个样

4、本点,由此构成了自变量与因变量的数据表X= x1 ,., xp 和.Y= y1,., yq 。偏最小二乘回归分别在 X 与 Y 中提取出成分 t1 和u (也就是说, t1 是 x1 ,., xp 的线形组合, u1是 y1 ,., yq 的线形组合).在提取这1两个成分时 ,为了回归分析的需要 ,有下列两个要求 :(1) t1 和 u1 应尽可能大地携带他们各自数据表中的变异信息 ;(2) t1 与 u1 的相关程度能够达到最大。这两个要求表明，t 和 u1 应尽可能好的代表数据表 X 和 Y,同时自变量的成分1t 对因变量的成分 u1 又有最强的解释能力。1在第一个成分t 和 u1被提取后

5、，偏最小二乘回归分别实施X 对 t1 的回归1以及 Y 对u 的回归。如果回归方程已经达到满意的精度，则算法终止；否则 ,1将利用 X 被 t1 解释后的残余信息以及 Y 被 t1 解释后的残余信息进行第二轮的成分提取。如此往复，直到能达到一个较满意的精度为止。若最终对 X 共提取了 m 个成分t1， ， t m ， 偏最小二乘回归将通过实施y 对 t1k， ， t m ， 的回归 ,然后再表达成y 关于原变量 xkm ， 的回归方程 ,k=1,2, ,q 。1 ， ， x1.2 计算方法推导为了数学推导方便起见,首先将数据做标准化处理。 X 经标准化处理后的数据矩阵记为E =( E001，

6、， E0 )pn ， Yj 经标准化处理后的数据矩阵记为pF =( F01 ， ， F0q )0n 。p第一步记t 是1E 的第一个成分， w1 是0E 的第一个轴，它是一个单位向量，0既|w1 |=1。记 1 。 cu 是 F0 的第一个成分， u1= F0 c1 是1F 的第一个轴，并且|c1 |=1。02如果要t1，u1能分别很好的代表 X 与 Y 中的数据变异信息，根据主成分分析原理，应该有Var( u1) maxVar(t1 ) max另一方面，由于回归建模的需要，又要求 t1 对u 有很大的解释能力，有典型相关1分析的思路， t1与u 的相关度应达到最大值，既1r（t1 ，u ）

7、max1因此，综合起来，在偏最小二乘回归中，我们要求t 与u1的协方差达到最大，既1Cov( t1，u )= Var ( ) ( ) r( t1t1 Var u11，u ) max1正规的数学表述应该是求解下列优化问题，既max E w F c,w1,c0 1 011s.tw c1 1 wc11112 =1 和|c 2 =1 的约束条件下，去求 ( w因此，将在 |w | 1 |1 E10 F0 c1 )的最大值。如果采用拉格朗日算法，记s=w1 E1 1 ( w1 w11)2 (c1 c1 1)对 s 分别求关于w1，c1 ， 1和2 的偏导并令之为零，有sw1= E2 w1 =0 (1 -

8、2)1 1sc1= F0 E0 w12 c1 =0 (1-3) 2s1=( w1 w1 1)=0 (1-4)3s2=( c 1 1)=0 (1-5)1 c由式(1-2)(1-5),可以推出 2 1 2 w E F c E w ,F c1 02 0 1 0 1 0 1记 1 2 2 w E F c ,所以, 1正是优化问题的目标函数值 .1 2 0 11 0把式(1-2)和式(1-3)写成E (1-6) F c w00 1 1 1F (1-7) E w c00 1 1 1将式(1-7)代入式(1-6),有2E (1-8) F F E w w0 00 0 1 1 1同理,可得 2F (1-9) E

9、 E F c c0 00 0 1 1 1 F F E可见,w1 是矩阵 E 0 0 的特征向量 ,对应的特征值为0021 . 1 是目标函数值 ,它要求取最大值 ,所以,w 是对应于 E 0 0 矩阵最大特征值的单位特征向量 .而另0 F F E1 0 E E F一方面, c1 是对应于矩阵 F 0 0 最大特征值0021 的单位特征向量 .求得轴w 和c1后,即可得到成分1t1 E w0 1u1 F c0 1然后,分别求E 和F0 对t1 ,u1的三个回归方程0E0 t p E (1-10)11 1F0 u q F (1-11)1 11F0 t r F (1-12)11 1式中,回归系数向量

10、是4E t01p (1-13)1 |t |2 1F u0 1q (1-14)1 |u |21 F t0 1r (1-15)1 |t |21而E1, F 1 ,F1分别是三个回归方程的残差矩阵 .第二步 用残差矩阵E 和F1 取代1E 和 F0 ,然后,求第二个轴 w2 和c2 以及第0二个成分t ,u2 ,有2t = E1 w22u = F1 c22 2 t ,u w E F c2 2 1 22 1w 是对应于矩阵 E 1 1 最大特征值 F F E12 122 的特征值 , c2 是对应于矩阵F 最大特征值的特征向量 .计算回归系数 E E F 1 11 1E t12p2 |t |2 2 F

11、 t1 2r2 |t |22因此,有回归方程E1 t p E22 2F1 t r F22 2如此计算下去 ,如果 X 的秩是 A,则会有 E0 t （1-16）ptA p1 A1 F0 t r tAr F （1-17）1 A1 A由于,t1, ,tA 均可以表示成 E01, ,E0 p 的线性组合 ,因此,式(1-17)还可以还原5成* 关于yk F0kxj * E0 的回归方程形式，即k* *yk * x x F k=1,2, ,qpk1 1 kp AkF 是残差距阵FA 的第 k 列。Ak1.3 交叉有效性下面要讨论的问题是在现有的数据表下 ,如何确定更好的回归方程。在许多情形下 ,偏最小

12、二乘回归方程并不需要选用全部的成分 t1, ,tA 进行回归建模 ,而是 可 以 象 在 主 成 分 分 析 一 样 , 采 用 截 尾 的 方 式 选择前 m 个 成 分(m A, A 秩( X) ,仅用这 m 个后续的成分 t1, ,tm 就可以得到一个预测性较好的模型。 事实上 ,如果后续的成分已经不能为解释F 提供更有意义的信息时,采用0过多的成分只会破坏对统计趋势的认识,引导错误的预测结论。在多元回归分析一章中 ,我们曾在调整复测定系数的内容中讨论过这一观点。下面的问题是怎样来确定所应提取的成分个数。在多元回归分析中 ,曾介绍过用抽样测试法来确定回归模型是否适于预测应用。我们把手中的

13、数据分成两部分 :第一部分用于建立回归方程 ,求出回归系数估计量b ,拟合值y?B 以及残差均方和B2?B ;再用第二部分数据作为实验点 ,代入刚才所求得的回归方程 ,由此求出2y?T 和 ?T 。一般地 ,若有2?T2?B ,则回归方程会有更好的预测效果。若2?T2?B ,则回归方程不宜用于预测。在偏最小二乘回归建模中 ,究竟应该选取多少个成分为宜 ,这可通过考察增加一个新的成分后 ,能否对模型的预测功能有明显的改进来考虑。采用类似于抽样测试法的工作方式 ,把所有 n 个样本点分成两部分 :第一部分除去某个样本点 i 的所有样本点集合 (共含 n-1 个样本点 ),用这部分样本点并使用 h个

14、成分拟合一个回归方程 ;第二部分是把刚才被排除的样本点 i 代入前面拟合的回归方程 ,得到 yj 在样本点 i 上的拟合值y? 。对于每一个 i =1,2, ,n,重复上述测试 ,则可以定义hj( i )yj的预测误差平方和为PRESS ,有hjn2PRESShj ( y y? ( ) (1-18)ij hj i )i 16定义 Y 的预测误差平方和为PRESS ,有hpPRESSh PRESShjj 1(1-19)显然,如果回归方程的稳健性不好 ,误差就很大 ,它对样本点的变动就会十分敏感 ,这种扰动误差的作用 ,就会加大PRESS 的值。h另外,再采用所有的样本点 ,拟合含 h 个成分的回

15、归方程。这是 ,记第i 个样本点的预测值为y? ,则可以记hjiy 的误差平方和为 SShj ,有jn 2SShj (y y? ) (1-20)ij hjii 1定义 Y 的误差平方和为SS ,有hpSSh SShjj 1(1-21)一般说来,总是有PRESS 大于 SSh ,而SSh 则总是小于 SSh 1。下面比较 SSh 1和hPRESS 。 SSh 1 是用全部样本点拟合的具有 h-1 个成分的方程的拟合误差 ;hPRESS 增加了一个成分 th ,但却含有样本点的扰动误差。 如果 h 个成分的回归方h程的含扰动误差能在一定程度上小于 (h-1)个成分回归方程的拟合误差 ,则认为增加一

16、个成分t ,会使预测结果明显提高。因此我们希望 ( PRESSh / SSh 1 ) 的比值能h越小越好。在 SIMCA-P 软件中,指定(PRESSh / SSh21) 0.95即 PRESSh 0.95 SSh 1 时 , 增加 成分 th 就是 有益 的 ; 或 者反 过 来说 , 当PRESSh 0.95 SSh 时,就认为增加新的成分 th ,对减少方程的预测误差无明显1的改善作用 .另有一种等价的定义称为交叉有效性。对每一个变量y ,定义kPRESS2 1 hk (1-22) Qhk SS( h 1)k7对于全部因变量 Y,成分t 交叉有效性定义为hqPRESShkPRESS2 1

17、 k 1 1 hQ (1-23)h SS( h 1)SS(h 1) k用交叉有效性测量成分t 对预测模型精度的边际贡献有如下两个尺度。h2 2(1) 当 Q (1 0.95 ) 0.0975时, th 成分的边际贡献是显著的。显而易h2见, Q 0. 0975与h2(PRESSh / SSh 1) 0.95 是完全等价的决策原则。(2) 对于 k=1,2, ,q,至少有一个 k,使得2Qh1.4这时增加成分t ,至少使一个因变量 yk 的预测模型得到显著的改善,因此 ,也h可以考虑增加成分t 是明显有益的。h明确了偏最小二乘回归方法的基本原理、方法及算法步骤后，我们将做实证分析。附 录func

18、tion w=maxdet(A)%求矩阵的最大特征值v,d=eig(A);n,p=size(d);d1=d*ones(p,1);d2=max(d1);8i=find(d1=d2);w=v(:,i);%function c,m,v=norm1(C)%对数据进行标准化处理n,s=size(C);for i=1:nfor j=1:sc(i,j)=(C(i,j)-mean(C(:,j)/sqrt(cov(C(:,j);endendm=mean(C);for j=1:sv(1,j)=sqrt(cov(C(:,j);end%function t,q,w,wh,f0,FF=fun717(px,py,C)%

19、px 自变量的输入个数% py 输入因变量的个数。% C 输入的自变量和因变量组成的矩阵% t 提取的主成分% q 为回归系数。% w 最大特征值所对应的特征向量。9% wh 处理后的特征向量% f0 回归的标准化的方程系数% FF 原始变量的回归方程的系数c=norm1(C); %norm1 为标准化函数y=c(:,px+1:px+py); % 截取标准化的因变量E0=c(:,1:px);F0=c(:,px+1:px+py);A=E0*F0*F0*E0;w(:,1)=maxdet(A); % 求最大特征向量t(:,1)=E0*w(:,1); % 提取主成分E(:,1:px)=E0-t(:,1

20、)*(E0*t(:,1)/(t(:,1)*t(:,1);% 获得回归系数p(:,1:px)=(E0*t(:,1)/(t(:,1)*t(:,1);for i=0:px-2B(:,px*i+1:px*i+px)=E(:,px*i+1:px*i+px)*F0*F0*E(:,px*i+1:px*i+px);w(:,i+2)=maxdet(B(:,px*i+1:px*i+px);% maxdet 为求最大特征值的函数t(:,i+2)=E(:,px*i+1:px*i+px)*w(:,i+2);p(:,px*i+px+1:px*i+2*px)=(E(:,px*i+1:px*i+px)*t(:,i+2)/(t

21、(:,i+2)*t(:,i+2);E(:,px*i+px+1:px*i+2*px)=E(:,px*i+1:px*i+px)-t(:,i+2)*(E(:,px*i+1:px*i+px)*t(:,i+2)/(t(:,i+2)*t(:,i+2);endfor s=1:px10q(:,s)=p(1,px*(s-1)+1:px*s);endn,d=size(q);for h=1:pxiw=eye(d);for j=1:h-1iw=iw*(eye(d)-w(:,j)*q(:,j);endwh(:,h)=iw*w(:,h);endfor j=1:pyzr(j,:)=(regress1(y(:,j),t);

22、% 求回归系数endfor j=1:pxfori=1:py %生成标准化变量的方程的系数矩阵w1=wh(:,1:j);zr1=(zr(i,1:j);f0(i,:,j)=(w1*zr1);endnormxy,meanxy,covxy=norm1(C); %normxy 标准化后的数据矩阵11%meanxy每一列的均值%covxy每一列的方差ccxx=ones(py,1)*meanxy(1,1:px);ccy=(covxy(1,px+1:px+py)*ones(1,px);ccx=ones(py,1)*(covxy(1,1:px);ff=ccy.*f0(:,:,j)./ccx;fff=-(sum(

23、ccy.*ccxx.*f0(:,:,j)./ccx)-meanxy(1,px+1:px+py);FF(:,:,j)=fff,ff; % 生成原始变量方程的常数项和系数矩阵end%function r,Rdyt,RdYt,RdYtt,Rdytt,VIP=fun8y(px,py,c)X=c(:,1:px);Y=c(:,px+1:px+py);x=norm1(X);y=norm1(Y);t,q,w=fun717(px,py,X,Y);r1=corrcoef(y,t);r=r1(py+1:px+py,1:py);Rdyt=r.2;RdYt=mean(Rdyt)for m=1:pxRdYtt(1,m)=

24、sum(RdYt(1,1:m);endfor j=1:pyfor m=1:pyRdytt(j,m)=sum(Rdyt(j,1:m);12endendfor j=1:pxfor m=1:pxRd(j,m)=RdYt(1,1:m)*(w(j,1:m).2);endendfor j=1:pxVIP(j,:)=sqrt(px*ones(1,px)./RdYtt).*Rd(j,:);end%function r,Rdxt,RdXt,RdXtt,Rdxtt=fun8x(px,py,c)X=c(:,1:px);Y=c(:,px+1:px+py);x=norm1(X);y=norm1(Y);t,q,w=fun717(px,py,X,Y);r1=corrcoef(x,t);r=r1(px+1:px+px,1:px);Rdxt=r.2;RdXt=mean(Rdxt);for m=1:pxRdXtt(1,m)=sum(RdXt(1,1:m);endfor j=1:pxfor m=1:pxRdxtt(j,m)=sum(Rdxt(j,1:m);endend13% for j=1:px% for m=1:px% Rd(j,m)=RdXt(1,