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

1 MCMC方法概述

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

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

1.2 馬爾科夫鏈

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

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

1.3 蒙特卡羅方法的基本思路與實(shí)質(zhì)

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

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

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

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

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

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

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

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

2.2 回歸分析任務(wù)

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

表1 26例糖尿病患者血樣中有關(guān)指標(biāo)的測定結(jié)果

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

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

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

2.4 采用貝葉斯回歸分析方法建?！付ㄏ闰?yàn)分布為N～(0，106)

(1)基于自變量篩選結(jié)果并借助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過程中的語句和選項(xiàng)的含義比較復(fù)雜，因篇幅所限，此處從略。

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

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梁整體有較明顯扭轉(zhuǎn)，病害的主要表現(xiàn)為：

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

(1)基于自變量篩選結(jié)果并借助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)多重線性回歸分析的輸出結(jié)果

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，標(biāo)準(zhǔn)化均方誤差NMSE=1-R2=0.292140。

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

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

注：先驗(yàn)分布指正態(tài)分布、無信息先驗(yàn)實(shí)際上是取均勻分布作為先驗(yàn)分布；NMSE為標(biāo)準(zhǔn)化均方誤差

3 小 結(jié)

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

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

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

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

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

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

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

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

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

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

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

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