上篇文章介紹了多狀態(tài)模型，今天我們基于R軟件講解多狀態(tài)模型的具體實(shí)現(xiàn)。

1. 數(shù)據(jù)材料

本篇實(shí)例數(shù)據(jù)來(lái)自EBMT（歐洲血液和骨髓組織）提供的血癌患者骨髓移植后的生存數(shù)據(jù)，數(shù)據(jù)集包括1985年至1998年的2279例血癌患者，數(shù)據(jù)集中描述患者狀態(tài)的變量除了最常見(jiàn)的死亡（Death）和復(fù)發(fā)(Relapse)之外，還包括，骨髓移植后是否康復(fù)(Recovery, Rec)，骨髓移植后是否發(fā)生不良事件(Adverse event, AE)，以及骨髓移植后先康復(fù)后發(fā)生不良事件(AE and Rec)，共計(jì)五個(gè)結(jié)局狀態(tài)。同時(shí)數(shù)據(jù)集中還包括影響預(yù)后的四個(gè)協(xié)變量（分類(lèi)屬性變量），供體-受體是否性別匹配（match）,是否預(yù)防（proph）,骨髓移植的年份（year）,骨髓移植時(shí)患者的年齡（agecl），四個(gè)協(xié)變量的具體代表意義和簡(jiǎn)單統(tǒng)計(jì)描述見(jiàn)下表：

表1. 數(shù)據(jù)集中各協(xié)變量的信息表

本例中的多狀態(tài)模型關(guān)系圖：

模型中的六個(gè)狀態(tài)具體解釋如下：

（1）      骨髓移植后患者活著，病情還在緩解期，患者沒(méi)有康復(fù)或發(fā)生不良事件（Alive and in remission, no recovery or adverse event）；

（2）     患者存活于緩解期，從治療中康復(fù)（Alive in remission, recovered from the treatment）；

（3）      患者存活于緩解期，發(fā)生了不良事件（Alive in remission, occurrence of the adverse event）；

（4）      患者活著，康復(fù)和不良事件具有發(fā)生（Alive, both recovered and adverse event occurred）；

（5）      活著，但發(fā)生了復(fù)發(fā)事件（骨髓移植失?。?/span>Alive, in relapse (treatment failure)）；

（6）      死亡（治療失敗）（Dead (treatment failure)）。

R中數(shù)據(jù)的格式如下：

對(duì)于這類(lèi)數(shù)據(jù)，臨床研究中可以從各角度出發(fā)去展開(kāi)分析研究，例如：

（1）      年齡這一協(xié)變量是如何影響不同階段的患者疾病康復(fù)過(guò)程的？

（2）      骨髓移植的2個(gè)月后發(fā)生的的不良反應(yīng)是否改變了患者10年生存率？

（3）      對(duì)于一個(gè)康復(fù)后的病人，應(yīng)監(jiān)測(cè)哪些風(fēng)險(xiǎn)因素？

1. 實(shí)現(xiàn)方法

（1）      使用軟件

本篇中的數(shù)據(jù)案例采用R軟件中的mstate包實(shí)現(xiàn)。

（2）      數(shù)據(jù)集整理和統(tǒng)計(jì)描述

R語(yǔ)言分析之前，需要對(duì)原始數(shù)據(jù)重新整理，首先對(duì)生存時(shí)間和生存狀態(tài)的數(shù)據(jù)格式整理。

對(duì)圖1中的六個(gè)狀態(tài)的多狀態(tài)模型的描述，我們可以借用轉(zhuǎn)移矩陣來(lái)描述，我們這個(gè)模型中有六個(gè)狀態(tài)，可以建立一個(gè)六乘以六的矩陣（六行六列），矩陣中的每一個(gè)元素代表一條轉(zhuǎn)移路徑例如X（g,h）代表從g狀態(tài)轉(zhuǎn)移到h狀態(tài)的路徑，我們圖一中總共有12條路徑，如果X（g,h）=1，則代表從g狀態(tài)轉(zhuǎn)移到h狀態(tài)的路徑為第一條。

Mstate包中transMat( )函數(shù)可以實(shí)現(xiàn)上述轉(zhuǎn)移矩陣，具體代碼如下：

tmat <- transmat(x="list(c(2," 3,="" 5,="" 6),="" c(4,5,="" 6),="" c(4,="" 5,="" 6),="" c(5,="">

+ c(), c()), names = c('Tx','Rec', 'AE', 'Rec+AE', 'Rel','Death'))

結(jié)果：

上圖為轉(zhuǎn)移矩陣，第一行第二列為1，代表從TX狀態(tài)轉(zhuǎn)移至Rec狀態(tài)所走的路徑為第一條路徑，NA代表兩個(gè)狀態(tài)之間沒(méi)有轉(zhuǎn)移。

（3）      從R語(yǔ)言的原始數(shù)據(jù)可以看出，每一條觀測(cè)都有rec,ae,recae,rel等五條結(jié)局時(shí)間時(shí)間及狀態(tài)，患者的各種狀態(tài)路徑最多有12種，這在分析中有點(diǎn)混亂，因此R語(yǔ)言在處理前，對(duì)數(shù)據(jù)做了進(jìn)一步的轉(zhuǎn)換，具體將數(shù)據(jù)轉(zhuǎn)換成一行為一個(gè)狀態(tài)轉(zhuǎn)移路徑，具體格式如下：

上述格式即為，ID號(hào)為1的患者的所有狀態(tài)轉(zhuǎn)移情況。一號(hào)患者剛開(kāi)始狀態(tài)為1（Tx）可能轉(zhuǎn)移為狀態(tài)2（Rec）,3(AE),5(Rel),6(Death),患者在第22個(gè)月時(shí)從狀態(tài)1轉(zhuǎn)移到了狀態(tài)2，從第22個(gè)月開(kāi)始，患者可能轉(zhuǎn)移的狀態(tài)為4（Rec+AE）,5(Rel),6(Death),直到第995個(gè)月患者發(fā)生刪失，這期間患者都在狀態(tài)2未發(fā)生狀態(tài)轉(zhuǎn)移。數(shù)據(jù)格式中的id代表患者id，from 指的是起始狀態(tài)，to指的時(shí)轉(zhuǎn)移后的狀態(tài)，trans 指的是轉(zhuǎn)移矩陣中的第幾條路徑，Tstart和Tstop指的起始和終止時(shí)間，time指的是兩個(gè)時(shí)間的差值，status在當(dāng)前這個(gè)狀態(tài)中是否發(fā)生結(jié)局事件，后面的變量proph,year,agecl是指協(xié)變量。

（4）      數(shù)據(jù)集中從當(dāng)前狀態(tài)轉(zhuǎn)移到下一個(gè)狀態(tài)的人數(shù)頻率，占比，Mstate包中events( )函數(shù)都可給出：

結(jié)果中可以看到每個(gè)狀態(tài)轉(zhuǎn)移的人數(shù)，占比，發(fā)生結(jié)局事件數(shù)，占比情況等。

（5）      對(duì)于協(xié)變量proph,year,agecl，我們假設(shè)每個(gè)協(xié)變量對(duì)每個(gè)狀態(tài)的影響是不同的。下節(jié)分析時(shí)，需要不同狀態(tài)不同對(duì)待，因此對(duì)每個(gè)狀態(tài)需要指定相應(yīng)的協(xié)變量。這里我們通過(guò)設(shè)置啞變量（指示變量）的方式實(shí)現(xiàn)。例如對(duì)變量Z被分解到每個(gè)狀態(tài)中分析，在i到g的狀態(tài)轉(zhuǎn)移過(guò)程中Z_ig=Z，其他狀態(tài)則對(duì)應(yīng)為Z_ig=0.

例如對(duì)協(xié)變量year的具體設(shè)置如上面的形式。

其他協(xié)變量使用同樣的方式處理。至此數(shù)據(jù)整理部分完成。

1. 分析過(guò)程及結(jié)果解釋

（1）      第一步分析，先不考慮協(xié)變量的情況下，分析狀態(tài)之間的轉(zhuǎn)移強(qiáng)度和累計(jì)轉(zhuǎn)移風(fēng)險(xiǎn)。

考慮狀態(tài)g轉(zhuǎn)移到狀態(tài)h過(guò)程中，某時(shí)刻的瞬間轉(zhuǎn)移強(qiáng)度，公式為：            

這條公式假設(shè)多狀態(tài)模型是個(gè)馬氏過(guò)程，即當(dāng)前的轉(zhuǎn)移強(qiáng)度只與當(dāng)前狀態(tài)有關(guān)，與它前面的狀態(tài)沒(méi)有關(guān)系。

那么累計(jì)轉(zhuǎn)移風(fēng)險(xiǎn)即為對(duì)轉(zhuǎn)移強(qiáng)度求積分（類(lèi)似累加）：

這個(gè)累計(jì)轉(zhuǎn)移風(fēng)險(xiǎn)我們先用用survival包中的coxph( )函數(shù)實(shí)現(xiàn)分層分析:

這里，只要用到上面整理好的數(shù)據(jù)中的生存時(shí)間和生存狀態(tài)，對(duì)轉(zhuǎn)移狀態(tài)變量trans做分層分析即可，不加任何協(xié)變量。

在mstate中msift( )函數(shù)以coxph( )函數(shù)的輸出結(jié)果為輸入變量，進(jìn)行二次處理，可以計(jì)算出轉(zhuǎn)移風(fēng)險(xiǎn)和協(xié)方差。

上面的結(jié)果是msfit輸出的結(jié)果中的每個(gè)狀態(tài)上發(fā)生結(jié)局事件的時(shí)間點(diǎn)的累計(jì)轉(zhuǎn)移風(fēng)險(xiǎn)值。注意這里我們的時(shí)間time是用年做單位，對(duì)時(shí)間是先處理后分析，具體處理：

可以用plot( )函數(shù)將所有的轉(zhuǎn)移狀態(tài)的累計(jì)風(fēng)險(xiǎn)函數(shù)用不同顏色畫(huà)出來(lái)。

Probtrans（）函數(shù)可以預(yù)測(cè)出每個(gè)狀態(tài)之間的累計(jì)轉(zhuǎn)移概率

用plot()函數(shù)畫(huà)出堆積的轉(zhuǎn)移概率圖

兩條相鄰曲線之間的距離表示處于對(duì)應(yīng)轉(zhuǎn)移狀態(tài)的概率。這里的圖的狀態(tài)順序是可以人為設(shè)定的。上圖是以tx為基準(zhǔn)。

下面兩張圖可用AE的結(jié)果。

 （2）      下一步考慮將協(xié)變量納入分析,考慮模型：

此模型是狀態(tài)g到h的比例風(fēng)險(xiǎn)模型，Z_gh是狀態(tài)g到h的模型的協(xié)變量。

首先第一步需要檢驗(yàn)：是否協(xié)變量對(duì)所有狀態(tài)的影響是相同的？這個(gè)過(guò)程可以通過(guò)likelihood ratio test 實(shí)現(xiàn)，

我們先建立全模型，即考慮所有狀態(tài)，即狀態(tài)下的協(xié)變量的cox模型，使用survival 包里的coxph函數(shù)建模：

cfull<- coxph(surv(tstart,="" tstop,="" status)="" ~="" match.1="">

                 + match.2 + match.3 + match.4+ match.5 + match.6

                 + match.7 + match.8 + match.9+ match.10 + match.11

                 + match.12 + proph.1 + proph.2+ proph.3 + proph.4

                 + proph.5 + proph.6 + proph.7+ proph.8 + proph.9

                 + proph.10 + proph.11 +proph.12 + year1.1 + year1.2

                 + year1.3 + year1.4 + year1.5+ year1.6 + year1.7

               + year1.8 + year1.9 + year1.10 +year1.11 + year1.12

                 + year2.1 + year2.2 + year2.3 +year2.4 + year2.5

                 + year2.6 + year2.7 + year2.8+ year2.9 + year2.10

                 + year2.11 + year2.12 +agecl1.1 + agecl1.2 + agecl1.3

                 + agecl1.4 + agecl1.5 +agecl1.6 + agecl1.7 + agecl1.8

                 + agecl1.9 + agecl1.10 +agecl1.11 + agecl1.12 + agecl2.1

                 + agecl2.2 + agecl2.3 +agecl2.4 + agecl2.5 + agecl2.6

                 + agecl2.7 + agecl2.8 +agecl2.9 + agecl2.10 + agecl2.11

                 +agecl2.12 + strata(trans), data = msebmt, method = 'breslow')

下表是所有狀態(tài)下，各模型協(xié)變量的cox回歸系數(shù)，其中黑體加粗是指P值小于0.05.通過(guò)表格可以發(fā)現(xiàn)，同一個(gè)協(xié)變量對(duì)不同狀態(tài)下的系數(shù)是不同的。

表中數(shù)字是回歸系數(shù)，括號(hào)里數(shù)字是回歸系數(shù)的標(biāo)準(zhǔn)誤。對(duì)每一個(gè)協(xié)變量，估計(jì)的系數(shù)有正有負(fù)，代表狀態(tài)轉(zhuǎn)移時(shí)的風(fēng)險(xiǎn)高低。通過(guò)這張表，可以很好的看出，一個(gè)狀態(tài)向另外一個(gè)狀態(tài)轉(zhuǎn)移風(fēng)險(xiǎn)的高低情況。

（未完待續(xù)）

本公眾號(hào)精彩歷史文章：

04：如何在R軟件中求一致性指數(shù)( Harrell'concordance index:C-index)？

05：Nomogram 繪制原理及R&SAS實(shí)現(xiàn).

06  : Lasso方法簡(jiǎn)要介紹及其在回歸分析中的應(yīng)用

07  : 最優(yōu)模型選擇中的交叉驗(yàn)證（Cross validation）方法

08  : 用R語(yǔ)言進(jìn)行分位數(shù)回歸(Quantile Regression)

09  : 樣本數(shù)據(jù)中異常值（Outliers）檢測(cè)方法及SPSS & R實(shí)現(xiàn)

10  : 原始數(shù)據(jù)中幾類(lèi)缺失值（Missing Data）的SPSS及R處理方法

11  :  [Survival analysis] Kaplan-Meier法之SPSS實(shí)現(xiàn)

12  :  [Survival analysis] COX比例風(fēng)險(xiǎn)回歸模型在SPSS中的實(shí)現(xiàn)

13  :  用R繪制地圖：以疾病流行趨勢(shì)為例

14  :  數(shù)據(jù)挖掘方法：聚類(lèi)分析簡(jiǎn)要介紹 及SPSS&R實(shí)現(xiàn)

15  :  醫(yī)學(xué)研究中的Logistic回歸分析及R實(shí)現(xiàn)

16  :  常用的非參數(shù)檢驗(yàn)(Nonparametric Tests)總結(jié)

17  :  高中生都能看懂的最小二乘法原理

18  :  R語(yǔ)言中可實(shí)現(xiàn)的常用統(tǒng)計(jì)假設(shè)檢驗(yàn)總結(jié)（側(cè)重時(shí)間序列）

19  :  如何根據(jù)樣本例數(shù)、均數(shù)、標(biāo)準(zhǔn)差進(jìn)行T-Test和ANOVA

20  :  統(tǒng)計(jì)學(xué)中自由度的理解和應(yīng)用

21  :  ROC和AUC介紹以及如何計(jì)算AUC

22  :  支持向量機(jī)SVM介紹及R實(shí)現(xiàn)

23  :  SPSS如何做主成分分析?

24  : Bootstrap再抽樣方法簡(jiǎn)介

25  :  定量測(cè)量結(jié)果的一致性評(píng)價(jià)及 Bland-Altman 法的應(yīng)用 

26  :  使用R繪制熱圖及網(wǎng)絡(luò)圖  

27  :  幾種常用的雙坐標(biāo)軸圖形繪制 

28  :  遺失的藝術(shù)—諾謨圖（Nomogram） 

29  :  Nomogram 繪制原理及R&SAS實(shí)現(xiàn)（二） 

30  :  WOE:信用評(píng)分卡模型中的變量離散化方法 

31  :  結(jié)構(gòu)方程模型（SEM）簡(jiǎn)介及教程下載  

32  :  重復(fù)測(cè)量的多因素方差分析SPSS實(shí)現(xiàn)操作過(guò)程