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

1 MCMC方法概述

1.1 基于貝葉斯統(tǒng)計思想建立回歸模型的基本思路

若單純基于貝葉斯統(tǒng)計思想來建立多重線性回歸模型，需要進行“共軛先驗下的貝葉斯推斷”和(或)“廣義先驗下的貝葉斯推斷”求出多重線性回歸模型中的“估計回歸系數(shù)矩陣”和“估計誤差協(xié)方差矩陣”[1]。然而，隨著統(tǒng)計學的發(fā)展，統(tǒng)計學家又提出了新的思路：基于隨機抽樣進行統(tǒng)計模擬并結合貝葉斯統(tǒng)計思想進行計算，能夠得到基于數(shù)萬次、甚至數(shù)十萬次抽樣樣本條件下更加穩(wěn)定的回歸系數(shù)估計結果，誤差很小，結果的精確度很高。這個方法叫做馬爾科夫鏈蒙特卡羅方法，簡稱為MCMC(Markov chain Monte Carlo)方法[2-3]。

1.2 馬爾科夫鏈

馬爾科夫鏈是一種具有馬爾科夫性的特別隨機過程。所謂隨機過程，用通俗的語言來表述，就是以“時間”為自變量的一個變量(即時間的函數(shù))。例如我國自2000年以來各年某種疾病的發(fā)病率資料，每個發(fā)病率數(shù)據(jù)都與特定的年份有關。這樣的一串數(shù)被稱為“時間序列”，也被稱為一個隨機過程。

那么，什么樣的隨機過程可以被稱為馬爾科夫鏈呢？用通俗的語言來表述，在隨機過程或時間序列中，任何一個時間點的取值僅受其前一個時間點上取值的影響，而不受其他時間點上的數(shù)值影響，這一特性被稱為“馬爾科夫性”，并且，人們稱具有馬爾科夫性的隨機過程或時間序列為“馬爾科夫鏈”。

1.3 蒙特卡羅方法的基本思路與實質

蒙特卡羅是摩納哥國的世界聞名賭城，其實，蒙特卡羅方法與賭城“蒙特卡羅”之間沒有任何關聯(lián)！事實上，美國在研究原子彈的過程中，遇到了復雜的有關核反應的計算問題，波蘭裔美國數(shù)學家斯塔尼斯拉夫·吳拉姆(Stanislaw Ulam)創(chuàng)造性地提出了隨機模擬方法，用來計算當時遇到的復雜計算問題，并由科學家馮·紐曼(Von Neumann)在計算機上實現(xiàn)。為了保密，他們就將該方法稱為“蒙特卡羅(Monte Carlo)”法[3-4]。

蒙特卡羅方法的基本思想是模擬從總體中反復抽取一批又一批樣本，然后利用所抽取的多批模擬樣本進行估計、假設檢驗等統(tǒng)計推斷。所以，蒙特卡羅方法的實質就是基于有限的樣本資料，通過反復抽取隨機樣本，使“原先的樣本”接近其所對應的“總體”情況。換一句話說，蒙特卡羅方法的實質是基于“有限樣本所提供的信息”去洞察“總體的真實情況”。

1.4 馬爾科夫鏈蒙特卡羅方法

前面介紹蒙特卡羅方法時提到了一個關鍵詞，即“隨機抽樣”。隨機抽樣不等于“胡亂”或“隨便”操作，而是依據(jù)一定的“規(guī)則”實施抽樣，這里所說的“規(guī)則”在數(shù)學上就叫做“抽樣算法”。

非常著名的抽樣算法有：梅切波利斯(Metropolis,1953)算法、由哈斯廷斯(Hastings,1970)對梅切波利斯算法改進的算法(即MH算法)、吉曼(Geman,1984)兄弟在研究數(shù)字圖像恢復問題時提出的吉布斯抽樣法(Gibbs sampling)以及其他一些模擬退火方法(simulated annealing)[5]，現(xiàn)在人們將所有這些方法和各種各樣的推廣統(tǒng)稱為馬爾科夫鏈蒙特卡羅(MCMC)方法。

2 實例及基于貝葉斯統(tǒng)計思想實現(xiàn)多重線性回歸分析

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

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

2.2 回歸分析任務

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

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

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

2.3 采用經(jīng)典統(tǒng)計思想篩選自變量

實際資料中往往包含了很多自變量，而其中有些自變量對因變量的影響作用無統(tǒng)計學意義，故需要對自變量進行篩選；在SAS中，并非所有用于回歸分析的SAS過程都具有篩選自變量的功能。在處理多重線性問題時，一般都用REG過程來實現(xiàn)對自變量的篩選。這部分內容參見本期第二篇文章《基于經(jīng)典統(tǒng)計思想實現(xiàn)多重線性回歸分析》，對例1資料篩選的結果為：保留X2、X3、X4三個自變量。

2.4 采用貝葉斯回歸分析方法建?！付ㄏ闰灧植紴镹～(0，106)

(1)基于自變量篩選結果并借助SAS中MCMC過程創(chuàng)建回歸模型所需要的SAS程序

proc mcmc data=cra1 seed =12345 nbi=10000

nmc=100000 thin=10 diag=all OUTPOST=POST STATS=ALL MCHISTORY=DETAILED dic;

parms beta0 0 beta1 100 beta2 0 beta3 10 ;

parms sigma2 20;

prior beta0 beta1 beta2 beta3～ normal(mean=0, var=1e6);

prior sigma2 ～ igamma(shape=0.001, scale=0.001);

mu = beta0+beta1*x2+beta2*x3+beta3*x4;

run;

【說明】MCMC過程中的語句和選項的含義比較復雜，因篇幅所限，此處從略。

(2)多重線性回歸分析的輸出結果

PosteriorSummariesParameterNMeanStandardDeviationPercentiles25%50%75%beta0100004．92742．23983．44564．90536．3989beta1100000．43890．14170．34540．43880．5333beta210000-0．30070．1009-0．3666-0．3000-0．2346beta3100000．81200．21430．67440．81030．9547sigma2100003．22861．07662．48053．02663．7333

多重線性回歸方程為：

相對誤差絕對值的平均值為0.1149813505，殘差平方和64.4418，決定系數(shù)R2=0.70787。

橋梁改建于2009年10月開工，2012年7月建成通車。2015年8月發(fā)現(xiàn)橋梁第六聯(lián)(第19～22跨)T梁整體有較明顯扭轉，病害的主要表現(xiàn)為：

2.5 采用貝葉斯回歸分析方法建?！獰o信息先驗分布

(1)基于自變量篩選結果并借助SAS中MCMC過程創(chuàng)建回歸模型所需要的SAS程序

proc mcmc data =cra1 seed =12345 nbi=10000

nmc=100000 thin=10 diag=all OUTPOST =POST STATS = ALL MCHISTORY=DETAILED dic;

parms beta0 0 beta1 100 beta2 0 beta3 10;

parms sigma2 20;

prior beta0 beta1 beta2 beta3～ general(0);

priorsigma2=-log(sigma2);

prior sigma2 ～ general(priorsigma2);

mu = beta0+ beta1 * x2 + beta2 * x3+beta3*x4;

model y ～ normal(mu,var = sigma2);

run;

(2)多重線性回歸分析的輸出結果

PosteriorSummariesParameterNMeanStandardDeviationPercentiles25%50%75%beta0100004．92222．24983．43164．90646．3549beta1100000．43610．14060．34430．43590．5278beta210000-0．30150．1012-0．3688-0．3009-0．2349beta3100000．81350．21500．67490．81650．9525sigma2100003．19651．06892．45372．98723．7111

多重線性回歸方程為：

相對誤差絕對值的平均值為0.11495609，殘差平方和64.4434，決定系數(shù)R2=0.70786，標準化均方誤差NMSE=1-R2=0.292140。

基于經(jīng)典統(tǒng)計思想和貝葉斯統(tǒng)計思想(給定兩種先驗分布)建立多重線性回歸分析模型，對模型擬合效果進行評價的結果見表2。

表2 基于經(jīng)典統(tǒng)計與貝葉斯統(tǒng)計建模模型評價結果的比較

注：先驗分布指正態(tài)分布、無信息先驗實際上是取均勻分布作為先驗分布；NMSE為標準化均方誤差

3 小 結

數(shù)據(jù)本身質量(即自變量間不存在多重共線性關系、無嚴重的異常點)較好時，即經(jīng)典統(tǒng)計思想建模方法所要求的假定條件基本上能夠得到滿足，此時，經(jīng)典統(tǒng)計建模的效果就比較好；如果有足夠準確的先驗信息，采用經(jīng)典統(tǒng)計思想進行多重線性回歸分析和采用貝葉斯統(tǒng)計思想回歸分析的結果相差不大。

【說明】感興趣的讀者可以參照本期第二篇文章《基于經(jīng)典統(tǒng)計思想實現(xiàn)多重線性回歸分析》中的做法，產(chǎn)生派生變量，并基于自變量篩選得到Z4、Z5、Z6三項，再將這三項分別代入本文中兩段調用MCMC過程的SAS程序中去計算，會得到相應的多重線性回歸分析結果。因篇幅所限，此處從略。

值得一提的是：文獻[6-9]中的資料都可以采用本文介紹的方法進行多重線性回歸分析，其結果和結論更具有可信性。

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

[2] 茆詩松, 王靜龍, 濮曉龍. 高等數(shù)理統(tǒng)計[M]. 北京: 高等教育出版社, 2010: 440-459.

[3] 黃長全. 貝葉斯統(tǒng)計及其R實現(xiàn)[M]. 北京: 清華大學出版社, 2017: 114-138.

[4] 高惠璇. 統(tǒng)計計算[M]. 北京: 北京大學出版社, 2017: 173-223.

[5] 康崇祿. 蒙特卡羅方法理論和應用[M]. 北京: 科學出版社, 2017: 86-440.

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

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

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

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