胡良平

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

1　概　　述

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

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

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

(1)

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

(2)

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

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

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

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

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

此法的關(guān)鍵在于使估計出來的回歸系數(shù)比較穩(wěn)定，其實(shí)質(zhì)就是設(shè)法修改“普通最小二乘法”，使構(gòu)造出來的正規(guī)方程組對“異常點(diǎn)”不敏感，再通過類似于“迭代再加權(quán)最小二乘法”等方法求解正規(guī)方程組，從而獲得各回歸系數(shù)相對穩(wěn)定的估計值。具體的方法有多種，例如L估計、R估計、M估計、S估計和MM估計等[1-2]。其中有些估計方法還可以作進(jìn)一步細(xì)分，例如M估計可進(jìn)一步分為“Huber估計”“Tukey估計”和“中位數(shù)估計”。由于這些估計方法涉及很深的數(shù)學(xué)知識，在文獻(xiàn)[1]中用了8頁篇幅介紹了前述提及的估計方法，感興趣的讀者可參閱有關(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”，這個“負(fù)值”表明：體重越重的人收縮壓(SBP)越低，這似乎不符合臨床專業(yè)知識。盡管計算出來的因變量的預(yù)測值在專業(yè)上都成立，且模型殘差的方差=122.32418、R2=0.9931，這些結(jié)果都提示所構(gòu)建的多重線性回歸模型很好。但畢竟存在回歸系數(shù)的正負(fù)號不符合專業(yè)知識的“嚴(yán)重瑕疵”，這是一個需要徹底解決的“疑難問題”！

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個派生變量 */

proc reg data=a2;

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

quit;

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

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

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

quit;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(3)

y=100+e

(4)

比較式(3)與式(4)可知：y的前900個數(shù)據(jù)中的每一個都在基數(shù)“10”基礎(chǔ)上再加上三項(xiàng)并不大的數(shù)值，其平均值大約為“10+5+3=18”；而y的后100個數(shù)據(jù)中的每一個都在基數(shù)“100”基礎(chǔ)上再加上一個隨機(jī)誤差，其平均值大約為100。由此可知：表2的1000行數(shù)據(jù)中，對因變量y而言，后100個y值明顯大于前900個y值，故屬于“異常值”，它們所對應(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)健估計結(jié)果。

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

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

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

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

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

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

2.3　小結(jié)

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

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

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

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