谷恒明，胡良平,2(.軍事醫(yī)學科學院生物醫(yī)學統(tǒng)計學咨詢中心，北京 00850；2.世界中醫(yī)藥學會聯(lián)合會臨床科研統(tǒng)計學專業(yè)委員會，北京 00029 通信作者：胡良平，E-mail：lphu82@sina.com)

1 基于三種統(tǒng)計思想建立多重線性回歸模型的概述[1-3]

1.1 經(jīng)典多重線性回歸分析建模概述

多重線性回歸分析是用回歸方程定量地刻畫一個因變量與多個自變量之間的線性依存關系。其中，因變量是連續(xù)型變量，自變量是相互獨立的連續(xù)型變量(也常包括少量分類變量)。

經(jīng)典多重線性回歸分析的內(nèi)容包括對自變量的篩選和回歸診斷(含多重共線性診斷和異常點診斷)、對回歸模型和模型中全部參數(shù)的假設檢驗、對模型擬合效果的評價以及利用求得的回歸模型對因變量進行預測。

對自變量的篩選有八種方法：①前進法；②后退法；③逐步法；④最大R2增量法；⑤最小R2增量法；⑥R2選擇法；⑦修正R2選擇法；⑧Mallow’s Cp統(tǒng)計量選擇法。

當自變量之間存在較多的、嚴重的多重共線性關系時，通?？刹捎脦X回歸分析或者主成分回歸分析。

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

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

(1)

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

(2)

1.2 貝葉斯回歸分析建模概述

Yi=μi+εi,εi～N(0,σ2),i=1,2,…,n,

μi=β0+βiXi,

(3)

這里要為各個參數(shù)指定一個先驗分布，例如：

π(β0)=φ(0,var=le6)

π(β1)=φ(0,var=le6)

π(σ2)=fiΓ(shape=3/10,scale=10/3)

(4)

經(jīng)典多重線性回歸分析假定自變量前的回歸系數(shù)是固定的，而貝葉斯回歸分析認為參數(shù)是隨機的?；谪惾~斯統(tǒng)計思想建立回歸模型時，要為各個自變量前的參數(shù)(即回歸系數(shù))和殘差指定一個先驗分布。可以依靠經(jīng)驗或預試驗的結果指定各自合適的先驗分布；如果沒有辦法給定先驗分布，可以使用無信息先驗(相當于均勻分布)代替。貝葉斯回歸分析中沒有自變量篩選功能，因此要借助經(jīng)典多重線性回歸分析中篩選方法篩選出來的自變量來建立回歸模型。

貝葉斯統(tǒng)計建模可以參考SAS軟件的STAT模塊中MCMC過程來實現(xiàn)。MCMC方法即馬氏鏈蒙特卡羅方法，默認算法是使用正態(tài)分布隨機游動Metropolis算法。MCMC的抽樣方法有Gibbs抽樣、Metropolis抽樣、獨立性抽樣、隨機游動Metropolis抽樣等[2]。

1.3 機器學習回歸分析建模概述

有別于經(jīng)典統(tǒng)計思想和貝葉斯統(tǒng)計思想，機器學習統(tǒng)計思想則另辟蹊徑，它不依賴于概率分布知識、也不依賴于先驗分布知識，而是通過基于訓練樣本的學習獲取知識和經(jīng)驗，再用測試樣本來驗證。屬于機器學習的具體方法很多，通常包括決策樹法、支持向量機法、神經(jīng)網(wǎng)絡法、隨機森林法和集成學習法等[3]。

1.4 三類方法回歸建模效果

1.4.1 三類回歸建模效果的評價

情形一，樣本量較少時：分別使用兩種方法建立回歸模型，用相對誤差絕對值的均值(abserror)，殘差平方和(SSress)與決定系數(shù)(R2)作為評價指標。

情形二，樣本量較多時：①全部數(shù)據(jù)用來建立模型并比較，評價指標同樣本量較少時。②K-Fold交叉驗證，即全部數(shù)據(jù)拆分為K份，其中(K-1)份用作建立模型的訓練集，剩下一份當做測試集。訓練集擬合效果使用相對誤差絕對值的均值(abserror)，殘差均方(MSE)與決定系數(shù)(R2)作為評價指標；測試集使用相對誤差絕對值的均值(abserror)，殘差均方(MSE)與標準化均方誤差(NMSE)作為評價指標。③K-Fold交叉驗證中，K取值分別為10、7、4和2。當K取定一個數(shù)值后，分別重復抽取10次，即進行10次重復建模。

1.4.2 評價指標的具體公式

(5)

在數(shù)值上，NMSE等于1-R2，這里的R2是回歸的決定系數(shù)，但是，對于測試集來說，其NMSE與測試集回歸的R2沒有什么關系。交叉驗證主要關心測試集的NMSE。

(6)

N為樣本量，n為自變量個數(shù)。從上述的公式中可以看出，殘差均方MSE是殘差平方和與自由度的比值。交叉驗證中K取值不同，建立模型的訓練集和預測使用的測試集樣本量是不同的，直接基于殘差平方和比較不夠合理，因此，需要除以自由度。

說明：BP 神經(jīng)網(wǎng)絡建模效果的評價指標采用公式(5)。

2 實例及基于經(jīng)典統(tǒng)計思想的回歸分析

2.1 問題與數(shù)據(jù)

【例1】26例糖尿病患者的血清總膽固醇(X1)、甘油三酯(X2)、空腹胰島素(X3)、糖化血紅蛋白(X4)、空腹血糖(Y)的測量值列于表1，試基于經(jīng)典統(tǒng)計思想建立血糖與其他幾項指標間的多重線性回歸方程，并完成其他有關的任務。

2.2 回歸分析任務

對例1表1中的數(shù)據(jù)，設因變量為Y，自變量為X1、X2、X3、X4，試建立因變量依賴自變量的多重線性回歸模型，并做相應的假設檢驗。

表1 26例糖尿病患者血樣中有關指標的測定結果

注：詳細數(shù)據(jù)見本期第一篇文章《多重線性回歸分析的核心內(nèi)容與關鍵技術概述》

2.3 采用經(jīng)典統(tǒng)計思想實現(xiàn)多重線性回歸分析的方法

(1)不產(chǎn)生派生變量并采用三種篩選自變量的方法建模

data cra1;

input id x1-x4 y @@;

cards;

……

;

run;

proc reg data=cra1;

model y=x1-x4/selection=stepwise sle=0.5 sls=0.05;

run;

proc reg data=cra1;

model y=x1-x4/selection=forward sle=0.05;

run;

proc reg data=cra1;

model y=x1-x4/selection=backward sls=0.05;

run;

【程序說明】“cards”語句后的省略號代表表1中26行6列數(shù)據(jù)，其中，第1列為“編號”；REG過程被調(diào)用了3次，分別采用逐步法、前進法和后退法篩選自變量；“sle=0.5”代表選變量進入回歸模型的顯著性水平，其概率值選用0.5是非常大的，以便有較多的變量有機會進入回歸模型與其他變量進行組合，可以較好地保證單個作用不大但與某些自變量同時存在時作用明顯增大的自變量不會被排斥在回歸模型之外，這叫做“寬進”；“sls=0.05”代表已進入回歸模型的自變量仍能被保留在回歸模型之中的顯著性水平，其概率值選用0.05是統(tǒng)計學上被公認的顯著性水平，這叫做“嚴出”。

【主要輸出結果】本例資料采用上述三種篩選自變量方法所得結果完全相同，其主要結果如下：

方差分析來源自由度平方和均方F值Pr>F模型3156．1500252．0500117．77<．0001誤差2264．441132．92914校正合計25220．59115

變量參數(shù)估計值標準誤差II型SSF值Pr>FIntercept4．914802．1491915．318005．230．0322x20．437960．1342531．1716810．640．0036x3-0．299490．0969927．927119．530．0054x40．812670．2062445．4791815．530．0007

【輸出結果解釋】第1部分表明：總的回歸模型具有統(tǒng)計學意義(F=17.77、P<0.0001)；第2部分表明：自變量X1被淘汰掉了，其他3個自變量以及截距項均有統(tǒng)計學意義，得到的多重線性回歸方程為：

(2)不產(chǎn)生派生變量并采用逐步法篩選自變量且進行共線性診斷和殘差分析等方法建模

上面的SAS程序中的數(shù)據(jù)步程序不變，刪除第2個和第3個過程步程序，第1個過程步程序修改如下：

proc reg data=cra1;

model y= x1-x4/selection=stepwise sle=0.5

sls=0.05 collin collinoint vif tol r stb;

run;

【程序說明】“collin”選項是要求系統(tǒng)給出采用“方差比例”算法且未校正回歸模型中截距項影響的多重共線性診斷的結果；“collinoint”與前面選項的區(qū)別在于“校正了回歸模型中截距項影響”；“vif”選項是要求系統(tǒng)給出采用“方差膨脹因子”算法的多重共線性診斷的結果、“tol”等于“1/vif”，即要求系統(tǒng)給出采用“容許度”算法的多重共線性診斷的結果；“r”要求系統(tǒng)給出“殘差分析”的計算結果，有助于發(fā)現(xiàn)是否存在異常點；“stb”要求系統(tǒng)給出“標準化回歸系數(shù)”的計算結果。

【主要輸出結果及其解釋】逐步回歸分析的主要結果同上，此處從略?；?種方法進行共線性診斷的結果如下：

參數(shù)估計值變量自由度參數(shù)估計值標準誤差t值Pr>|t|標準化估計值容差方差膨脹Intercept14．914802．149192．290．03220．0x210．437960．134253．260．00360．382810．964301．03703x31-0．299490．09699-3．090．0054-0．374050．904871．10513x410．812670．206243．940．00070．485610．874261．14382

以上為第1部分輸出結果：倒數(shù)第3列為“標準化回歸系數(shù)”，其絕對值越大，表明所對應的自變量對因變量的貢獻就越大，由大到小依次為X4>X2>X3；倒數(shù)第2列和第1列分別為“容許度”與“方差膨脹因子”方法診斷共線性的結果，只需要看vif的數(shù)值是否大于10，大于10的那些行上的自變量間存在較嚴重的共線性。結果表明：3個自變量間不存在共線性關系。

共線性診斷個數(shù)特征值條件指數(shù)偏差比例Interceptx2x3x413．414191．000000．002000．025740．016570．0024120．377223．008490．001850．772910．169130．0005586330．194854．185920．017580．194120．598620．0436540．0137515．759980．978580．007230．215680．95338

以上為第2部分輸出結果：這是未對截距項進行校正且依據(jù)“方差比例”算法進行共線性診斷的結果。適用于“回歸分析模型中截距項無統(tǒng)計學意義”的場合，而本例截距項有統(tǒng)計學意義，故不需要看這部分輸出結果。

共線性診斷(截距已調(diào)整)個數(shù)特征值條件指數(shù)偏差比例x2x3x411．366571．000000．104390．241360．3119120．982611．179300．719610．245010．0005871330．650821．449060．176000．513620．68750

以上為第3部分輸出結果：這是對截距項進行校正后且依據(jù)“方差比例”算法進行共線性診斷的結果。適用于“回歸分析模型中截距項有統(tǒng)計學意義”的場合，而本例，截距項有統(tǒng)計學意義，故應該看這部分輸出結果。評判是否存在共線性的方法：看3個自變量列輸出結果的最后一行，當這些數(shù)值中有兩個或多個數(shù)值都很大且接近于1，那么，它們對應的自變量間存在共線性。本例3個自變量間不存在共線性關系。

殘差分析的輸出結果很多，此處從略。從學生化殘差結果來看，沒有取值的絕對值大于2的觀測點；從“Cook' s D”統(tǒng)計量計算結果來看，僅第25個觀測點的取值為0.698大于0.5，表明此觀測點是“可疑的異常點”。

(3)產(chǎn)生派生變量并采用三種篩選自變量的方法建模

所謂產(chǎn)生派生變量，就是除資料中已有的4個自變量外，再通過變量變換的方法，引入“新變量”。通?？梢砸胱宰兞康亩雾?，包括各自變量的平方項和任何兩個自變量的交叉乘積項。在上面第1段SAS程序基礎上，進行如下修改即可：

data cra2;

set cra1;

z1=x1*x1;z2=x2*x2;z3=x3*x3;z4=x4*x4;

z5=x1*x2;z6=x1*x3;z7=x1*x4;z8=x2*x3;

z9=x2*x4;z10=x3*x4;

run;

proc reg data=cra2;

model y=x1-x4 z1-z10/selection=stepwise sle=0.5sls=0.05;

run;

proc reg data=cra2;

model y=x1-x4 z1-z10/selection=forward sle=0.05;

run;

proc reg data=cra2;

model y=x1-x4 z1-z10/selection=backward sls=0.05;

run;

【程序說明】在已經(jīng)運行原先SAS程序的基礎上(即已經(jīng)創(chuàng)建了SAS數(shù)據(jù)集cra1)，再創(chuàng)建新數(shù)據(jù)集cra2，這就是前兩個語句的作用。新數(shù)據(jù)集中增加了z1-z10共10個新變量，它們是由原先的4個自變量產(chǎn)生的派生變量，代表10個二次項；3個過程步分別采用“逐步法”“前進法”和“后退法”篩選自變量，現(xiàn)在的自變量個數(shù)為14個。

以下是“逐步法”和“前進法”篩選自變量的結果：

方差分析來源自由度平方和均方F值Pr>F模型3163．1870154．3956720．85<0．0001誤差2257．404142．60928校正合計25220．59115

變量參數(shù)估計值標準誤差II型SSF值Pr>FIntercept6．633831．4014958．4617622．410．0001x30．728440．2845917．094866．550．0179z6-0．185080．0514233．7988912．950．0016z70．126960．02044100．6624738．58<0．0001

以下是“后退法”篩選自變量的結果：

方差分析來源自由度平方和均方F值Pr>F模型3172．9698357．6566126．64<0．0001誤差2247．621322．16461校正合計25220．59115

變量參數(shù)估計值標準誤差II型SSF值Pr>FIntercept8．669291．02421155．0855471．65<0．0001z40．045470．0087458．5510627．05<0．0001z50．051790．0117042．3827619．580．0002z6-0．054800．0156726．4665612．230．0020

考察以上兩個不同的建模結果可以發(fā)現(xiàn)：后退法的建模結果稍好一些，因為其總模型的假設檢驗結果的F=26.64較大(較小者為F=20.85)且模型中所含的項數(shù)相同(均為4項)，總模型和各項假設檢驗結果均具有統(tǒng)計學意義。

結合前面未引入派生變量時得到的多重線性回歸模型，其總模型的F=17.77，小于現(xiàn)在獲得的最好的模型對應的F=26.64，故擬合本例資料最好的多重線性回歸方程如下：

(4)產(chǎn)生派生變量并采用后退法篩選自變量且進行共線性診斷和殘差分析等方法建模

在運行了前面第1個數(shù)據(jù)步程序(創(chuàng)建數(shù)據(jù)集cra1)和第2個數(shù)據(jù)步程序(創(chuàng)建數(shù)據(jù)集cra2)的基礎上，再運行下面的過程步程序：

proc reg data=cra2;

model y= X1-X4 Z1-Z10/selection=backward

sls=0.05 collin collinoint vif tol r stb;

run;

輸出結果較多，此處從略。結果表明：Z4、Z5、Z6三個項之間不存在共線性關系；殘差分析的結果與前面陳述的結果基本相同，但第25個觀測點對應的Cook' s D=0.441<0.5，說明它已經(jīng)不屬于“可疑異常點”了。

3 小 結

本文介紹了基于經(jīng)典統(tǒng)計思想構建多重線性回歸模型的主要內(nèi)容，此法可用于很多多因素臨床試驗研究和觀察性研究中，其主要的應用條件是結果變量為“計量變量”，例如主要評價指標為“評分”。值得注意的是：許多實際工作者在對此類資料進行統(tǒng)計分析時，喜歡選用單因素差異性分析，見文獻[4-7]。其實，選用多重線性回歸分析方法來實現(xiàn)多因素分析，效果更佳。

雖然，在本文中所介紹的實例中，自變量都是計量變量，而在實際使用中，自變量可以是計量的、計數(shù)的、定性的。但是，當自變量為多值名義變量時，需要產(chǎn)生啞變量后才能引入多重線性回歸模型之中[8]，否則，可能得出錯誤的結論。

[1] Kleinbaum DG, Kupper LL, Muller KE, et al. Applied regression analysis and other multivariable methods[M]. 3版.北京: 機械工業(yè)出版社, 2003: 111-159.

[2] 劉金山, 夏強. 基于MCMC算法的貝葉斯統(tǒng)計方法[M]. 北京: 科學出版社, 2017: 118-174.

[3] 吳喜之. 復雜數(shù)據(jù)統(tǒng)計方法——基于R的應用[M]. 3版. 北京: 中國人民大學出版社, 2015: 41-56.

[4] 任傳波, 姜季妍, 董黎明. 青少年抑郁障礙患者心理社會學特征[J]. 四川精神衛(wèi)生, 2017, 30(5): 455-457.

[5] 徐華麗, 孫崇勇, 高悅. 大學生人格特質(zhì)對手機成癮傾向的影響[J]. 四川精神衛(wèi)生, 2017, 30(5): 458-462.

[6] 李青青, 張倩. 醫(yī)學院校大學生情緒狀態(tài)與學業(yè)拖延的關系[J]. 四川精神衛(wèi)生, 2017, 30(5): 463-465.

[7] 魏國英, 曾麗娟, 周桂成, 等. 精神科護士職業(yè)倦怠與工作壓力的相關性[J]. 四川精神衛(wèi)生, 2017, 30(5): 466-469.

[8] 胡良平. 醫(yī)學統(tǒng)計學——運用三型理論進行現(xiàn)代回歸分析[M]. 北京: 人民軍醫(yī)出版社, 2010: 9-33.