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

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

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

3、能 十分强大。在一次偏最小二乘回归分析计算后， 不但可以得到多因变量对多自变 量的回归模型，而且可以在平面图上直接观察两组变量之间的相关关系，以及观察样本点间的相似性结构。这种高维数据多个层面的可视见性， 可以使数据系统 的分析内容更加丰富，同时又可以对所建立的回归模型给予许多更详细深入的实 际解释。偏最小二乘回归的建模策略原理方法1.1建模原理设有q个因变量，yq和p自变量x_,Xp。为了研究因变量和自变量的统计关 系,我们观测了 n个样本点，由此构成了自变 量与因变量的数据表 X= Xi,，Xp和.Y=，yq。偏最小二乘回归分别在X与丫中提取出成分h和u1 （也就是说，匕是Xi，., Xp

4、的线形组合,u1是y-，yq的线形组合）.在提取这 两个成分时，为了回归分析的需要，有下列两个要求：（1）ti和Ui应尽可能大地携带他们各自数据表中的变异信息；（2）ti与Ui的相关程度能够达到最大。这两个要求表明，ti和Ui应尽可能好的代表数据表 X和丫,同时自变量的成分ti对因变量的成分Ui又有最强的解释能力。在第一个成分ti和Ui被提取后，偏最小二乘回归分别实施X对ti的回归以及丫对ui的回归。如果回归方程已经达到满意的精度，贝U算法终止；否则 , 将利用X被ti解释后的残余信息以及丫被ti解释后的残余信息进行第二轮的成分提取。如此往复，直到能达到一个较满意的精度为止。若最终对X共提取了

5、 m个成分ti . tm偏最小二乘回归将通过实施yk对ti . tm 的回归,然后再表达成y k关于原变量X i . X m 的回归方程，k=i,2，,q 。1.2计算方法推导为了数学推导方便起见，首先将数据做标准化处理。X经标准化处理后的数 据矩阵记为Eo=（Eoi，Eop）观，匕经标准化处理后的数据矩阵记为 F0 =（ F 0i ，F oq ） “ p。第一步 记ti是Eo的第一个成分，Wi是Eo的第一个轴，它是一个单位向量， 既 |wi | = io记Ui是Fo的第一个成分，Ui = Fo oci是Fo的第一个轴，并且llcFi。3如果要ti Ui能分别很好的代表X与丫中的数据变异信息，

6、根据主成分分析原理，应该有Var(uJ； maxVagi )max另一方面，由于回归建模的需要，又要求ti对ui有很大的解释能力，有典型相关分析的思路，ti与U1的相关度应达到最大值，既r(ti，Ui) 'max因此，综合起来，在偏最小二乘回归中，我们要求 t1与u1的协方差达到最大，既COV(tiUi)=、：Var(t)Var i 丿 r(tiUi) T max正规的数学表述应该是求解下列优化问题，既maxwe=ii _iC iC i 一 丨因此，将在|wi=i和Ilcil|2=i的约束条件下，去求(Wi E0 FoCi)的最大值。如果采用拉格朗日算法，记s=Wi E 0 F0 Ci

7、 ' i (Wi wi i) ' 2 (Ci Ci i)对s分别求关于wi Ci i和2的偏导并令之为零，有=E o F0 Ci 2九 i wi=0:w(1-2)-S = Fo E0 wi 2 ' 2 Ci=O.'Ci(1-3)-s=(Wi Wi i)=0(1-4)5js然后,分别求Eo和Fo对t1t= E 0W1U1= F0C1U1的三个回归方程E 0 = t1 p 1 E1(1-10)F 0 = U1 q 1 ' F 1(1-11)F 0= t1r 1F1(1-12)=(C1C1 T)=°(1-5)- 2由式(1-2)(1-5)，可以推出2

8、 爲=2 J = w 1 E o F° Ci =: E o wi, Foci -记二=21 =2 J =w1 EoFoS所以正是优化问题的目标函数值.把式(1-2)和式(1-3)写成E o F0C1 二 T1W1(1-6)F o EoW - 二心(1-7)将式(1-7)代入式(1-6),有' ' .2E o Fo F o EoW - ,TI1 w1(1-8)同理，可得' ' .2F o E°E oFoC - 1 C1(1-9)可见，W1是矩阵E'oF°F 'oE°的特征向量，对应的特征值为二2.已是目标函数

9、值，它要 求取最大值，所以,W1是对应于E oFoF 0 E 0矩阵最大特征值的单位特征向量.而另 一方面,C1是对应于矩阵F'oEoE'oFo最大特征值R2的单位特征向量.求得轴W1和C1后,即可得到成分6#式中，回归系数向量是#qiriE otilltl |2F ou12II Ui IIF o t1II ti |2(i-i3)(i-i4)(i-i5)而Ei,F i ,Fi分别是三个回归方程的残差矩阵第二用残差矩阵Ei和Fi取代Eo和Fo然后,求第二个轴W2和C2以及第 二个成分t2 ,U2 ,有t2 = Ei W2 u 2 = Fi c2d2 =： 12, u 2=w2Ei

10、Fic2W是对应于 矩阵E'iFiF'iEi最大特征值-22的特征值,C2是对应于矩 阵F 'i EiE'iFi最大特征值的特征向量.计算回归系数P2E it 2II t2 II27#r2F it22II t2 II因此，有回归方程Ei - t2 p 2 E 2Ft = t2r 2 F2如此，如果X的秩是A，则会有Eo 二 ti p i 'tA p aFo *r'i tAr - F由于,ti, ,tA均可以表示成Eoi,，Eop的线性组合，因此，式(i-i7)还可以还原*成yk =F0k关于Xj* =E°k的回归方程形式，即y=“Xi

11、川川七kpX*P FAkk=1,2,qFAk是残差距阵Fa的第k列。1.3交叉有效性下面要讨论的问题是在现有的数据表下，如何确定更好的回归方程。在许多 情形下，偏最小二乘回归方程并不需要选用全部的成分ti,，tA进行回归建模，而是可以象在主成分分析一样，采用截尾的方式选择前m个成分 （m ： A,A =秩（X）,仅用这m个后续的成分 切,tm就可以得到一个预测性较好 的模型。事实上，如果后续的成分已经不能为解释Fo提供更有意义的信息时，采用 过多的成分只会破坏对统计趋势的认识，引导错误的预测结论。在多元回归分析 一章中，我们曾在调整复测定系数的内容中讨论过这一观点。下面的问题是怎样来确定所应提

12、取的成分个数。在多元回归分析中，曾介绍过用抽样测试法来确定回归模型是否适于预测应用。我们把手中的数据分成两部分：第一部分用于建立回归方程，求出回归系数估计量 bB,拟合值?B以及残差均方和?B；再用第二部分数据作为实验点，代入刚才所求 得的回归方程，由此求出?T和哮0 一般地若有塔' 吭，则回归方程会有更好的预 测效果。若 汐；冷应B,则回归方程不宜用于预测。在偏最小二乘回归建模中，究竟应该选取多少个成分为宜，这可通过考察增加 一个新的成分后，能否对模型的预测功能有明显的改进来考虑。采用类似于抽样 测试法的工作方式，把所有n个样本点分成两部分：第一部分除去某个样本点i的 所有样本点集合

13、（共含n-1个样本点）,用这部分样本点并使用h个成分拟合一个回 归方程;第二部分是把刚才被排除的样本点i代入前面拟合的回归方程，得到y j在 样本点i上的拟合值？叭丄）。对于每一个i=1,2，,n,重复上述测试，则可以定义yj的预测误差平方和为PRESS N，有PRESShj2-Whj ( _4)(1-18)定义丫的预测误差平方和为PRESS h，有pPRESS h =' PRESS hj(1-19)j 4显然,如果回归方程的稳健性不好，误差就很大，它对样本点的变动就会十分敏感， 这种扰动误差的作用，就会加大PRESS h的值。另外,再采用所有的样本点，拟合含h个成分的回归方程。这是，

14、记第i个样本点的预测值为y?hji ,则可以记yj的误差平方和为SSh ,有n2SShj= V Wij - ?hji )(1-20)i -4定义丫的误差平方和为SSh ,有pSSh = 6 SShj(1-21)j 土一般说来,总是有PRESS h大于SSh ,而SSh则总是小于SSh。下面比较SSh和PRESS h。SSh丄是用全部样本点拟合的具有h-1个成分的方程的拟合误差；PRESS h增加了一个成分th,但却含有样本点的扰动误差。如果h个成分的回归方 程的含扰动误差能在一定程度上小于(h-1)个成分回归方程的拟合误差，则认为增 加一个成分.会使预测结果明显提高。因此我们希望 (PRESS

15、 h / SSh丄)的比值能 越小越好。在SIMCA-P软件中，指定(PRESS h / SSh)_ 0.95 2即.PRESS h乞0.95、. SSh丄时,增加成分th就是有益的；或者反过来说，当.PRESS： 095 .SSh时，就认为增加新的成分th,对减少方程的预测误差无明显 的改善作用.另有一种等价的定义称为交叉有效性。对每一个变量yk，定义2hk(1-22)PRESS hk=1 SS(h 4)k10对于全部因变量丫,成分th交叉有效性定义为11#(1-23)q' PRESS hkkPRESS h=1 -SS( h 2)#(h l)k用交叉有效性测量成分t h对预测模型精度

16、的边际贡献有如下两个尺度(1) 当Q；亠(1 0.95 2) = 0.0975时，th成分的边际贡献是显著的。显而易 见，Q： _ 0.0975与(PRESS h / SSh ) :：: 0.95 2是完全等价的决策原则。 对于k=1,2,q,至少有一个k,使得2Q h - 0.0975这时增加成分th,至少使一个因变量yk的预测模型得到显著的改善，因此也 可以考虑增加成分th是明显有益的。明确了偏最小二乘回归方法的基本原理、方法及算法步骤后，我们将做实证分析。附 录fun cti on w=maxdet(A)%求矩阵的最大特征值v,d=eig(A);n ,p=size(d);d1=d*o n

17、es(p,1);d2=max(d1);i=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)% px自变量的输入个数% py输入因变量的个数。% C输入的自变量和因变量组成的矩阵% t提取的主成分% q为回归

18、系数。% w最大特征值所对应的特征向量。% 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)*(E0'*t(:,1)/(t(:,1)'*t(:,1)'% 获得回归系数 p(:,1:px

19、)=(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)Jw(:,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(:,i+2)'*t(

20、:,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:pxq(:,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,:)=(regr

21、ess1(y(:,j),t)' %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 标准化后的数据矩阵%mea nx每一列的均值%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);

22、 ff=ccy.*f0(:,:,j)./ccx;fff=-(sum(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)&

23、#39;Rdyt=r.A2;RdYt=mean(Rdyt)for m=1:pxRdYtt(1,m)=sum(RdYt(1,1:m)');endfor j=1:pyfor m=1:pyRdytt(j,m)=sum(Rdyt(j,1:m)');endendfor j=1:pxfor m=1:pxRd(j,m)=RdYt(1,1:m)*(w(j,1:m)A2)');endendfor j=1:pxVIP(j,:)=sqrt(px*ones(1,px)./RdYtt).*Rd(j,:);end%function r,Rdxt,RdXt,RdXtt,Rdxtt=fun8x(px,

24、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.A2;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)');endend% for j=1:px % for m=1:px%Rd(j,m)=RdXt(1,1:m

25、)*(w(j,1:m)A2)');% end% end% for j=1:px% VIP(j,:)=sqrt(px*ones(1,px)./RdYtt).*Rd(j,:);% end%function t,u=TU(px,py,C)%t 提取的自变量的主成分% u 提取的因变量的主成分c=norm1(C);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);B=F0'*E0*E0'*F0;cc(:,1)=maxdet(B); u(:,1)=F0*cc(:,1);% function drew(px,py,c)X=c(:,1:px);Y=c(:,px+1:px+py);line,l=size(Y);t,q,w,wh,f0,FF=fun717(px,py,c);YY=X*FF(:,2:px+1,3)'+ones(line,1)*FF(:,1,3)'subplot(1,1,1,1)bar(f0(:,:,3)title('直方图 ')legend('SG','TZBFB','FHL','J