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

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

MCMC過程可應用于基于貝葉斯統(tǒng)計思想的生存資料的Cox比例風險回歸模型中，該方法的使用應基于數(shù)據(jù)滿足比例風險假定的前提上［1］。但在實際研究當中，往往會遇到生存資料不滿足該假定的情況，此時，可借助SAS軟件提供的多種過程步實現(xiàn)相關回歸模型的構建。本文將通過MCMC過程對效應隨時間變化的時間依賴型變量（即依時協(xié)變量）［2］構建基于貝葉斯統(tǒng)計思想的生存資料Cox非比例風險回歸模型，并對相關的主要內容加以說明。

1 含依時協(xié)變量的擴展Cox模型

經(jīng)典Cox比例風險回歸模型見式（1）：

βj是第j個協(xié)變量的回歸系數(shù)，h0(t)是基線風險函數(shù)。該模型有兩個重要特征：①基線風險函數(shù)依賴于t，但是與協(xié)變量無關且未知；②風險率依賴于協(xié)變量，但與時間t無關。然而，在實際的生存資料中，某些協(xié)變量Zj值會隨著時間而變化，即不滿足“比例風險假設”，如患者的心率。此時，我們可以在傳統(tǒng)的Cox比例風險回歸模型中加入該變量與時間的交互項，以描述其對基線風險函數(shù)的影響。帶有依時協(xié)變量的Cox回歸模型見式（2）：

需使用局部似然法來估計β，計算公式見式（3）：

Z（t）是時間與協(xié)變量的交互項。此時，關于HR（風險率）的統(tǒng)計推斷與經(jīng)典的Cox比例風險回歸模型類似，唯一的不同是風險率的主要部分exp[β'Z(t)]會隨著時間而變化，這就是“Cox非比例風險回歸模型”。

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

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

利用多發(fā)性骨髓瘤研究［3］的數(shù)據(jù)創(chuàng)建數(shù)據(jù)集Myeloma，所包含變量及解釋此處不再贅述，本文所涉及變量見表1。為了簡化計算，本例僅將表1中前三個定量自變量納入考慮，即將它們視為“依時協(xié)變量”。

表1 變量賦值表

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

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

SAS程序如下：

【程序說明】MODEL語句等號左邊定義生存時間和生存結局變量（括號內為截尾數(shù)據(jù)標識），右邊為各協(xié)變量（即自變量）以及時間與各協(xié)變量的交互項。BAYES語句指定回歸模型使用貝葉斯分析，并設定隨機數(shù)生成器種子seed=1；為了使初始值對后驗推斷的影響最小化，需要在Markov Chain達到目標分布后棄掉先前的部分樣本，因此，nmc用于設定棄掉先前的部分樣本后的迭代次數(shù)=10000。

【主要輸出結果及解釋】

這是輸出結果的“第1部分”，其中LogBUN、HGB、Platelet分別為樣本中的三個協(xié)變量，z2、z3、z4分別為各協(xié)變量與時間交互項。表中第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(3.2266×LogBUN-0.1395×z2)

圖1顯示LogBUN回歸系數(shù)的均值在3.0左右波動，隨著迭代次數(shù)的增加，擺動的幅度基本保持不變，所以有理由認為Markov Chain已經(jīng)收斂。圖2顯示變量LogBUN回歸系數(shù)不存在自相關，后驗樣本獨立。圖3顯示后驗密度核密度圖成鐘形，左右兩邊近似對稱。因篇幅所限，其他協(xié)變量的相關結果從略。

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

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

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

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

2.3.1 對觀測值進行轉置處理

【說明】因為Zj（ti）取決于ti，所以每個時期之和是當前時間ti和風險集中所有觀測結合的結果。風險集Ri包含所有生存時間大于等于ti的觀測。修改輸入的數(shù)據(jù)集，以使每一行不僅包含當前觀察值，而且還包含相應風險集中的所有觀察值。這樣，當為每個觀測值構建對數(shù)似然函數(shù)時，將得到所有相關數(shù)據(jù)。為此，首先需創(chuàng)建在不同行具有不同風險集的新數(shù)據(jù)集［4］，添加變量stop，此變量是指當前觀察值風險集中的觀察值數(shù)量。其余變量是整個數(shù)據(jù)集中模型協(xié)變量的轉置值。由于篇幅所限，SAS程序從略。

【主要輸出結果及解釋】

此處輸出結果（第3～6列）與前面輸出結果的第1部分（第3～6列）基本相同，各列的含義相同，此處從略。其中beta1～beta6分別對應各協(xié)變量及其與時間的交互項，即前面輸出結果中第1部分的“第1列”。因篇幅所限，馬爾可夫鏈迭代軌跡圖等從略。

2.3.2 對觀測值不進行轉置處理

對不獨立的數(shù)據(jù)建模，如果不想對每個觀測的數(shù)據(jù)進行轉置處理，還可以使用JOINTMODEL選項進行分析。因SAS程序過于復雜，此處從略。

3 討論與小結

3.1 討論

當生存資料經(jīng)比例風險假設檢驗結果判定為不滿足PH假定，即數(shù)據(jù)含有依時協(xié)變量時，研究者可以選用非比例風險回歸模型進行統(tǒng)計分析，進而得到協(xié)變量隨時間變化的趨勢性信息。

PHREG過程可以用來構建Cox非比例風險回歸模型，BAYES語句的應用則可實現(xiàn)回歸參數(shù)的貝葉斯估計。MCMC過程則是假設輸入的觀測值是獨立的，并且聯(lián)合對數(shù)似然函數(shù)是各個對數(shù)似然函數(shù)的總和。因為此過程統(tǒng)計推斷的原理是基于貝葉斯統(tǒng)計思想，所以在使用時要為數(shù)據(jù)指定似然函數(shù)，并為參數(shù)指定先驗分布。如果觀測值彼此不獨立，則該求和將產(chǎn)生錯誤的對數(shù)似然值。因此，對不滿足比例風險假定的數(shù)據(jù)建模，在MCMC過程中有兩種構建模型的方法：一是對觀測值進行轉置后，在MODEL語句中使用GENERAL函數(shù)，此函數(shù)可以運行不需要響應變量的蒙特卡洛模擬；二是不對觀測值進行轉置，則可以使用MCMC過程中的JOINTMODEL選項，其基本思想是將所有必需的數(shù)據(jù)集變量存儲在數(shù)組中，并且僅使用數(shù)組來構造整個數(shù)據(jù)集的對數(shù)似然函數(shù)。此時，MCMC過程將不再在模擬過程中逐步歷遍輸入數(shù)據(jù)，所以不再利用數(shù)據(jù)集變量構造對數(shù)似然函數(shù)，而是將數(shù)據(jù)集存儲在數(shù)組中，并用數(shù)組而不是數(shù)據(jù)集變量來計算對數(shù)似然值。

從結果來看，PHREG和MCMC過程中的三種建模方法均可實現(xiàn)基于貝葉斯統(tǒng)計思想構建生存資料Cox非比例風險回歸模型，并且所得結果差別不大，只是結果列出形式稍有不同；但從SAS程序上來看，使用MCMC過程所需要編寫的SAS程序冗長、復雜，而在PHREG過程中的“bayes語句”（其作用相當于一個很復雜的子程序）大大簡化了用戶的編程過程。

3.2 小結

對于不滿足比例風險假設的生存資料，可以選擇SAS軟件中的PHREG過程或MCMC過程構建基于貝葉斯統(tǒng)計思想的Cox非比例風險回歸模型，前者利用BAYES語句實現(xiàn)貝葉斯分析，后者則根據(jù)是否對觀測值進行轉置處理，提供兩種方法實現(xiàn)貝葉斯分析。三種方法建模所得結果基本一致，但基于PHREG過程來實現(xiàn)，SAS程序簡單得多。