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

1、偏最小二乘回歸是一種新型的多元統(tǒng)計數(shù)據(jù)分析方法，它與1983年由伍德 和阿巴諾等人首次提出。 近十年來， 它在理論、 方法和應(yīng)用方面都得到了迅速的發(fā)展。密西根大學(xué)的弗耐爾教授稱偏最小二乘回歸為第二代回歸分析方法。偏最小二乘回歸方法在統(tǒng)計應(yīng)用中的重要性主要的有以下幾個方面： (1)偏最小二乘回歸是一種多因變量對多自變量的回歸建模方法。 (2)偏最小二乘回歸可以較好地解決許多以往用普通多元回歸無法解決的問題。 在普通多元線形回歸的應(yīng)用中，我們常受到許多限制。最典型的問題就是自變量之間的多重相關(guān)性。如果采用普通的最小二乘方法，這種變量多重相關(guān)性就會嚴(yán)重危害參數(shù)估計，擴(kuò)大模型誤差，并破壞模型的穩(wěn)定性。

2、變量多重相關(guān)問題十分復(fù)雜，長期以來在理論和方法上都未給出滿意的答案，這一直困擾著從事實(shí)際系統(tǒng)分析的工作人員。在偏最小二乘回歸中開辟了一種有效的技術(shù)途徑，它利用對系統(tǒng)中的數(shù)據(jù)信息進(jìn)行分解和篩選的方式，提取對因變量的解釋性最強(qiáng)的 綜合變量，辨識系統(tǒng)中的信息與噪聲，從而更好地克服變量多重相關(guān)性在系統(tǒng)建模中的不良作用。 (3)偏最小二乘回歸之所以被稱為第二代回歸方法， 還由于它可以實(shí)現(xiàn)多種數(shù)據(jù)分析方法的綜合應(yīng)用。 偏最小二乘回歸=多元線性回歸分析+典型相關(guān)分析+主成分分析 由于偏最小二乘回歸在建模的同時實(shí)現(xiàn)了數(shù)據(jù)結(jié)構(gòu)的簡化，因此，可以在二 維平面圖上對多維數(shù)據(jù)的特性進(jìn)行觀察，這使得偏最小二乘回歸分析的

3、圖形功能十分強(qiáng)大。在一次偏最小二乘回歸分析計算后，不但可以得到多因變量對多自變量的回歸模型，而且可以在平面圖上直接觀察兩組變量之間的相關(guān)關(guān)系，以及觀 察樣本點(diǎn)間的相似性結(jié)構(gòu)。這種高維數(shù)據(jù)多個層面的可視見性，可以使數(shù)據(jù)系統(tǒng)的分析內(nèi)容更加豐富，同時又可以對所建立的回歸模型給予許多更詳細(xì)深入的實(shí)際解釋。 偏最小二乘回歸的建模策略原理方法 1.1 建模原理 設(shè)有q個因變量yi,.,yq和p自變量x1,.,Xp。為了研究因變量和自變量的統(tǒng)計關(guān)系，我們觀測了n個樣本點(diǎn)，由此構(gòu)成了自變量與因變量的數(shù)據(jù)表X=Xi,.,Xp和.Y=yi,.,yq 。偏最小二乘回歸分別在X與Y中提取出成分 ti和 Ui(也就是說

4、， ti是 Xi，.,Xp的線形組合，Ui是 yi,.,yq的線形組合).在提取這兩個成分時，為了回歸分析的需要，有下列兩個要求： (1) ti和 Ui應(yīng)盡可能大地攜帶他們各自數(shù)據(jù)表中的變異信息； (2) ti與 Ui的相關(guān)程度能夠達(dá)到最大。 這兩個要求表明，ti和 Ui應(yīng)盡可能好的代表數(shù)據(jù)表X和Y,同時自變量的成分 ti對因變量的成分 Ui又有最強(qiáng)的解釋能力。 在第一個成分 ti和 Ui被提取后，偏最小二乘回歸分別實(shí)施X對 ti的回歸 以及Y對 Ui的回歸。如果回歸方程已經(jīng)達(dá)到滿意的精度，則算法終止；否則，將利用X被ti解釋后的殘余信息以及Y被ti解釋后的殘余信息進(jìn)行第二輪的成分提取。 如此

5、往復(fù)，直到能達(dá)到一個較滿意的精度為止。若最終對X共提取 了m個成分 titm偏最小二乘回歸將通過實(shí)施yk對 titm的 ， 回歸然后再表達(dá)成yk關(guān)于原變量Xi.Xm的回歸方程，k=i,2,qo 1.2 計算方法推導(dǎo) 為了數(shù)學(xué)推導(dǎo)方便起見，首先將數(shù)據(jù)做標(biāo)準(zhǔn)化處理。X經(jīng)標(biāo)準(zhǔn)化處理后的數(shù)據(jù)矩陣記為 E0=(E01,，Eop)n冷，Yj經(jīng)標(biāo)準(zhǔn)化處理后的數(shù)據(jù)矩陣記為 Fo=(Foi,，F(xiàn)oq)nxp0 第一步記 ti是Eo的第一個成分，Wi是Eo的第一個軸，它是一個單位向量, 既|阿|二1。 記 Ui是 Fo的第一個成分，Ui=FoCi。Ci是 Fo的第一個軸，并且11cl11=1。如果要 tiUi能分

6、別很好的代表X與Y中的數(shù)據(jù)變異信息，根據(jù)主成分分 析原理，應(yīng)該有 Var(u1)max Var(ti)，max 另一方面，由于回歸建模的需要，又要求 ti對 Ui有很大的解釋能力，有典型相關(guān)分析的思路，3 與 ui的相關(guān)度應(yīng)達(dá)到最大值，既 r(tiui)tmax 因此，綜合起來，在偏最小二乘回歸中，我們要求 ti與 Ui的協(xié)方差達(dá)到最大，既 UiqVar(tJVar(uJr(ti,ui)tmax 正規(guī)的數(shù)學(xué)表述應(yīng)該是求解下列優(yōu)化問題，既 max wiE。wi,FCi s.t 因此，將在|wi|2=i和11cl|2=i的約束條件下，去求(w,EF。)的最大 值。 如果采用拉格朗日算法，記 s=w

7、iE。F。一九1(Wiwi1)九2(CiCi-1) 對s分別求關(guān)于 wiCi人和九2的偏導(dǎo)并令之為零，有 -s, (wiwi1) -iCov(ti s mEF。Ci2九iwi-。 (1-2) 丁。&也-22C1=Q (1-3) =-(oci-1)=0(1-5) 由式(1-2)(1-5)，可以推出 一一一一 2.12,.2=w1E0F0C1：.E0wi,FoC2 記&=2%=2 九2=w1E0F0C1,所以,由正是優(yōu)化問題的目標(biāo)函數(shù)值. 把式(1-2)和式(1-3)寫成一, E0F0C1二1w(1-6) 一,F0E0W1=u1C1(1-7) 將式(1-7)代入式(1-6)，有LL-2 E0F0F

8、0E0w1=w1(1-8) 同理，可得 一一一一.2 F0E0E0F0C1-川G(1-9) 可見,W1是矩陣 E0F0F0E0的特征向量，對應(yīng)的特征值為 32.3 是目標(biāo)函數(shù)值，它要 求取最大值，所以，W1是對應(yīng)于 E0F0F0E0矩陣最大特征值的單位特征向量.而另 24 一萬面，C1是對位于矩陣 F0E0E0F0最大特征值備的單位特征向量. 求得軸 W1和 g 后，即可得到成分 t-E0W1 U1=F0C1 然后,分別求 E。和 F。對吊 3 的三個回歸方程 _1._E =t1P1E1(1-10) F0=uq1F1(1-11) F0=F1(1-12) 式中，回歸系數(shù)向量是 E0ti Pi2

9、l|ti|2 qi二 F0Ui 2 llull2 ri 一 F0tl 2 lltill2 而 Ei,F*i,Fi分別是三個回歸方程的殘差矩陣 第二步用殘差矩陣 Ei和 Fi取彳tE0和 Fo，然后，求第二個軸 個成分 t2,u2,有 t2=E1w2 U2=F1C2 2=；12,u2=w2E1F1c2 一，._9 W2是對位于矩陣 EiF1FiEi最大特征值外的特征值 .、. FiEiEiFi最大特征值的特征向量.計算回歸系數(shù) (i-i3) (i-i4) (i-i5) w2和 c2以及第二 C2是對應(yīng)于矩陣 p2= Eit2 , Fit2 2 11t2ll2 因此，有回歸方程 Ei=t2P2E2

10、 Fi=t2r2F2 如此計算下去，如果X的秩是A,則會有 E=tiPi-tAPA F0=tiitArA-FA (i-i6) (i-i7) 由于，ti,tA均可以表示成 E i，E0P的線性組合，因止匕，式(i-i7)還可以還原 成yk=Fok關(guān)于 Xj*=Eok的回歸萬程形式，即 yk*=kiXi,一：kpXpFAkk=1,2,q FAk是殘差距陣 FA的第k歹I。 1.3交叉有效性 下面要討論的問題是在現(xiàn)有的數(shù)據(jù)表下， 如何確定更好的回歸方程。 在許多情形下，偏最小二乘回歸方程并不需要選用全部的成分 ti,，tA進(jìn)行回3建模，而 是可以象在主成分分析一樣，采用截尾的方式選擇前m個成分(mA

11、,A=秩(X),僅用這m個后續(xù)的成分L,，tm就可以得到一個預(yù)測性較好的模型。事實(shí)上，如果后續(xù)的成分已經(jīng)不能為解釋 Fo提供更有意義的信息時，采用過多的成分只會破壞對統(tǒng)計趨勢的認(rèn)識，引導(dǎo)錯誤的預(yù)測結(jié)論。在多元回歸分析一章中，我們曾在調(diào)整復(fù)測定系數(shù)的內(nèi)容中討論過這一觀點(diǎn)。 下面的問題是怎樣來確定所應(yīng)提取的成分個數(shù)。 在多元回歸分析中，曾介紹過用抽樣測試法來確定回歸模型是否適于預(yù)測應(yīng)用。我們把手中的數(shù)據(jù)分成兩部分：第一部分用于建立回歸方程，求出回歸系數(shù)估計量 bB,擬合值？B以及殘差均方和屯；再用第二部分?jǐn)?shù)據(jù)作為實(shí)驗點(diǎn)，代入剛才所求得的回歸方程，由此求出外和封。一般地,若有解宋虎，則回歸方程會有更

12、好的預(yù)測效果。若夕;夕；,則回歸方程不宜用于預(yù)測。 在偏最小二乘回歸建模中，究竟應(yīng)該選取多少個成分為宜，這可通過考察增加一個新的成分后，能否對模型的預(yù)測功能有明顯的改進(jìn)來考慮。采用類似于抽樣測試法的工作方式，把所有n個樣本點(diǎn)分成兩部分：第一部分除去某個樣本點(diǎn)i的所有樣本點(diǎn)集合(共含n-1個樣本點(diǎn))，用這部分樣本點(diǎn)并使用h個成分?jǐn)M合一個回歸方程;第二部分是把剛才被排除的樣本點(diǎn)i代入前面擬合的回歸方程將到 yj在樣本點(diǎn)i上的擬合值(與。對于每一個i=1,2,n,重復(fù)上述測試，則可以定義 yj的預(yù)測誤差平方和為 PRESS#,有 n PRESSj八(yj-？川(4)2(1-18) i1Qhk=1 P

13、RESS S&hJ)k 定義Y的預(yù)測誤差平方和為 PRESSh，有 p PRESS=PRESShj(1-19) j4 顯然,如果回歸方程的穩(wěn)健性不好，誤差就很大，它對樣本點(diǎn)的變動就會十分敏感， 這種擾動誤差的作用，就會加大 PRESSh的值。 另外,再采用所有的樣本點(diǎn)，擬合含h個成分的回歸方程。這是，記第i個樣本 點(diǎn)的預(yù)測值為編，則可以記 yj的誤差平方和為 SShj,有n SSj八(yij-yhji)2(1-20) i4 定義Y的誤差平方和為 SSh，有 p SSSSj(1-21) j1 一般說來,總是有 PRESSh大于 SSh,而 SSh則總是小于 SSh。下面比較 SSh和 PRESS

14、hS&是用全部樣本點(diǎn)擬合的具有h-1個成分的方程的擬合誤差 PRESSh增加了一個成分 th，但卻含有樣本點(diǎn)的擾動誤差。如果h個成分的回歸方程的含擾動誤差能在一定程度上小于(h-1)個成分回歸方程的擬合誤差， 則認(rèn)為增加一個成分th,會使預(yù)測結(jié)果明顯提高。因此我們希望(PRESSh/S5)的比值能越小越女?0在SIMCA-P軟件中指定 (PRESSh/SShj)0.952 即JPRESS0.95：SS、時，增加成分 th就是有益的；或者反過來說，當(dāng)PRESSA0.95VSS7時,就認(rèn)為增加新的成分 th,對減少方程的預(yù)測誤差無明顯的改善作用. 另有一種等價的定義稱為交叉有效性。對每一個變量 y

15、一定義 對于全部因變量Y,成分 th交叉有效性定義為 用交叉有效性測量成分 th對預(yù)測模型精度的邊際貢獻(xiàn)有如下兩個尺度。 當(dāng) Q：（10.952）=0.0975 時，th成分的邊際貢獻(xiàn)是顯著的。顯而易 見，Q2圭 0.0975 與（PRESSh/SSh）0.952是完全等價的決策原貝限 （2）對于k=1,2,q,至少有一個k,使得 Q；_0.0975 這時增加成分 th,至少使一個因變量 yk的預(yù)測模型得到顯著的改善， 因此， 也可以考慮增加成分 th是明顯有益的。 明確了偏最小二乘回歸方法的基本原理、方法及算法步驟后，我們將做 實(shí)證分析。 附錄 functionw=maxdet(A) %求矩

16、陣的最大特征值 v,d=eig(A); n,p=size(d); d1=d*ones(p,1); d2=max(d1); “PRESShk :jPRESS .二SS；h）k SS（h） (1-23) i=find(d1=d2); w=v(：,i); % % functionc,m,v=norm1(C) %對數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理 n,s=size(C); fori=1:n forj=1:s c(i,j)=(C(i,j)-mean(C(:,j)/sqrt(cov(C(:,j); end end m=mean(C); forj=1:s v(1,j)=sqrt(cov(C(:,j); end % % f

17、unctiont,q,w,wh,f0,FF=fun717(px,py,C) %px自變量的輸入個數(shù) %py輸入因變量的個數(shù)。 %C輸入的自變量和因變量組成的矩陣 %t提取的主成分 %q為回歸系數(shù)。 %w最大特征值所對應(yīng)的特征向量。 %wh處理后的特征向量 %f0回歸的標(biāo)準(zhǔn)化的方程系數(shù) %FF原始變量的回歸方程的系數(shù)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(

18、:,1)/(t(:,1)*t(:,1); 得回歸系數(shù) p(:,1:px)=(E0*t(:,1)/(t(:,1)*t(:,1); fori=0:px-2 B(:,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為求最大特征值的函數(shù)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)*

19、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); %norm1為標(biāo)準(zhǔn)化函數(shù) %截取標(biāo)準(zhǔn)化的因變量 %求最大特征向量 %提取主成分 end fors=1:px q(:,s)=p(1,px*(s-1)+1:px*s); end n,d=size(q); forh=1:px iw=eye(d); forj=1:h-1 iw=iw*(eye(d)-w(:,j)*q(:,j); end wh(:,h)=iw*w(:,h); en

20、d forj=1:py zr(j,:)=(regress1(y(:,j),t); end forj=1:px fori=1:py% 生成標(biāo)準(zhǔn)化變量的方程的系數(shù)矩陣 w1=wh(:,1:j); zr1=(zr(i,1:j); f0(i,:,j)=(w1*zr1); end normxy,meanxy,covxy=norm1(C);%n ormxy標(biāo)準(zhǔn)化后的數(shù)據(jù)矩陣 %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

21、(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; %求回歸系數(shù) 生成原始變量方程的常數(shù)項和系數(shù)矩陣 end % % functionr,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(p

22、y+1:px+py,1:py),; Rdyt=r.A2; RdYt=mean(Rdyt) form=1:px RdYtt(1,m)=sum(RdYt(1,1:m); end forj=1:py form=1:py Rdytt(j,m)=sum(Rdyt(j,1:m); endend forj=1:px form=1:px Rd(j,m)=RdYt(1,1:m)*(w(j,1:m).A2); end end forj=1:px VIP(j,:)=sqrt(px*ones(1,px)./RdYtt).*Rd(j,:); end % % functionr,Rdxt,RdXt,RdXtt,Rdxtt

23、=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);RdXtt(1,m)=sum(RdXt(1,1:m); end forj=1:px form=1:px Rdxtt(j,m)=sum(Rdxt(j,1:m); end end %forj=1:px %form=1:px %Rd(j,m)=RdXt(1,1:m)*(w(j

24、,1:m).A2); %end %end %forj=1:px %VIP(j,:)=sqrt(px*ones(1,px)./RdYtt).*Rd(j,:); %end % % functiont,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); % % functiondrew(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(