劉媛媛 ，李長平 ，2，胡良平

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

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

1 基于“貝葉斯統(tǒng)計思想”構建Cox比例風險回歸模型

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

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

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

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

2 構建Cox比例風險回歸模型

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

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

表1 變量賦值表

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

2.2 采用“PHREG過程”且基于貝葉斯統(tǒng)計思想構建Cox比例風險回歸模型

可以利用PHREG過程的BAYES語句擬合Cox比例風險回歸模型。

SAS程序如下：

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

【主要輸出結果及解釋】

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

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

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

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

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

2.3 采用“MCMC過程”且基于貝葉斯統(tǒng)計思想構建Cox比例風險回歸模型

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

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

SAS程序如下：

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

【主要輸出結果及解釋】

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

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

2.3.2 利用JOINTMODEL選項擬合模型

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

3 討論與小結

3.1 討論

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

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

3.2 小結

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