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

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

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

1-D(1-Ci=0(1-s1「=FoEoWICi=0(1-s1=-(wiWi—1)=04)(1-一責(zé)料分字一(1-5)由式(1-2)?(1-5),可以推出wSFcEw,Fc1°010101，所以，】正是優(yōu)化問(wèn)題的目標(biāo)函數(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)。的特征向量，對(duì)蒞的特征值為1?】是目標(biāo)函數(shù)值，它要。的特征向量，對(duì)蒞的特征值為1?】是目標(biāo)函數(shù)值，它要可見(jiàn),W］是矩陣E00求取最大值,所以，w是對(duì)應(yīng)于E=0七矩陣最大特征值的單位特征向量.而另一方面，C】是對(duì)應(yīng)于矩陣EE,一方面，C】是對(duì)應(yīng)于矩陣00最大特征值求得軸w和c】后，即可得到成分1t】％W]然后,分別求E和F°對(duì)ti,u】的三個(gè)回歸財(cái)呈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而分別是三個(gè)回歸方程的殘差矩陣?0=<E和Fi取E和氏，然后，求第二個(gè)軸W2和c2以及第二步用殘差矩陣二個(gè)成分t,u2,W是對(duì)應(yīng)于矩陣t=ElW22=u=F]c2

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

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

Q'?號(hào)這時(shí)增加成分至少使一個(gè)因變量％的預(yù)測(cè)模型得到顯著的善因此，也n可以考慮增加成分％是明顯有益的。明確了偏最小二乘回歸方法的基本原理方法及算法步驟后我們將做實(shí)證分析附錄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ù)進(jìn)行標(biāo)準(zhǔn)化處理[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自變量的輸入個(gè)數(shù)%py輸入因變量的個(gè)數(shù)c%C輸入的自變量和因變量組成的矩陣%t提取的主成分%q為回歸系數(shù)。%w最大特征值所對(duì)應(yīng)的特征向量。一費(fèi)料分耳一%wh處理后的特征向量%fO回歸的標(biāo)準(zhǔn)化的方程系數(shù)%FF原始變量的回歸方程的系數(shù)c=norml(C);%norml為標(biāo)準(zhǔn)化函數(shù)y=c(:,px+l:px+py);截取標(biāo)準(zhǔn)化的因變量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回歸的標(biāo)準(zhǔn)化的方程系數(shù)%FF原始變量的回歸方程的系數(shù)c=norml(C);%norml為標(biāo)準(zhǔn)化函數(shù)y=c(:,px+l:px+py);截取標(biāo)準(zhǔn)化的因變量求最大特征向量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%生成標(biāo)準(zhǔn)化變量的方程的系數(shù)矩陣wl=wh(:,l-j)；zrl=(zr(i,l:j))';fO(i,:,j)=(wl*zrl)';end[normxy,meanxy,covxy]=norm1(C);rmxy標(biāo)準(zhǔn)化后的數(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ù)項(xiàng)和系數(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'