偏最小二乘回歸=多元線性回歸分析+典型相關分析+主成分分析

偏最小二乘回歸是一種新型的多元統(tǒng)計數(shù)據(jù)分析方法，它與1983年由伍德和阿巴諾等人首次提出。近十年來，它在理論、方法和應用方面都得到了迅速的發(fā)展。密西根大學的弗耐爾教授稱偏最小二乘回歸為第二代回歸分析方法。偏最小二乘回歸方法在統(tǒng)計應用中的重要性主要的有以下幾個方面：偏最小二乘回歸是一種多因變量對多自變量的回歸建模方法。偏最小二乘回歸可以較好地解決許多以往用普通多元回歸無法解決的問題。在普通多元線形回歸的應用中，我們常受到許多限制C最典型的問題就是自變量之間的多重相關性。如果采用普通的最小二乘方法，這種變量多重相關性就會嚴重危害參數(shù)估計，擴大模型誤差，并破壞模型的穩(wěn)定性。變量多重相關問題十分復雜，長期以來在理論和方法上都未給出滿意的答案，這一直困擾著從事實際系統(tǒng)分析的工作人員。在偏最小二乘回歸中開辟了一種有效的技術途徑，它利用對系統(tǒng)中的數(shù)據(jù)信息進行分解和篩選的方式，提取對因變量的解釋性最強的綜合變量，辨識系統(tǒng)中的信息與噪聲，從而更好地克服變量多重相關性在系統(tǒng)建模中的不良作用。偏最小二乘回歸之所以被稱為第二代回歸方法，還由于它可以實現(xiàn)多種數(shù)據(jù)分析方法的綜合應用。偏最小二乘回歸=多元線性回歸分析+典型相關分析+主成分分析由于偏最小二乘回歸在建模的同時實現(xiàn)了數(shù)據(jù)結(jié)構的簡化，因此，可以在二維平面圖上對多維數(shù)據(jù)的特性進行觀察，這使得偏最小二乘回歸分析的圖形功能十分強大。在一次偏最小二乘回歸分析計算后，不但可以得到多因變量對多自變量的回歸模型，而且可以在平面圖上直接觀察兩組變量之間的相關關系，以及觀察樣本點間的相似性結(jié)構。這種高維數(shù)據(jù)多個層面的可視見性，可以使數(shù)據(jù)系統(tǒng)的分析內(nèi)容更加豐富，同時又可以對所建立的回歸模型給予許多更詳細深入的實際解釋。一、偏最小二乘回歸的建模策略原理方法

1.1建模原理設有q個因變量｛y】,...，｝和p自變量｛xi,...,xP｝。為了研究因變量和自變量yq的統(tǒng)計關就們9!n個樣本獨此構成了自變量與因變量的數(shù)據(jù)表X=｛X1,...,X］,｝和.Y=｛yi,...,yj。偏最小二乘回歸分鐘X與Y中提取出成分ti和U］(也就是說t】是x】，...,xp的線形繾u】是貝，...，無的線形捶在提取這兩個成分時，為了回歸分析的需，有下列兩個要求：t】和U】應盡可能大地攜帶他們各自數(shù)據(jù)表中的變篇息；t】與U】的相關程度能夠達到最大。這兩個要求表明，t和Ui應盡可能好的代表數(shù)據(jù)表X和Y,同時自變量的成、對因變量的成分U】又有最強的解釋能力在第一個成分t］和U】被提取后，偏最小二乘回歸分荊施對t】的回歸以及Y對u的回歸。如果回歸方程已經(jīng)達到滿意的精度，則算法遂屈將利用X被t】解釋后的殘余信息以及Y被t】解釋后的殘余信息進行第二前成分提取。如此往復，直到能達到一個較滿意的精度為止。鼓咨X共提取了m個成分t,?-,、，偏最小二乘回歸將通規(guī)y對，…，,1%的回歸，然后再表達成v關于原變童xlu,的回歸方k程，k=l,2，…,qo1,xX1.2計算方法推導*為了數(shù)學推導方便先將數(shù)據(jù)做標準化處理。X經(jīng)標準化處理后的數(shù)據(jù)矩隱e=(e°，，E。)n,Yj經(jīng)標準化處理后的數(shù)據(jù)矩隔0IpPF=(F(n,nEQ的第一個成分，W］是E伊第一個軸，它是一個單位向量,第一步記是］「手)°P既I|wi|EQ的第一個成分，W］是E伊第一個軸，它是一個單位向量,第一步記是］u是Fo的第一個成分，U1=Foc】點。2Cov(tu)=fv^-()iFt/t與Cov(tu)=fv^-()iFt/t與u】的協(xié)方差達到最大,年U)T，max析原理，應該有Var(Ui)->maxVar(tiImax另一方面，由于回歸建模的需要，又要求t】對U］有很大的解釋能力，有典型相關分析的思路，t】與％的相關度應達到最大值，既r(tl,因此，綜合起來，在偏最小二乘回歸中，我們要求值。正規(guī)的數(shù)學表述應亥是求偷下列優(yōu)化問題)既maxcw：1s.twclw11'c11C12=1和I|c因此，將在I|wI|如果采用拉格朗日算條，記?8.對S分別求關于W「C1,e1,as=EW12=1(~的約京條伸下，去求W1MEW1(Wi1-1)-2的偏導并令之為零，有22)Wi=000%)的最F大E(clc

1-D(1-Ci=0(1-s1「=FoEoWICi=0(1-s1=-(wiWi—1)=04)(1-一責料分字一(1-5)由式(1-2)?(1-5),可以推出wSFcEw,Fc1°010101，所以，】正是優(yōu)化問題的目標函數(shù)=0把式(1-2)和式(1-3)寫成=0(1-6)FEw°01(1-7)將式(1-7)代入式(1-6),有E'FF°0(1-8)Ew°0同理，可得(1-9)F'EE°0(1-9)。的特征向量，對蒞的特征值為1?】是目標函數(shù)值，它要。的特征向量，對蒞的特征值為1?】是目標函數(shù)值，它要可見,W］是矩陣E00求取最大值,所以，w是對應于E=0七矩陣最大特征值的單位特征向量.而另一方面，C】是對應于矩陣EE,一方面，C】是對應于矩陣00最大特征值求得軸w和c】后，即可得到成分1t】％W]然后,分別求E和F°對ti,u】的三個回歸財呈Eotp

1

9F°uqiIFotI*】式中，回歸系數(shù)向量是4E(1-10)1F](1-11)F(1-12)1°1(1-14)lluII乙

1°1(1-15)IItII

乙1而分別是三個回歸方程的殘差矩陣?0=<E和Fi取E和氏，然后，求第二個軸W2和c2以及第二步用殘差矩陣二個成分t,u2,W是對應于矩陣t=ElW22=u=F]c2

乙EFc1121:最大雋征值FFE11的特征值，5是對應于矩陣最大特征值的特征向量.計算回歸系數(shù)FP因此，有回歸方程IItII

乙2E】Fi如此計算下去，如果X的秩是A,則會有EotP(1-16)tAPAFotr1tAraFAA(1-17)由于,ti,,tA均可以表示成Eoi,,E°p的線性組合，因此，式(1-17)還可以還原成制關Xj*=E?的回歸方程形式，即ykFf=a+…+a+**yk*XXFk=l,2,…,qkl11&AkF曲是殘差距隔的第k列。Ak1.3交叉有效性…下面要討論的問題是在現(xiàn)有的數(shù)據(jù)表，如何確定更好的回歸方程。在誨情形了，偏最小二乘回歸方程并不需要選用全部的成分t】，，tA進行回歸建模，而是可以象在主成分分析一樣，采用截尾的方式選捕m個成分(mA,A秩(X)),僅用這m個后續(xù)的成分t】，，捻就可以得到一個預測性枚再用的模型。事實上，如果后續(xù)的成分已經(jīng)不能為解F提街更有意義的信過多的成分只會破壞對統(tǒng)藉,引導錯誤的預測結(jié)論C在多元回歸分析一章中，我們曾在調(diào)整復測定系數(shù)的內(nèi)容中討論過這一觀。再用下面的問題是怎樣來確定所應提取的成分個數(shù)。在多元回歸分析中，曾介紹過點抽樣測試法來確定回歸模型是否適預測應用。我們把手中的數(shù)據(jù)分成兩部分O：第一部分用于建立回歸方程，求出回歸系數(shù)估量擬臺簿以零碧割方和？:;再用第二部分數(shù)據(jù)作為實賽，代入剛才所求得的回歸方程，由此求出殉和*。f地，若有4?B，則回歸方程會有更好的預測效果。若？'？乂TB，則回歸方程不宜用于預測。在偏最小二乘回歸建模中，究竟應該選取多少個成分為宜，這可通過考察增一個新的成分后，能否對模型的預測功能有明顯的改進來考慮。采用新瞞測試法的工作方式，把所有n個樣本點分成兩部分：第一部分除去某個樣本點i的所有樣本點集合(共含n'l個樣本點)，用這部分樣本點并用h個成分擬合一個回歸方程；第二部分是把剛才被排除的樣本點i代入前面擬合的回歸方程，得到y(tǒng)3在樣本點i上的擬合值街(阿〒每艾個i三L2,二，n,重復上述測試，則可以定義yj的預測誤差弄和為PRESS，有(1-18)PRESSS(yi11](1-18)PRESSh定義Y的預測誤差平方和為PRESS^(1-19)顯然,如果回歸方程的穩(wěn)健性不好，誤差就很大，它對樣本點的變動就會十分敏感，這種擾動誤差的作用，就會加大PRESS.的值。PRESSh另外,再采用所有的樣本點，擬合含h個成分的回歸方程。這是，記第i個樣本點的預測值為反而］則可以記y的誤差平方和為，，坷，-n2SShj-（y如.）（1-20）=々hji定義Y的誤差平方和為SS」,有h--_PSShSg.（1-21）j1'一般說來,總是有PRESS大于SSh,而SSh則總是小于SSh10下面比較SSh1和，PRESq。SS-是用全部樣本點擬合的具有h-1個成分的方程的擬合誤差PRESSS加了一個成分加,但卻含有樣本，秘勺擾動誤差。如果h個成分的回歸方h一程的/含擾動誤差能在#^荏度上小于（h-1）個成分回歸方程的擬合誤差，則認為增羽一個成分〉t，會罪預四結(jié)果明顯提高。因此我們希望（PRESS、/SSh1）的比傕能一越小越好。在SIMCA-P軟件中，指定（PRESS】】（PRESS】】/SS】】z0.95PRESS】】0.95即PRESSh0.95SSh時，增加曲分捉就是有益的；或者反過來說，當時,就認為增加新的成分h,對減少方程的預測誤差無明顯PRESS】】0.95另有一種等價的定義稱為交叉有效性。對每一個變量尸卜,定義PRESS2122)Qhk(1-hkSS、(hl)k7對于全部因變量Y,成分%交叉有效性定光去PRESS=nl<21^1-Q對于全部因變量Y,成分%交叉有效性定光21^1-QhSS>一=(hl)kPRESS1一h(1-23)SS(h1)用交叉有效性測量成分'對預測模型精度的邊際貢獻有如南。n22⑴當(10.95)0.0975時，也成分的邊際貢獻是顯著的顯而易見Qu0.0975與(PRESSh/SShJ0.95是完全等價的決策原則。n(2)對于k=l,2,…，q,至少有一個虹使得214

Q'?號這時增加成分至少使一個因變量％的預測模型得到顯著的善因此，也n可以考慮增加成分％是明顯有益的。明確了偏最小二乘回歸方法的基本原理方法及算法步驟后我們將做實證分析附錄functionw=maxdet(A)探矩陣的最大特值[v,d]=eig(A);[n,p]=size(d);dl=d*ones(p,l)；d2=max(dl);8i=find(dl==d2);w=v(:,i)；%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function[c,m,v]=norml(C)渤數(shù)據(jù)進行標準化處理[n,s]=size(C);fori=l:nforj=l:sc(i,j)=(C(i,j)-mean(C(:,j)))/sqrt(cov(C(:,j)));endendm=mean(C);forj=l:sv(l,j)=sqrt(cov(C(:,j)));end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function[t,q,w,wh,fO,FF]=fun717(px,py,C)%px自變量的輸入個數(shù)%py輸入因變量的個數(shù)c%C輸入的自變量和因變量組成的矩陣%t提取的主成分%q為回歸系數(shù)。%w最大特征值所對應的特征向量。一費料分耳一%wh處理后的特征向量%fO回歸的標準化的方程系數(shù)%FF原始變量的回歸方程的系數(shù)c=norml(C);%norml為標準化函數(shù)y=c(:,px+l:px+py);截取標準化的因變量EO=c(:,l:px);FO=c(:,px+1:px+py);A=EO,*FO*FO,*EO;求最大特征向量w(:,l)=maxdet(A);提取主成分E(:,l:px)=EO-t(:,l)*(ECT*t(:,l)/(t(:,l)?*t(:,l)))?；%%wh處理后的特征向量%fO回歸的標準化的方程系數(shù)%FF原始變量的回歸方程的系數(shù)c=norml(C);%norml為標準化函數(shù)y=c(:,px+l:px+py);截取標準化的因變量求最大特征向量fori=0:px-2B(:?px*i+l:px*i+px)=E(:,px*i+l:px*i+px),*FO*FO'*E(:,px*i+l:px*i+px)w(:,i+2)=maxdet(B(:,px*i+l:px*i+px));%maxdet為求最大特征值的函數(shù)

t(:,i+2)=E(:,px*i+l:px*i+px)*w(:,i+2);p(:,px*i+px+l:px*i+2*px)=(E(:,px*i+l:px*i+px),*t(:,i+2)/(t(:,i+2)'*t(E(:,px*i+px+l:px*i+2*px)=E(:,px*i+l:px*i+px)-t(:,i+2)*(E(:,px*i+l:px*i+px)'*t(:,i+2)/(t(:,i+2)'*t(:,i+2)))‘;endfors=l:px求回歸系數(shù)%noq(:,s)=p(l,px*(s-l)+l:px*s)r;求回歸系數(shù)%noend[n,d]=size(q);forh=l:pxiw=eye(d);forj=l:h-liw=iw*(eye(d)-w(:,j)*q(：,j)')；endwh(:,h)=iw*w(:,h);endforj=l：pyzr(j,:)=(regress1(y(:,j),t)),；%endforj=l:pxfori=l:py%生成標準化變量的方程的系數(shù)矩陣wl=wh(:,l-j)；zrl=(zr(i,l:j))';fO(i,:,j)=(wl*zrl)';end[normxy,meanxy,covxy]=norm1(C);rmxy標準化后的數(shù)據(jù)矩陣%meanxy每一列的均值%covxy每一列的方差ccxx=ones(py,1)*meanxy(1,l:px);ccy=(covxy(l,px+1:px+py))'*ones(1,px);ccx=ones(py,l)*(covxy(1,1:px));ff=ccy.*fO(:,:,j)./ccx;fff=-(sum((ccy.*ccxx?*fO(:,:,j)./ccx)')-meanxy(l,px+l:px+py))';FF(:,:,j)=[fff,ff];%生成原始變量方程的常數(shù)項和系數(shù)矩陣end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function[r,Rdyt,RdYt,RdYtt,Rdytt?VIP]=fun8y(px,py,c)X=c(:,l:px);Y=c(:,px+1:px+py);x=norml(X);y=norml(Y);[t,q,w]=fun717(px,py,[X,Y]);rl=corrcoef([y,t]);r=rl(py+1:px+py,1:py)1；Rdyt=r.A2;RdYt=mean(Rdyt)form=l:pxRdYtt(1,m)=sum(RdYt(l,1endforj=l：pyform=l:pyRdytt(j,m)=sum(Rdyt(j,endendforj=l:pxform=l:pxRd(j,m)=RdYt(l,l:m)*((w(j,l:m).A2)')；endendforj=l:pxVIP0,:)=sqrt((px*ones(1,px)./RdYtt).*Rd0,:));end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function[r,Rdxt,RclXt,RdXtt,Rdxtt]=fun8x(px,py,c)X=c(:,l:px);Y=c(:,px+1:px+py);x=norml(X);y=norml(Y);[t,q,w]=fun717(px,py,[X,Y]);rl=corrcoef([x4]);r=rl(px+1:px+px,1:px)';Rdxt=r.A2;RclXt=mean(Rdxt);form=l:pxRclXtt(1,m)=sum(RdXt(1,1：m)')；endforj=l:pxform=l:pxRdxtt(j,m)=sum(Rdxt(j,1:m));endend%forj=l:px%form=l:px%Rd(j,m)=RdXt(l,l:m)*((w(j,1:m).A2)')；%end%end%forj=l:px%VIP(j,:)=sqrt((px*ones(l,px)./RdYtt).*Rd(j,:));%end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function[t,u]=TU(px,py,C)%t提取的自變量的主成分%U提取的因變量的主成分c=norml(C);y=c(:,px+l:px+py);EO=c(:,l:px);FO=c(:,px+1:px+py);A=EO'*FO*FO'*EO;w(:,1)=maxdet(A);t(:,l)=EO*w(:,l)；B=FO'*EO*EO'*FO;cc(:,l)=maxdet(B);u(:,l)=FO*cc(:,l)；%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%functiondrew(px,py,c)X=c(:,l:px);Y=c(:,px+1:px+py);[line,l]=size(Y);[t,q,w,wh,fO,FF]=fun717(px,py,c);YY=X*FF(:,2:px+l,3)'+ones(line,l)*FF(:,l,3)';14subplot(l,1,1,1)bar(f0(:,:,3))titlef直方圖')legend('SG','TZBFB,'FHL','JK','HPZD','JPZD'