嶺回歸分析

胡良平

(1.軍事科學(xué)院研究生院，北京　100850；2.世界中醫(yī)藥學(xué)會聯(lián)合會臨床科研統(tǒng)計(jì)學(xué)專業(yè)委員會，北京　100029

1　概　　述

1.1　何為嶺回歸分析

嶺回歸分析是在構(gòu)建多重線性回歸模型時，對基于“最小二乘原理”推導(dǎo)出的估計(jì)回歸系數(shù)的計(jì)算公式作一下校正，使回歸系數(shù)更穩(wěn)定。

1.2　嶺回歸分析應(yīng)用的場合

當(dāng)自變量之間存在較強(qiáng)的多重共線性時，求得的多重線性回歸模型很不穩(wěn)定；尤其是某些自變量前回歸系數(shù)的正負(fù)號與實(shí)際問題的專業(yè)背景不吻合。此時，嶺回歸分析有可能較好地解決前述提及的問題。

1.3　嶺回歸分析的原理

多重線性回歸方程的回歸系數(shù)可以表示為

β=(X'X)-1X'Y

(1)

β(k)=(X'X+kIm)-1X'Y

(2)

即在矩陣X'X的主對角線元素上加上一個非負(fù)因子k，其中Im為m階單位矩陣，k>0稱為嶺參數(shù)或偏參數(shù)。如果k取與試驗(yàn)數(shù)據(jù)Y無關(guān)的常數(shù)，則β(k)為線性估計(jì)，否則β(k)為非線性估計(jì)。取不同的k，得到不同的嶺估計(jì)。故上式實(shí)際定義了一個很大的估計(jì)類。特別當(dāng)k=0時，β(k)=(X'X)X'Y就是β的最小二乘估計(jì)。當(dāng)k∈[0,+∞)時，對于每個i，β(k)的第i個分量βi(k)的值均為k的函數(shù)，在直角坐標(biāo)系中，點(diǎn)[k,βi(k)]所構(gòu)成的變化軌跡，稱為嶺跡。隨著k的增大，βi(k)的絕對值均趨于不斷變小[由于自變量之間可能存在相關(guān)性，個別βi(k)可能有小范圍的向上波動或改變正、負(fù)號]，若k→+∞，則β(k)→0。

1.4　如何實(shí)施嶺回歸分析

如果在進(jìn)行多重線性回歸分析時，從專業(yè)上或通過共線性診斷得知自變量間存在多重共線性，那么可以考慮用嶺回歸分析進(jìn)行參數(shù)估計(jì)，在進(jìn)行嶺估計(jì)時通常采用以下幾步：

第1步，嶺回歸分析通常要先對X變量作中心化和標(biāo)準(zhǔn)化處理，以使不同自變量處于同樣數(shù)量級上而便于比較。

第2步，確定k值。

(1)嶺跡法

嶺跡法主要是通過將β(k)的分量βi(k)的嶺跡畫在同一幅圖上，從圖中選擇盡可能小的k值，使得各回歸系數(shù)的嶺估計(jì)大體穩(wěn)定，即各分量在圖上的嶺跡曲線趨于平行于X軸。選擇k值的一般原則主要有：①各回歸系數(shù)的嶺估計(jì)基本穩(wěn)定；②用最小二乘估計(jì)時符號不合理的回歸系數(shù)，其嶺估計(jì)的符號將變得合理；③回歸系數(shù)的大小要與實(shí)際相符，即從專業(yè)上講對因變量影響較大的自變量其系數(shù)的絕對值也較大；④均方誤差增大不太多。

(2)方差膨脹因子法

方差膨脹因子cjj度量了多重共線性的嚴(yán)重程度，一般當(dāng)cjj>10時，模型就有嚴(yán)重的多重共線性，如果計(jì)算嶺估計(jì)βi(k)的協(xié)方差陣為：

上式中矩陣Cij(k)的對角元素cjj(k)就是嶺估計(jì)的方差膨脹因子，不難看出，cjj(k)隨著k的增大而減小。應(yīng)用方差膨脹因子選擇k的經(jīng)驗(yàn)做法是：選擇k使所有方差膨脹因子cjj(k)≤10。

此外還有Cp準(zhǔn)則法、Hoerl-Kennad公式法、Mcdorard-Gararneau法、雙h公式法等，因SAS的REG過程主要可以運(yùn)用嶺跡法及方差膨脹因子法，所以其他的方法在此不作介紹，有興趣的讀者可以參考相關(guān)的文獻(xiàn)。

第3步，根據(jù)嶺跡圖進(jìn)行變量篩選及重新確定k值。

嶺跡圖不僅能夠?qū)做出確定，還可以根據(jù)自變量的嶺跡曲線對自變量進(jìn)行篩選，也就是說可以根據(jù)自變量的嶺跡曲線來判斷該變量是否可以進(jìn)入回歸方程。把嶺跡應(yīng)用于回歸分析中自變量的選擇，其基本原則為：

(1)去掉嶺回歸系數(shù)比較穩(wěn)定且絕對值比較小的自變量。這里嶺回歸系數(shù)可以直接比較大小，因?yàn)樵O(shè)計(jì)陣X是假定已經(jīng)中心標(biāo)準(zhǔn)化了的。

(2)去掉嶺回歸系數(shù)不穩(wěn)定但隨著k值的增加迅速趨于零的自變量。

(3)去掉一個或若干個具有不穩(wěn)定嶺回歸系數(shù)的自變量。如果不穩(wěn)定的嶺回歸系數(shù)很多，究竟去掉幾個，去掉哪幾個，并無一般原則可遵循。這要結(jié)合已找出的復(fù)共線性關(guān)系以及去掉后重新進(jìn)行嶺回歸分析的效果來決定。

第4步，對模型進(jìn)行表達(dá)及作出專業(yè)結(jié)論。

在進(jìn)行嶺估計(jì)后，應(yīng)根據(jù)所估計(jì)的參數(shù)寫出回歸方程，并結(jié)合專業(yè)知識判斷方程中各自變量的系數(shù)及正負(fù)號是否符合實(shí)際情況。最后根據(jù)回歸系數(shù)的大小來判斷各自變量對因變量影響的大小及根據(jù)所求得的回歸方程進(jìn)行預(yù)測。

2　基于嶺回歸分析解決實(shí)際問題

2.1　問題與數(shù)據(jù)結(jié)構(gòu)

沿用文獻(xiàn)[2]中的“問題與數(shù)據(jù)”，并基于派生變量得到的“最優(yōu)回歸模型”所決定的“數(shù)據(jù)集”，來提出下面的“新問題”：即“weight”的回歸系數(shù)為“-88.00801”，這個“負(fù)值”表明體重越重的人收縮壓(SBP)越低，這似乎不符合臨床專業(yè)知識。盡管計(jì)算出來的因變量的預(yù)測值在專業(yè)上都成立，而且，模型的殘差方差=122.32418、R2=0.9931，這些結(jié)果都提示所構(gòu)建的多重線性回歸模型很好。但畢竟存在回歸系數(shù)的正負(fù)號不符合專業(yè)知識的“嚴(yán)重瑕疵”，這是一個需要徹底解決的“疑難問題”！

2.2　所需要的SAS程序

解決上述提及的“疑難問題”的一個有效方法就是使用“嶺回歸分析”。所需要的SAS程序如下：

data a1;

input id age height weight bmi sbp;

cards;

(此處輸入文獻(xiàn)[2]表1中50行6列數(shù)據(jù))

;

run;

/* 以上程序?yàn)榱藙?chuàng)建數(shù)據(jù)集a1 */

data a2;

set a1;

x1=age*age; x2=age*height; x3=age*weight;

x4=age*bmi; x5=height*height; x6=height*weight;

x7=height*bmi; x8=weight*weight; x9=weight*bmi;

x10=bmi*bmi;

run;

/* 以上程序是在數(shù)據(jù)集a1基礎(chǔ)上創(chuàng)建數(shù)據(jù)集a2，它增添了10個派生變量 */

proc reg data=a2;

model sbp=age height weight bmi x1-x10/noint selection=backward sls=0.05;

quit;

/* 以上程序是基于數(shù)據(jù)集a2擬合文獻(xiàn)[2]中那個“最佳”回歸模型 */

symbol1 v=x c=blue;

symbol2 v=circle c=black;

symbol3 v=square c=red;

symbol4 v=triangle c=green;

symbol5 v=dot c=yellow;

symbol6 v=# c=orange;

symbol7 v=% c=purple;

symbol8 v=$ c=blue;

legend1 mode=protect position=(bottom right inside)

across=3 cborder=black offset=(0,0) label=(color=blue position=(top center) 'independent variables') cframe=white;

/* 以上程序?yàn)榱嗽O(shè)置繪圖的基本條件 */

proc standard data=a2 m=0 s=1 out=a3;

run;

/* 以上程序?yàn)榱藢?shù)據(jù)集a2進(jìn)行標(biāo)準(zhǔn)化變換 */

/*因前面用后退法篩選自變量，故僅保留有統(tǒng)計(jì)學(xué)意義的自變量，列在下面的model語句*/

proc reg data=a3 outest=b1 ridge=0 to 0.1 by 0.01;

model sbp=age weight x3 x6-x10/noint;

plot /ridgeplot vref=0 lvref=1 nomodel legend=legend1 nostat;

quit;

/* 以上程序是基于標(biāo)準(zhǔn)化變換后的數(shù)據(jù)集a3進(jìn)行嶺回歸分析并繪制出嶺跡圖 */

proc print data=b1;

run;

/* 以上程序?yàn)榱溯敵鰩X回歸分析的計(jì)算結(jié)果 */

proc reg data=a2 outest=b2 ridge=0 to 0.1 by 0.01;

model sbp=age weight x3 x6-x10/noint;

quit;

/* 以上程序是基于數(shù)據(jù)集a2進(jìn)行嶺回歸分析以便獲得與原變量對應(yīng)的嶺回歸分析結(jié)果 */

/* 其中，“ridge=0 to 0.1 by 0.01”就是讓嶺參數(shù)k取一系列數(shù)值代入建模 */

proc print data=b2;

run;

/* 以上程序是輸出與原變量對應(yīng)的嶺回歸分析結(jié)果 */

【SAS程序說明】在以上的SAS程序中，都用“/* ……… */”注釋語句作了說明。

2.3　SAS輸出結(jié)果及其解釋

以上SAS程序輸出的結(jié)果非常多，因篇幅所限，此處從略。下面僅摘要給出最主要的結(jié)果。

由繪制出的嶺跡圖(圖略)可知：當(dāng)嶺參數(shù)k=0.01時，全部8個自變量的回歸系數(shù)就已經(jīng)趨向于穩(wěn)定狀態(tài)，但結(jié)合數(shù)據(jù)集b1的輸出結(jié)果可從數(shù)值清楚地看出：當(dāng)嶺參數(shù)k=0.08時，全部8個自變量的回歸系數(shù)的取值波動都在小數(shù)點(diǎn)后第2位上，因?yàn)閗的取值越大，自變量的回歸系數(shù)數(shù)值波動就會越小，但殘差方差會越大。故要求不高時，本例可取k=0.01；要求稍高一點(diǎn)時，可取k=0.08。

(1)標(biāo)準(zhǔn)化條件下嶺回歸分析所得到的全部自變量的回歸系數(shù)

由數(shù)據(jù)集b1可獲得標(biāo)準(zhǔn)化條件下嶺回歸分析所得到的全部自變量的回歸系數(shù)：

k=0.01和k=0.08時，全部自變量的標(biāo)準(zhǔn)化回歸系數(shù)：

kageweightx3x6x7x8x9x100.011.361760.370-1.248120.6060.3183-0.0352-0.1420.10000.080.577370.228-0.269360.2830.19050.0213-0.048-0.0693

(2)基于原變量的嶺回歸分析所得到的全部自變量的回歸系數(shù)

kageweightx3x6x7x8x9x100.011.400700.1400-0.0067670.002870.01370-0.000034-0.001100.005920.080.593880.0865-0.0014600.001340.008200.000021-0.00037-0.00410

(3)基于原變量的常規(guī)多重線性回歸分析所得到‘最佳結(jié)果’的全部自變量的回歸系數(shù)

ageweightx3x6x7x8x9x101.82182-88.00801-0.00970.645694.32456-0.058350.78530-2.62458

【說明】比較上面的“(2)”與“(3)”中各自變量前的回歸系數(shù)可知，引入派生變量并構(gòu)建出“最優(yōu)回歸模型”后，weight的回歸系數(shù)為“-88.00801”[見上面的“(3)”標(biāo)題之下第2行]，該負(fù)值不符合臨床專業(yè)知識；經(jīng)過嶺回歸分析后，得出：weight的回歸系數(shù)為“0.0865”[見上面的“(2)”標(biāo)題之下k=0.08的那一行]，這個系數(shù)是正值，與臨床專業(yè)知識相符。由此可知，在創(chuàng)建多重線性回歸模型的過程中，引入派生變量并采取多種方法進(jìn)行自變量篩選，在獲得“最佳多重線性回歸模型”后，再采用嶺回歸分析進(jìn)一步優(yōu)化多重線性回歸模型，這樣在獲得了較小殘差方差的多重線性回歸模型的基礎(chǔ)上，又能使少數(shù)自變量回歸系數(shù)的正負(fù)號不符合專業(yè)要求的問題得到合理解決。