穩(wěn)健回歸分析

胡良平

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

1　概　　述

1.1　何為穩(wěn)健回歸分析

設(shè)用Y代表因變量，X1、X2、…、Xm分別代表m個(gè)自變量，則多重線性回歸模型可以表示為：

Y=β0+β1X1+β2X2+…+βmXm+ε

(1)

式中β0為總體截距，β1、β2、…、βm分別為各個(gè)自變量所對(duì)應(yīng)的總體偏回歸系數(shù)，ε為隨機(jī)誤差，常假定其服從正態(tài)分布。偏回歸系數(shù)βi(i=1,2,…,m)表示在其他自變量固定不變的情況下，Xi每改變一個(gè)測量單位時(shí)所引起的因變量Y的平均改變量。多重線性回歸模型的樣本回歸方程可以表示為：

(2)

如何求出模型(1)中的參數(shù)(包括截距項(xiàng)和回歸系數(shù))呢？當(dāng)資料滿足一些前提條件(例如模型的誤差項(xiàng)服從正態(tài)分布、自變量互相獨(dú)立、不存在嚴(yán)重的異常點(diǎn))時(shí)，只需要采取普通的最小二乘法(簡稱OLS估計(jì)法，也叫做最小平方法)來構(gòu)造求解回歸系數(shù)的正規(guī)方程組，然后解此方程組，就可獲得全部參數(shù)的估計(jì)值。但是，當(dāng)資料中存在嚴(yán)重異常點(diǎn)時(shí)，就需要采用“穩(wěn)健回歸分析方法”來給出參數(shù)的估計(jì)值。

所謂穩(wěn)健回歸分析，就是在構(gòu)建多重線性回歸模型(1)時(shí)，不以普通的最小二乘法來構(gòu)造求解回歸系數(shù)的正規(guī)方程組，而是依據(jù)某些改進(jìn)的做法來構(gòu)造求解回歸系數(shù)的正規(guī)方程組。其目的就是依據(jù)推導(dǎo)出來的用于估計(jì)回歸參數(shù)的計(jì)算公式使參數(shù)的估計(jì)結(jié)果具有盡可能好的穩(wěn)定性，即盡可能降低或消除異常點(diǎn)對(duì)回歸分析結(jié)果的影響。

1.2　穩(wěn)健回歸分析應(yīng)用的場合

當(dāng)擬做多重線性回歸分析的原始數(shù)據(jù)存在較大比例的“異常點(diǎn)”且自變量間不存在嚴(yán)重多重共線性時(shí)，采用此方法構(gòu)建多重線性回歸模型，可以最大限度地消除“異常點(diǎn)”對(duì)建模結(jié)果造成的影響。

1.3　穩(wěn)健回歸分析的原理

此法的關(guān)鍵在于使估計(jì)出來的回歸系數(shù)比較穩(wěn)定，其實(shí)質(zhì)就是設(shè)法修改“普通最小二乘法”，使構(gòu)造出來的正規(guī)方程組對(duì)“異常點(diǎn)”不敏感，再通過類似于“迭代再加權(quán)最小二乘法”等方法求解正規(guī)方程組，從而獲得各回歸系數(shù)相對(duì)穩(wěn)定的估計(jì)值。具體的方法有多種，例如L估計(jì)、R估計(jì)、M估計(jì)、S估計(jì)和MM估計(jì)等[1-2]。其中有些估計(jì)方法還可以作進(jìn)一步細(xì)分，例如M估計(jì)可進(jìn)一步分為“Huber估計(jì)”“Tukey估計(jì)”和“中位數(shù)估計(jì)”。由于這些估計(jì)方法涉及很深的數(shù)學(xué)知識(shí)，在文獻(xiàn)[1]中用了8頁篇幅介紹了前述提及的估計(jì)方法，感興趣的讀者可參閱有關(guān)文獻(xiàn)，故此處從略。

2　基于穩(wěn)健回歸分析解決實(shí)際問題

2.1　此法可否解決多重共線性問題

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

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

2.1.2所需要的SAS程序

嘗試采用“穩(wěn)健回歸分析方法”解決上述提及的“疑難問題”。所需要的SAS程序如下：

data a1;

input id age height weight bmi sbp;

cards;

(此處輸入文獻(xiàn)[3]表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個(gè)派生變量 */

proc reg data=a2;

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

quit;

/* 以上程序是調(diào)用REG過程并基于數(shù)據(jù)集a2擬合文獻(xiàn)[3]中那個(gè)‘最佳’回歸模型 */

proc robustreg data=a2 method=m seed=100;

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

quit;

/* 以上程序是調(diào)用ROBUSTREG過程并基于數(shù)據(jù)集a2和M估計(jì)方法擬合文獻(xiàn)[3]中那個(gè)‘最佳’回歸模型 */

/*接下去，將上面SAS過程步中的關(guān)鍵詞“method=m”依次修改為：method=lts、method=s和method=mm，就是調(diào)用ROBUSTREG過程并基于數(shù)據(jù)集a2且分別用LTS估計(jì)法、S估計(jì)法和MM估計(jì)法來擬合文獻(xiàn)[3]中那個(gè)“最佳”回歸模型 */

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

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

以下將上面SAS程序中的5個(gè)過程步的輸出結(jié)果以濃縮的方式呈現(xiàn)出來，見表1。

【說明】比較上述由“REG過程”與基于“ROBUSTREG過程”并分別采用“M估計(jì)法”“LTS估計(jì)法”“S估計(jì)法”和“MM估計(jì)法”對(duì)同一個(gè)具有嚴(yán)重共線性問題的資料進(jìn)行多重線性回歸分析的結(jié)果可知：它們估計(jì)的“參數(shù)值”比較接近，但仍沒有消除多重共線性對(duì)回歸系數(shù)的嚴(yán)重影響(尤其是weight前的回歸系數(shù)為負(fù)值且其絕對(duì)值還比較大，不符合臨床專業(yè)知識(shí))。也就是說：SAS/STAT中的“ROBUSTREG過程”不能解決“多重共線性問題”。要想消除自變量間多重共線性的影響，常用的方法有兩種，第一種就是采用“主成分回歸分析”，見文獻(xiàn)[3]；第二種就是采用“嶺回歸分析”，參見本期中的《嶺回歸分析》一文。

2.2　此法可否解決資料中有較多異常點(diǎn)問題

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

【例2】假定有一個(gè)總樣本含量n=1000的數(shù)據(jù)集中包含10%異常點(diǎn)的資料，每組數(shù)據(jù)(即每個(gè)個(gè)體)包含三個(gè)變量(x1，x2，y)的觀測值，見表2。

表1　五種估計(jì)參數(shù)方法估計(jì)“例1資料”的主要結(jié)果

注：C與E分別代表“參數(shù)估計(jì)值”與“標(biāo)準(zhǔn)誤”；OLS、M、LTS、S、MM分別代表估計(jì)回歸模型中參數(shù)的方法依次為普通最小二乘法、M估計(jì)法、LTS估計(jì)法、S估計(jì)法和MM估計(jì)法；第7列上為空，因?yàn)長TS估計(jì)法給出的誤差項(xiàng)是“標(biāo)準(zhǔn)差”，與其他方法不一致(其他為“標(biāo)準(zhǔn)誤”)，故未呈現(xiàn)出來

表2　某資料中首尾各10組數(shù)據(jù)(n=1000)

注：此表省略編號(hào)為11到990之間的980行數(shù)據(jù)；在全部1000行數(shù)據(jù)中，最后100行數(shù)據(jù)為“異常點(diǎn)”，占10%

【特別說明】例2是人為構(gòu)造的，它來自SAS 9.3的ROBUSTREG過程中的“樣例”。三個(gè)變量“x1、x2和y”沒有實(shí)際的專業(yè)含義，僅為了造出一個(gè)樣本含量為1000且含10%異常點(diǎn)的數(shù)據(jù)集。

設(shè)定x1和x2及測量誤差e都是服從標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)變量(其均值為0、方差為1)，前900個(gè)y的數(shù)值按下面的模型(3)計(jì)算出來；最后100個(gè)y的數(shù)值按下面的模型(4)計(jì)算出來：

y=10+5*x1+3*x2+0.5*e

(3)

y=100+e

(4)

比較式(3)與式(4)可知：y的前900個(gè)數(shù)據(jù)中的每一個(gè)都在基數(shù)“10”基礎(chǔ)上再加上三項(xiàng)并不大的數(shù)值，其平均值大約為“10+5+3=18”；而y的后100個(gè)數(shù)據(jù)中的每一個(gè)都在基數(shù)“100”基礎(chǔ)上再加上一個(gè)隨機(jī)誤差，其平均值大約為100。由此可知：表2的1000行數(shù)據(jù)中，對(duì)因變量y而言，后100個(gè)y值明顯大于前900個(gè)y值，故屬于“異常值”，它們所對(duì)應(yīng)的那100行數(shù)據(jù)點(diǎn)就屬于“異常點(diǎn)”了。

【問題】試擬合表2中y依賴x1、x2變化的二重線性回歸模型。

2.2.2所需要的SAS程序

(1)先用下面的一段SAS數(shù)據(jù)步程序產(chǎn)生表2中的1000行3列數(shù)據(jù)，創(chuàng)建數(shù)據(jù)集aa。

data aa (drop=i);

do i=1 to 1000;

x1=rannor(1234);

x2=rannor(1234);

e=rannor(1234);

if i>900 then y=100 + e;

else y=10+5*x1+3*x2+0.5*e;

output;

end;

run;

/* 以上程序產(chǎn)生1000組數(shù)據(jù)(x1,x2,y)，其中，有10%的是異常值 */

(2)再用下面的SAS程序?qū)?shù)據(jù)集aa(即表2中的數(shù)據(jù))拷貝成數(shù)據(jù)集a，然后再調(diào)用SAS過程進(jìn)行建模(也可用前面例1中創(chuàng)建數(shù)據(jù)集a1的方法，此處從略)。

data a;

set aa;

proc reg data=a;

model y=x1 x2;

run;

proc robustreg data=a method=m seed=100;

model y=x1 x2;

run;

分別將上面的“method = m”后面的“m”替換成“l(fā)ts”“s”和“mm”，就可得到另三種回歸系數(shù)的穩(wěn)健估計(jì)結(jié)果。

【SAS過程步程序說明】第1個(gè)過程步調(diào)用REG過程創(chuàng)建二重線性回歸模型；從第2到第5個(gè)過程步都是調(diào)用ROBUSTREG過程構(gòu)建二重線性回歸模型，其區(qū)別在于它們分別采用M估計(jì)、LTS估計(jì)、S估計(jì)和MM估計(jì)方法來估計(jì)模型中的參數(shù)值。

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

以上SAS程序中5個(gè)SAS過程步輸出的主要結(jié)果列入表3中。

表3　五種估計(jì)參數(shù)方法估計(jì)“例2資料”的主要結(jié)果

注：C與E分別代表“參數(shù)估計(jì)值”與“標(biāo)準(zhǔn)誤”；OLS、M、LTS、S、MM分別代表估計(jì)回歸模型中參數(shù)的方法依次為“普通最小二乘法”、M估計(jì)法、LTS估計(jì)法、S估計(jì)法和MM估計(jì)法；第7列上為空，因?yàn)長TS估計(jì)法給出的誤差項(xiàng)是“標(biāo)準(zhǔn)差”，與其他方法不一致(其他為“標(biāo)準(zhǔn)誤”)，故未呈現(xiàn)出來

由表3可知：當(dāng)資料中存在10%異常點(diǎn)時(shí)，基于普通最小二乘法(OLS)給出的參數(shù)估計(jì)值偏離其真值(截距=10、x1前的斜率=5、x2前的斜率=3)很遠(yuǎn)，而基于“M估計(jì)法”“LTS估計(jì)法”“S估計(jì)法”和“MM估計(jì)法”給出的參數(shù)估計(jì)值都與其真值非常接近。再列出這5種算法對(duì)應(yīng)的復(fù)相關(guān)系數(shù)平方[即y與(x1和x2)的相關(guān)系數(shù)的平方]分別為：0.0234(OLS法)、0.7788(M估計(jì)法)、0.9932(LTS估計(jì)法)、0.9928(S估計(jì)法)和0.7520(MM估計(jì)法)。其中，OLS估計(jì)法的復(fù)相關(guān)系數(shù)平方非常小，而LTS估計(jì)法的復(fù)相關(guān)系數(shù)平方最大。

2.3　小結(jié)

用ROBUSTREG過程時(shí)應(yīng)選擇哪一種“穩(wěn)健估計(jì)”方法創(chuàng)建多重線性回歸模型很難一概而論。一般來說，在擬合優(yōu)度“Goodnes-of-Fit”評(píng)價(jià)的四個(gè)指標(biāo)中，R-Square取值越大、另三個(gè)評(píng)價(jià)指標(biāo)(因篇幅所限，在本文中省略了)取值越小，并且，回歸系數(shù)的“標(biāo)準(zhǔn)誤”越小越好，該估計(jì)方法的擬合效果就越好。就本例而言，總體來看，LTS估計(jì)法的效果最好。

最關(guān)鍵的問題在于：本例中的數(shù)據(jù)是基于上面的二重線性回歸模型(3)產(chǎn)生出來的。這個(gè)模型意味著：截距項(xiàng)為10、x1前的回歸系數(shù)為5、x2前的回歸系數(shù)為3，在此基礎(chǔ)上，加一個(gè)隨機(jī)誤差的二分之一。此處的“隨機(jī)誤差”是服從均值為0、方差為1的正態(tài)分布的隨機(jī)誤差?；诒?中的計(jì)算結(jié)果，可以寫出五個(gè)模型的具體表達(dá)式，見式(5)-式(9)。

【結(jié)論】以模型(3)為“金標(biāo)準(zhǔn)”，模型(5)偏離很遠(yuǎn)；模型(6)-(9)的質(zhì)量都很高。

若再結(jié)合“復(fù)相關(guān)系數(shù)平方”等評(píng)價(jià)指標(biāo)綜合評(píng)價(jià)，就本例而言，LTS估計(jì)法所得的結(jié)果最穩(wěn)健。