劉媛媛 ，李長(zhǎng)平 ，2，胡良平

（1.天津醫(yī)科大學(xué)公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計(jì)學(xué)教研室，天津 300070；2.世界中醫(yī)藥學(xué)會(huì)聯(lián)合會(huì)臨床科研統(tǒng)計(jì)學(xué)專業(yè)委員會(huì)，北京 100029；3.軍事科學(xué)院研究生院，北京 100850*通信作者：胡良平，E-mail：lphu927@163.com）

頻率學(xué)派和貝葉斯學(xué)派是統(tǒng)計(jì)學(xué)發(fā)展歷史上兩個(gè)重要的學(xué)派［1］。通常，前者認(rèn)為隨機(jī)事件的概率是客觀存在并假設(shè)固定不變的；而后者則認(rèn)為此“概率”是隨機(jī)的而不是固定不變的，并服從于某種分布。也即，經(jīng)典統(tǒng)計(jì)學(xué)認(rèn)為“參數(shù)”是固定不變的“常量”，而貝葉斯統(tǒng)計(jì)學(xué)認(rèn)為“參數(shù)”是隨機(jī)變量，這是兩者的根本分歧所在。1763年由Richard Price整理發(fā)表了貝葉斯的成果《An Essay towards solving a Problem in the Doctrine of Chances》提出了貝葉斯公式，并介紹了貝葉斯思想，其核心內(nèi)容就是對(duì)參數(shù)的估計(jì)并不是單純?nèi)Q于客觀數(shù)據(jù)，而是取決于客觀數(shù)據(jù)（包括總體和樣本信息）和先驗(yàn)信息的共同作用［2］。隨著計(jì)算機(jī)技術(shù)的發(fā)展和貝葉斯方法的改進(jìn)，特別是馬爾科夫蒙特卡洛（Markov chain Monte Carlo，MCMC）方法的提出和應(yīng)用，使得參數(shù)后驗(yàn)分布的模擬得以更方便地實(shí)現(xiàn)，從而體現(xiàn)出該法在處理小樣本數(shù)據(jù)時(shí)的優(yōu)勢(shì)［3］?，F(xiàn)在，越來越多的新的統(tǒng)計(jì)分析方法將經(jīng)典統(tǒng)計(jì)分析方法和貝葉斯思想有機(jī)地結(jié)合起來，例如，基于貝葉斯理論和生存分析相結(jié)合的貝葉斯生存分析在近年來越來越多地被應(yīng)用于不同的研究領(lǐng)域，尤其是醫(yī)學(xué)科學(xué)研究中［4-5］。因此，本文將介紹基于PHREG過程和MCMC過程分別構(gòu)建貝葉斯統(tǒng)計(jì)思想框架下生存資料的Cox比例風(fēng)險(xiǎn)回歸模型的相關(guān)內(nèi)容。

1 基于“貝葉斯統(tǒng)計(jì)思想”構(gòu)建Cox比例風(fēng)險(xiǎn)回歸模型

姚婷婷等［6］的文章已經(jīng)介紹了Cox比例風(fēng)險(xiǎn)回歸模型，見式（1）：

式（1）中，X1、X2、…、Xp為與生存時(shí)間可能有關(guān)的自變量（即影響因素）；h(t)為具有自變量X1、X2、…、Xp的個(gè)體在t時(shí)刻的風(fēng)險(xiǎn)函數(shù)；h0(t)為所有自變量為0時(shí)t時(shí)刻的基準(zhǔn)風(fēng)險(xiǎn)函數(shù)；β1、β2、…、βp為各自變量的偏回歸系數(shù)。偏回歸系數(shù)的估計(jì)需借助偏似然函數(shù)，用最大似然估計(jì)方法得到。偏似然函數(shù)的計(jì)算公式見式（2）：

式（2）中，qi為第i個(gè)死亡時(shí)點(diǎn)的條件死亡概率，其分子部分為第i個(gè)個(gè)體在ti（t1≤t2≤…≤ti≤…≤tk）死亡時(shí)點(diǎn)的風(fēng)險(xiǎn)函數(shù)h(ti)，分母部分為生存時(shí)間T≥ti的所有個(gè)體（包括死亡和刪失）的風(fēng)險(xiǎn)函數(shù)之和。

基于貝葉斯統(tǒng)計(jì)思想構(gòu)建生存資料回歸模型，即在原模型的基礎(chǔ)上利用貝葉斯方法的基本原理對(duì)回歸參數(shù)進(jìn)行估計(jì)的過程。所以，需要先對(duì)這些參數(shù)指定適當(dāng)?shù)南闰?yàn)分布，如果先驗(yàn)分布選擇不合適，則會(huì)對(duì)結(jié)果產(chǎn)生影響。故文獻(xiàn)［7］建議將Cox回歸模型系數(shù)βi的先驗(yàn)分布設(shè)定為正態(tài)分布，本研究也將按此進(jìn)行先驗(yàn)分布的設(shè)置。

2 構(gòu)建Cox比例風(fēng)險(xiǎn)回歸模型

2.1 實(shí)例與數(shù)據(jù)

本文所用數(shù)據(jù)來自一項(xiàng)骨髓瘤研究，研究者用烷基化劑治療65例患者，在隨訪期間內(nèi)，死亡48例，存活17例。變量賦值情況見表1。

表1 變量賦值表

【說明】完整數(shù)據(jù)來自文獻(xiàn)［8］。因篇幅所限，此處未呈現(xiàn)全部數(shù)據(jù)。

2.2 采用“PHREG過程”且基于貝葉斯統(tǒng)計(jì)思想構(gòu)建Cox比例風(fēng)險(xiǎn)回歸模型

可以利用PHREG過程的BAYES語句擬合Cox比例風(fēng)險(xiǎn)回歸模型。

SAS程序如下：

【程序說明】MODEL語句是PHREG過程的必需語句，等號(hào)左邊定義生存時(shí)間和生存結(jié)局變量（括號(hào)內(nèi)為截尾數(shù)據(jù)標(biāo)識(shí)），右邊為各協(xié)變量（即自變量）。使用BAYES語句則要求回歸模型的貝葉斯分析是采用Gibbs抽樣，同時(shí)設(shè)定seed為隨機(jī)數(shù)生成器種子，設(shè)為1；NMC為退火（退火是指為了使初始值對(duì)后驗(yàn)推斷的影響最小化，需要在Markov Chain達(dá)到目標(biāo)分布后棄掉先前的部分樣本）后的迭代次數(shù)，設(shè)為10000；OUTPOST選項(xiàng)將后驗(yàn)分布樣本保存在SAS數(shù)據(jù)集中以進(jìn)行以后的處理。

【主要輸出結(jié)果及解釋】

這是輸出結(jié)果的“第1部分”，第1列“參數(shù)”實(shí)際上是擬創(chuàng)建的回歸模型中的“自變量”；第2列指隨機(jī)重復(fù)抽樣一萬次；第3列“均值”實(shí)際上是各自變量前的回歸系數(shù)的估計(jì)值，而且，其中每個(gè)估計(jì)值都是一萬次隨機(jī)重復(fù)抽樣計(jì)算所得到結(jié)果的算術(shù)平均值；第4列為與“各均值”對(duì)應(yīng)的“標(biāo)準(zhǔn)差”；最后兩列為與“各均值”對(duì)應(yīng)的95%HPD（highest posterior density，HPD）區(qū)間，即95%最高后驗(yàn)密度置信區(qū)間。因此，根據(jù)此置信區(qū)間是否包含“0”（包含0時(shí)表明該變量對(duì)結(jié)果的影響無統(tǒng)計(jì)學(xué)意義），可得以下回歸方程：

h(t)=h0(t)exp(1.7610×LogBUN)

圖1 變量LogBUN回歸參數(shù)的馬爾可夫鏈迭代軌跡圖

圖2 變量LogBUN回歸參數(shù)的自相關(guān)函數(shù)圖

圖1顯示馬爾可夫鏈已經(jīng)收斂（因篇幅所限，其他協(xié)變量的相關(guān)結(jié)果此處從略）。

2.3 采用“MCMC過程”且基于貝葉斯統(tǒng)計(jì)思想構(gòu)建Cox比例風(fēng)險(xiǎn)回歸模型

圖3 變量LogBUN回歸參數(shù)的后驗(yàn)密度核密度圖

2.3.1 利用LAG函數(shù)擬合模型

SAS程序如下：

【程序說明】ARRAY語句用于將回歸系數(shù)的名稱與協(xié)變量、常量相關(guān)聯(lián)。PARMS語句給出模型中的參數(shù)名稱，并為其指定初始值。PRIOR語句指定模型參數(shù)的先驗(yàn)分布為正態(tài)分布。程序中的bZ為似然函數(shù)中的回歸項(xiàng)即每個(gè)觀測(cè)的風(fēng)險(xiǎn)集項(xiàng)。本例所用似然函數(shù)為Breslow似然函數(shù)。符號(hào)“l(fā)”為每個(gè)觀測(cè)計(jì)算對(duì)數(shù)似然的公式。IF-ELSE語句將所有的觀測(cè)分成三部分，并使用lag1函數(shù)來檢驗(yàn)兩個(gè)相鄰的生存時(shí)間time是否不同。符號(hào)“v”為生存狀態(tài)（Vstatus）的總和，因?yàn)閯h失數(shù)據(jù)不進(jìn)入似然公式的計(jì)算，所以需要將其去掉。MODEL語句用于在給定參數(shù)的似然函數(shù)的情況下指定數(shù)據(jù)的條件分布。

【主要輸出結(jié)果及解釋】

這里的輸出結(jié)果（特指第3列到第6列）與前面輸出結(jié)果的第1部分（特指第3列到第6列）是基本相同的，各列的含義相同，此處從略。

其中beta1～beta9分別對(duì)應(yīng)各協(xié)變量，實(shí)際上就是前面輸出結(jié)果中第1部分的“第1列”。

2.3.2 利用JOINTMODEL選項(xiàng)擬合模型

若利用PROC語句中的JOINTMODEL選項(xiàng)，則可使對(duì)數(shù)似然函數(shù)適用于整個(gè)數(shù)據(jù)集，而不只是單個(gè)的觀察值。但在使用此選項(xiàng)前，還要為數(shù)據(jù)風(fēng)險(xiǎn)集S指定包含其中的觀察的數(shù)量，為此需先創(chuàng)建一個(gè)stop變量。因篇幅所限，這部分的SAS程序和輸出結(jié)果從略。

3 討論與小結(jié)

3.1 討論

以Cox比例風(fēng)險(xiǎn)回歸模型為例，對(duì)于滿足PH假定的有刪失數(shù)據(jù)的生存資料來說，該模型能夠很好地利用偏似然（Partial Likelihood，PL）估計(jì)理論識(shí)別響應(yīng)變量的影響因素，但其應(yīng)用仍需一定的樣本量。貝葉斯思想的提出則有效地彌補(bǔ)了小樣本的不足。二者相結(jié)合的基于貝葉斯統(tǒng)計(jì)思想的Cox比例風(fēng)險(xiǎn)回歸模型很好地融合了其各自優(yōu)勢(shì)。

本研究通過實(shí)例，基于貝葉斯統(tǒng)計(jì)思想并借助PHREG過程和MCMC過程分別構(gòu)建了Cox比例風(fēng)險(xiǎn)回歸模型，并且分別利用MCMC過程中的LAG函數(shù)、JOINTMODEL選項(xiàng)擬合模型。通過對(duì)程序及結(jié)果的解釋比較，不難發(fā)現(xiàn)此三種方法均可得到后驗(yàn)樣本統(tǒng)計(jì)描述指標(biāo)（包括均值、標(biāo)準(zhǔn)差及95%最大后驗(yàn)密度置信區(qū)間）、有效樣本大小、馬爾可夫鏈迭代軌跡圖等主要結(jié)果，而且后驗(yàn)樣本的均值等指標(biāo)在數(shù)值上相差不大，MCMC過程中兩種方法的結(jié)果更是完全一致，只是在結(jié)果展示形式上稍有不同。由于PHREG過程本身是實(shí)現(xiàn)經(jīng)典統(tǒng)計(jì)學(xué)中Cox比例風(fēng)險(xiǎn)模型回歸分析的標(biāo)準(zhǔn)過程，這里只需加入BAYES語句用于指定進(jìn)行貝葉斯估計(jì)。因此，通過比較不同過程的具體語句可以發(fā)現(xiàn)，PHREG過程相對(duì)MCMC過程要簡(jiǎn)潔，并且可以提供普通的Cox比例風(fēng)險(xiǎn)模型的結(jié)果——回歸系數(shù)的最大似然估計(jì)的結(jié)果（包括回歸系數(shù)估計(jì)值、標(biāo)準(zhǔn)誤差及其95%置信區(qū)間），而MCMC過程只是提供了“平均意義”下的類似結(jié)果。

3.2 小結(jié)

在SAS中，可以借助PHREG過程或MCMC過程構(gòu)建基于貝葉斯統(tǒng)計(jì)思想的Cox比例風(fēng)險(xiǎn)回歸模型，兩種做法的主要結(jié)果相似，但前者程序語句相對(duì)簡(jiǎn)潔。研究者可根據(jù)具體情況選擇其一進(jìn)行回歸模型的構(gòu)建。