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

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

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

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

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

5、满意的精度为止。若最终对X 共提取了 m 个成分ti .tm偏最小二乘回归将通过实施yk对ti .tm的回归,然后再表达成yk关于原变量XiXm的回归方程，k=i,2,q。1.2 计算方法推导为了数学推导方便起见，首先将数据做标准化处理。X 经标准化处理后的数 据矩阵记为Eo=（Eoi，Ep）n邓，Yj经标准化处理后的数据矩阵记为F0=（F0i，F0q）n p。第一步 记ti是Eo的第一个成分，Wi是Eo的第一个轴，它是一个单位向量，既 |wi| = i。记Ui是Fo的第一个成分，Ui=FoCi。ci是Fo的第一个轴，并且|lCill=i。如果要t1U1能分别很好的代表 X 与丫中的数据变异信

6、息，根据主成分分析原理，应该有Var(ujr maxVar(ti) max另一方面，由于回归建模的需要，又要求t1对U1有很大的解释能力，有典型相关 分析的思路，1与 5 的相关度应达到最大值，既r (t1u1) t max因此，综合起来，在偏最小二乘回归中，我们要求t1与 5 的协方差达到最大，既Cov(t1uj=、；Var (tJVar (uJr(t1uJ tmax正规的数学表述应该是求解下列优化问题，既maxW1C；S2=-(CiCi-i)=o(i-5)由式(i-2)(i-5)，可以推出2 i=1 12 /. 2 wiEoFoC二:：EWZ|,FoC ?记斗=2打=2.辺=w iEoFo

7、*所以-正是优化问题的目标函数值.把式(i-2)和式(i-3)写成1EoFoc= ri(i-6)1FoEow =0iCi(i-7)将式(i-7)代入式(i-6),有.2EoFoFoEowi=$wi(i-8)同理，可得 .2FoEoEoFoC - iG(i-9)可见,W1是矩阵EoFoFoEo的特征向量，对应的特征值为宀2.宀是目标函数值，它要 求s.tW1W1 i1c1c1=1因此，将在值。如果采用拉格朗日算法，记S=W1E0FoC1_1(W1W1-1)_2对 s 分别求关于W1C11和2的偏导并令之为零，有(CiCi-1)-S小小 = E0F0C1-21W1=(1 -2):s7C=FoEo

8、W1-22C1=0(1-3)-sT-=- (W1W1-1)=0(1-4)取最大值，所以，Wi是对应于EoFOFOEO矩阵最大特征值的单位特征向量.而另 一方面,Ci是对应于矩阵FoEoEoFo最大特征值刊2的单位特征向量.求得轴Wi和Ci后,即可得到成分ti = E0W1Ui= FoCi然后,分别求Eo和Fo对ti,Ui的三个回归方程E。二tipiEi(i-i0)Fo二UiqiFi(i-ii)Fo二tiriFi(i-i2)式中，回归系数向量是Pl二EotiHtTii2qiFou12l|Ui|2Fot12lltlll2而Ei,Fi,Fi分别是三个回归方程的残差矩阵第二步 用残差矩阵Ei和Fi取代

9、Eo和Fo，然后，求第二个轴个成分t2,U2,有t2=EiW2U2=FiC2二2二:：12, u2二w2EFC2IIQW2是对应于 矩阵EiFiFiEi最大特征值匕 的特征值FiEiEiFi最大特征值的特征向量计算回归系数(i-i3)(i-i4)(i-i5)W2和C2以及第二C是对应于矩阵P2Eit2lKTFit2llt2ll2因此，有回归方程Ei弋卩2E2IF t2r2F2如此计算下去，如果 X 的秩是 A,则会有Eo=tiPi tAPAFo二tiri .tArA- FA(i-i6)(i-i7)由于，ti,，tA均可以表示成Eoi,，Eop的线性组合，因此，式(1-17)还可以还原成 y二

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

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

12、部分除去某个样本点 i 的 所有样本点集合（共含 n-1 个样本点）,用这部分样本点并使用 h 个成分拟合一个回 归方程;第二部分是把刚才被排除的样本点 i 代入前面拟合的回归方程，得到yj在 样本点 i 上的拟合值？町（书。对于每一个 i=1,2,n,重复上述测试，则可以定义yj的预测误差平方和为PRESS有nPRESShj八（yij-?hj（4）2（1-18）i =1-12hkPRES%SS；h -4)k(1-22)定义丫的预测误差平方和为PRESSh，有pPRESS 八 PRESShj(1-19)j4显然,如果回归方程的稳健性不好，误差就很大,它对样本点的变动就会十分敏感, 这种扰动误差

13、的作用，就会加大PRESSh的值。另外,再采用所有的样本点，拟合含 h 个成分的回归方程。这是，记第 i 个样本点的预测值为？hji，则可以记Yj的误差平方和为SShj，有nSShj八(Yij-Yhji)2(1-20)i二定义丫的误差平方和为SSh，有pSS,八 S0(1-21)j m一般说来,总是有PRESSh大于SSh,而SSh则总是小于SSh。下面比较SSh和PRESSh。SSh是用全部样本点拟合的具有h-1 个成分的方程的拟合误差PRESSh增加了一个成分th，但却含有样本点的扰动误差。如果 h 个成分的回归方 程的含扰动误差能在一定程度上小于(h-1)个成分回归方程的拟合误差，则认为

14、增 加一个成分th,会使预测结果明显提高。因此我们希望(PRESSh/SShG的比值能 越小越好。在 SIMCA-P 软件中，指定(PRESSh/SShJl0.952即PRESS 岂 0.95SSh 4时,增加成分th就是有益的；或者反过来说,当.PRESS 0.95 SSh时,就认为增加新的成分th,对减少方程的预测误差无明显 的改善作用.另有一种等价的定义称为交叉有效性。对每一个变量yk,定义对于全部因变量丫,成分th交叉有效性定义为qxPRESShk22hkPRESShQ2 =1 _- =1hSSg)用交叉有效性测量成分th对预测模型精度的边际贡献有如下两个尺度。(1)当Qh一(1 -0

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

16、d2);w=v(:,i);%function c,m,v=norm1(C)%对数据进行标准化处理n,s=size(C);for i=1:nfor j=1:s c(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 为回归系数。% w 最大特征值所对应的特征向量。% wh

17、 处理后的特征向量% f0 回归的标准化的方程系数% FF 原始变量的回归方程的系数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);E(:,1:px)=E0-t(:,1)*(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(:

18、,px*i+1:px*i+px);w(:,i+2)=maxdet(B(:,px*i+1:px*i+px);% maxdet 为求最大特征值的函数%norm1 为标准化函数%截取标准化的因变量%求最大特征向量%提取主成分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(:,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(:,

19、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,:)=(regress1(y(:,j),t);%求回归系数end for j=1:px fori=1:py % 生成标准化变量的方程的系数矩阵w1=wh(:,1:j);zr1=(zr(i,1:j);f0(i,:,j)=(w1*zr1);endn

20、ormxy,meanxy,covxy=norm1(C); %n ormxy 标准化后的数据矩阵%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(ccy.*ccxx.*f0(:,:,j)./ccx)-meanxy(1,px+1:px+py);FF(:,:,j)=fff,ff;生成原始变量方程的常数项和系数矩阵end%function

21、 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.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

22、,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,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:pxendRdX

23、tt(1,m)=sum(RdXt(1,1:m);for 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)*(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