偏最小二乘回歸方法及其應(yīng)用

1、偏最小二乘回歸方法及其應(yīng)用王惠文著國防工業(yè)出版社1999年版偏最小二乘回歸多元線性回歸分析典型相關(guān)分析主成分分析與傳統(tǒng)多元線性回歸模型相比，偏最小二乘回歸的特點是：（1）能夠在自變量存在嚴(yán)重多重相關(guān)性的條件下進(jìn)行回歸建模；（2）允許在樣本點個數(shù)少于變量個數(shù)的條件下進(jìn)行回歸建模；（3）偏最小二乘回歸在最終模型中將包含原有的所有自變量；（4）偏最小二乘回歸模型更易于辨識系統(tǒng)信息與噪聲（甚至一些非隨機性的噪聲）；（5）在偏最小二乘回歸模型中，每一個自變量的回歸系數(shù)將更容易解釋。在計算方差和協(xié)方差時，求和號前面的系數(shù)有兩種取法：當(dāng)樣本點集合是隨機抽取得到時，應(yīng)該取1/(n-1)；如果不是隨機抽取的，這

2、個系數(shù)可取1/n。多重相關(guān)性的診斷1 經(jīng)驗式診斷方法1、在自變量的簡單相關(guān)系數(shù)矩陣中，有某些自變量的相關(guān)系數(shù)值較大。2、回歸系數(shù)的代數(shù)符號與專業(yè)知識或一般經(jīng)驗相反；或者，它同該自變量與y的簡單相關(guān)系數(shù)符號相反。3、對重要自變量的回歸系數(shù)進(jìn)行t檢驗，其結(jié)果不顯著。特別典型的是，當(dāng)F檢驗?zāi)茉诟呔认峦ㄟ^，測定系數(shù)R2的值亦很大，但自變量的t檢驗卻全都不顯著，這時，多重相關(guān)性的可能性將很大。4、如果增加（或刪除）一個變量，或者增加（或刪除）一個觀測值，回歸系數(shù)的估計值發(fā)生了很大的變化。5、重要自變量的回歸系數(shù)置信區(qū)間明顯過大。6、在自變量中，某一個自變量是另一部分自變量的完全或近似完全的線性組合。7

3、、對于一般的觀測數(shù)據(jù)，如果樣本點的個數(shù)過少，樣本數(shù)據(jù)中的多重相關(guān)性是經(jīng)常存在的。但是，采用經(jīng)驗式方法診斷自變量系統(tǒng)中是否確實存在多重相關(guān)性，并不十分可靠，另一種較正規(guī)的方法是利用統(tǒng)計檢驗（回歸分析），檢查每一個自變量相對其它自變量是否存在線性關(guān)系。2 方差膨脹因子最常用的多重相關(guān)性的正規(guī)診斷方法是使用方差膨脹因子。自變量xj的方差膨脹因子記為（VIF）j，它的計算方法為（4-5） （VIF）j =（1-R j2）-1式中，R j2是以xj為因變量時對其它自變量回歸的復(fù)測定系數(shù)。所有xj變量中最大的（VIF）j通常被用來作為測量多重相關(guān)性的指標(biāo)。一般認(rèn)為，如果最大的（VIF）j超過10，常常表示

4、多重相關(guān)性將嚴(yán)重影響最小二乘的估計值。（VIF）j被稱為方差膨脹因子的原因，是由于它還可以度量回歸系數(shù)的估計方差與自變量線性無關(guān)時相比，增加了多少。不妨假設(shè)x1,x2,xp均是標(biāo)準(zhǔn)化變量。采用最小二乘法得到回歸系數(shù)向量B，它的精度是用它的方差來測量的。B的協(xié)方差矩陣為Cov(B)= 2 (X'X)-1式中，2是誤差項方差。所以，對于回歸系數(shù)b j，有Var(b j)= 2cjjcjj是(X'X)-1矩陣中第j個對角元素。可以證明，cjj =（VIF）j嶺回歸分析1 嶺回歸估計量嶺回歸分析是一種修正的最小二乘估計法，當(dāng)自變量系統(tǒng)中存在多重相關(guān)性時，它可以提供一個比最小二乘法更為穩(wěn)

5、定的估計，并且回歸系數(shù)的標(biāo)準(zhǔn)差也比最小二乘估計的要小。根據(jù)高斯馬爾科夫定理，多重相關(guān)性并不影響最小二乘估計量的無偏性和最小方差性。但是，雖然最小二乘估計量在所有線性無偏估計量中是方差最小的，但是這個方差卻不一定小。于是可以找一個有偏估計量，這個估計量雖然有微小的偏差，但它的精度卻能夠大大高于無偏的估計量。在應(yīng)用嶺回歸分析時，它的計算大多從標(biāo)準(zhǔn)化數(shù)據(jù)出發(fā)。對于標(biāo)準(zhǔn)化變量，最小二乘的正規(guī)方程為rXXb=ryX式中，rXX是X的相關(guān)系數(shù)矩陣，ryX是y與所有自變量的相關(guān)系數(shù)向量。嶺回歸估計量是通過在正規(guī)方程中引入有偏常數(shù)c（c0）而求得的。它的正規(guī)方程為+（4-8） （rXX+ cI） bR=ryX

6、所以，在嶺回歸分析中，標(biāo)準(zhǔn)化回歸系數(shù)為（4-9） bR =（rXX+ cI）-1 ryX2 嶺回歸估計量的性質(zhì)（1）嶺回歸系數(shù)是一般最小二乘準(zhǔn)則下回歸系數(shù)的線性組合，即（4-10） bR =（I+ crXX-1）-1b（2）記是總體參數(shù)的理論值。當(dāng)0時，可以證明一定存在一個正數(shù)c0，使得當(dāng)0< c< c0時，一致地有（4-11） E| bR -|2 E| b -|2（3）嶺回歸估計量的絕對值常比普通最小二乘估計量的絕對值小，即（4-12） | bR |<| b |嶺回歸估計量的質(zhì)量取決于偏倚系數(shù)c的選取。c的選取不宜過大，因為E（bR）=（I+ crXX-1）-1 E （b）

7、=（I+ crXX-1）-1關(guān)于偏倚系數(shù)c的選取尚沒有正規(guī)的決策準(zhǔn)則，目前主要以嶺跡和方差膨脹因子為依據(jù)。嶺跡是指p-1個嶺回歸系數(shù)估計量對不同的c值所描繪的曲線（c值一般在01之間）。在通過檢查嶺跡和方差膨脹因子來選擇c值時，其判斷方法是選擇一個盡可能小的c值，在這個較小的c值上，嶺跡中的回歸系數(shù)已變得比較穩(wěn)定，并且方差膨脹因子也變得足夠小。從理論上，最佳的c值是存在的，它可以使估計量的偏差和方差的組合效應(yīng)達(dá)到一個最佳水準(zhǔn)。然而，困難卻在于c的最優(yōu)值對不同的應(yīng)用而有所不同，對其選擇還只能憑經(jīng)驗判斷。其他補救方法簡介最常見的一種思路是設(shè)法去掉不太重要的相關(guān)性變量。由于變量間多重相關(guān)性的形式十分

8、復(fù)雜，而且還缺乏十分可靠的檢驗方法，刪除部分多重相關(guān)變量的做法常導(dǎo)致增大模型的解釋誤差，將本應(yīng)保留的系統(tǒng)信息舍棄，使得接受一個錯誤結(jié)論的可能和做出錯誤決策的風(fēng)險都不斷增長。另一方面，在一些經(jīng)濟(jì)模型中，從經(jīng)濟(jì)理論上要求一些重要的解釋變量必須被包括在模型中，而這些變量又存在多重相關(guān)性。這時采用剔除部分相關(guān)變量的做法就不符合實際工作的要求。另一種補救的辦法是增加樣本容量。然而，在實際工作中，由于時間、經(jīng)費以及客觀條件的限制，增大樣本容量的方法常常是不可行的。此外，還可以采用變量轉(zhuǎn)換的方式，來削弱多重相關(guān)性的嚴(yán)重性。一階差分回歸模型有可能減少多重相關(guān)性的嚴(yán)重性。然而，一階差分變換又帶來了一些其它問題。

9、差分后的誤差項可能不滿足總體模型中關(guān)于誤差項不是序列相關(guān)的假定。事實上，在大部分情形下，在原來的誤差項是不自相關(guān)的條件下，一階差分所得到的誤差項將會是序列相關(guān)的。而且，由于差分方法損失了一個觀察值，這在小樣本的情況下是極不可取的。另外，一階差分方法在截面樣本中是不宜利用的。1 主成分分析主成分分析的計算結(jié)果必然受到重疊信息的影響。因此，當(dāng)人為地采用一些無益的相關(guān)變量時，無論從方向上還是從數(shù)量上，都會扭曲客觀結(jié)論。在主成分分析之前，對變量系統(tǒng)的確定必須是慎之又慎的。2 特異點的發(fā)現(xiàn)第i個樣本點（樣本量為n）對第h主成分的貢獻(xiàn)率是（5-32） CTR(i)=Fh2(i)/(nh) （若遠(yuǎn)超過1/n

10、，為特異點）3 典型相關(guān)分析從某種意義上說，多元回歸分析、判別分析或?qū)?yīng)分析等許多重要的數(shù)據(jù)分析方法，都可以歸結(jié)為典型相關(guān)分析的一種特例，同時它還是偏最小二乘回歸分析的理論基石。典型相關(guān)分析，是從變量組X中提取一個典型成分F=Xa，再從變量組Y中提取一個成分G=Yb，在提取過程中，要求F與G的相關(guān)程度達(dá)到最大。在典型相關(guān)分析中，采用下述原則尋優(yōu)，即max<F,G>=aX'Yb a'X'Xa=1, b'Y'Yb=1其結(jié)果為，a是對應(yīng)于矩陣V11-1 V12 V22-1 V21最大特征值的特征向量，而b是對應(yīng)于矩陣V22-1 V21V11-1 V

11、12最大特征值的特征向量，這兩個最大特征值相同。其中，V11=X'X，V12=X'Y，V22=Y'Y。F與G之間存在著明顯的換算關(guān)系。有時只有一個典型成分還不夠，還可以考慮第二個典型成分。多因變量的偏最小二乘回歸模型1 工作目標(biāo)偏最小二乘回歸分析的建模方法設(shè)有q個因變量和p個自變量。為了研究因變量與自變量的統(tǒng)計關(guān)系，觀測了n個樣本點，由此構(gòu)成了自變量與因變量的數(shù)據(jù)表X和Y。偏最小二乘回歸分別在X與Y中提取出t和u，要求：（1）t和u應(yīng)盡可能大地攜帶它們各自數(shù)據(jù)表中的變異信息；（2）t和u的相關(guān)程度能夠達(dá)到最大。在第一個成分被提取后，偏最小二乘回歸分別實施X對t的回歸以及

12、Y對t的回歸。如果回歸方程已經(jīng)達(dá)到滿意的精度，則算法終止；否則，將利用X被t解釋后的殘余信息以及Y被t解釋后的殘余信息進(jìn)行第二輪的成分提取。如此往復(fù)，直到能達(dá)到一個較滿意的精度為止。若最終對X共提取了多個成分，偏最小二乘回歸將通過施行yk對X的這些成分的回歸，然后再表達(dá)成yk關(guān)于原自變量的回歸方程。2 計算方法首先將數(shù)據(jù)做標(biāo)準(zhǔn)化處理。X經(jīng)標(biāo)準(zhǔn)化處理后的數(shù)據(jù)矩陣記為E0=( E01,E0p)n×p，Y的相應(yīng)矩陣記為F0=( F01,F0q)n×q。第一步 記t 1是E0的第一個成分，t 1= E0w1，w1是E0的第一個軸，它是一個單位向量，即| w1|=1。記u 1是F0的

13、第一個成分，u 1= F0c1，c1是F0的第一個軸，并且| c1|=1。于是，要求解下列優(yōu)化問題，即（7-1）記1= w1'E0'F0c1，即正是優(yōu)化問題的目標(biāo)函數(shù)值。采用拉格朗日算法，可得（7-8） E0'F0F0'E0w1=12 w1（7-9） F0'E0E0'F0c1=12 c1所以，w1是對應(yīng)于E0'F0F0'E0矩陣最大特征值的單位特征向量，而c1是對應(yīng)于F0'E0E0'F0矩陣最大特征值12的單位特征向量。求得軸w1和c1后，即可得到成分t 1= E0w1u 1= F0c1然后，分別求E0和F0對t

14、1的回歸方程（7-10） E0= t 1 p1'+ E1（7-12） F0= t 1r1'+ F1式中，回歸系數(shù)向量是（7-13） p1= E0' t 1/| t 1|2（7-15） r1= F0' t 1/| t 1|2而E1和F1分別是兩個方程的殘差矩陣。第二步 用殘差矩陣E1和F1取代E0和F0，然后，求第二個軸w2和c2以及第二個成分t2，u2，有t 2= E1w2u 2= F1c22=< t2, u2>= w2'E1'F1c2w2是對應(yīng)于E1'F1F1'E1矩陣最大特征值的單位特征向量，而c2是對應(yīng)于F1&#

15、39;E1E1'F1矩陣最大特征值22的單位特征向量。計算回歸系數(shù)p2= E1' t 2/| t 2|2r2= F1' t 2/| t2|2因此，有回歸方程E1= t 2 p2'+ E2F1= t 2r2'+ F2如此計算下去，如果X的秩是A，則會有（7-16） E0= t 1 p1'+t A pA'（7-17） F0= t 1r1'+ +t A rA'+ FA由于t1，t A均可以表示成E01，E0p的線性組合，因此，式（7-17）還可以還原成yk*= F0k關(guān)于xj*= E0j的回歸方程形式，即yk*=k1 x1*+k

16、p xp*+ FAk， k=1,2,qFAk是殘差矩陣FA的第k列。3 交叉有效性如果多一個成分而少一個樣本的預(yù)測誤差平方和（所有因變量和預(yù)測樣本相加）除以少一個成分的誤差平方和（所有的因變量和樣本相加）小于0.952，則多一個成分是值得的。4 一種更簡潔的計算方法用下述原則提取自變量中的成分t 1，是與原則式（7-1）的結(jié)果完全等價的，即（7-24）（1）求矩陣E0'F0F0'E0最大特征值所對應(yīng)的單位特征向量w1，求成分t 1，得t 1= E0w1E1= E0-t 1 p1'式中， p1= E0' t 1/| t 1|2（2）求矩陣E1'F0F0&#

17、39;E1最大特征值所對應(yīng)的單位特征向量w2，求成分t2，得t 2= E1w2E2= E1-t 2 p2'式中， p2= E1' t 2/| t2|2（m）至第m步，求成分tm= Em-1wm，wm是矩陣Em-1'F0F0'Em-1最大特征值所對應(yīng)的單位特征向量.如果根據(jù)交叉有效性，確定共抽取m個成分t1，tm可以得到一個滿意的觀測模型，則求F0在t1，tm上的普通最小二乘回歸方程為F0= t 1r1'+ +t mrm'+ Fm偏最小二乘回歸的輔助分析技術(shù)1 精度分析定義自變量成分th的各種解釋能力如下（1）th對某自變量xj的解釋能力（8-1）

18、 Rd(xj; th)=r2(xj, th)（2）th對X的解釋能力（8-2） Rd(X; th)=r2(x1, th) + + r2(xp, th)/p（3）t1，tm對X的累計解釋能力（8-3） Rd(X; t1，tm)= Rd(X; t1) + + Rd(X; tm)（4）t1，tm對某自變量xj的累計解釋能力（8-4） Rd(xj; t1，tm)= Rd(xj; t1) + + Rd(xj; tm)（5）th對某因變量yk的解釋能力（8-5） Rd(yk; th)=r2(yk, th)（6）th對Y的解釋能力（8-6） Rd(Y; th)=r2(y1, th) + + r2(yq, t

19、h)/q（7）t1，tm對Y的累計解釋能力（8-7） Rd(Y; t1，tm)= Rd(Y; t1) + + Rd(Y; tm)（8）t1，tm對某因變量yk的累計解釋能力（8-8） Rd(yk; t1，tm)= Rd(yk; t1) + + Rd(yk; tm)2 自變量x j在解釋因變量集合Y的作用x j在解釋Y時作用的重要性，可以用變量投影重要性指標(biāo)VIP j來測度VIP j 2=pRd(Y; t1) w1j2+ + Rd(Y; tm) wmj2/Rd(Y; t1) + + Rd(Y; tm)式中，whj是軸wh的第j個分量。注意 VIP1 2+ + VIP p2=p3 特異點的發(fā)現(xiàn)定義第i個樣本點對第h成分th的貢獻(xiàn)率Thi2，用它來發(fā)現(xiàn)樣本點集合中的特異點，即（8-