|
上篇文章介绍了多状态模型，今天我们基于R软件讲解多状态模型的具体实现。
本篇实例数据来自EBMT（欧洲血液和骨髓组织）提供的血癌患者骨髓移植后的生存数据，数据集包括1985年至1998年的2279例血癌患者，数据集中描述患者状态的变量除了最常见的死亡（Death）和复发(Relapse)之外，还包括，骨髓移植后是否康复(Recovery, Rec)，骨髓移植后是否发生不良事件(Adverse event, AE)，以及骨髓移植后先康复后发生不良事件(AE and Rec)，共计五个结局状态。同时数据集中还包括影响预后的四个协变量（分类属性变量），供体-受体是否性别匹配（match）,是否预防（proph）,骨髓移植的年份（year）,骨髓移植时患者的年龄（agecl），四个协变量的具体代表意义和简单统计描述见下表： 表1. 数据集中各协变量的信息表
本例中的多状态模型关系图： 模型中的六个状态具体解释如下： （1） 骨髓移植后患者活着，病情还在缓解期，患者没有康复或发生不良事件（Alive and in remission, no recovery or adverse event）； （2） 患者存活于缓解期，从治疗中康复（Alive in remission, recovered from the treatment）； （3） 患者存活于缓解期，发生了不良事件（Alive in remission, occurrence of the adverse event）； （4） 患者活着，康复和不良事件具有发生（Alive, both recovered and adverse event occurred）； （5） 活着，但发生了复发事件（骨髓移植失败）（Alive, in relapse (treatment failure)）； （6） 死亡（治疗失败）（Dead (treatment failure)）。 R中数据的格式如下： 对于这类数据，临床研究中可以从各角度出发去展开分析研究，例如： （1） 年龄这一协变量是如何影响不同阶段的患者疾病康复过程的？ （2） 骨髓移植的2个月后发生的的不良反应是否改变了患者10年生存率？ （3） 对于一个康复后的病人，应监测哪些风险因素？
（1） 使用软件 本篇中的数据案例采用R软件中的mstate包实现。 （2） 数据集整理和统计描述 R语言分析之前，需要对原始数据重新整理，首先对生存时间和生存状态的数据格式整理。 对图1中的六个状态的多状态模型的描述，我们可以借用转移矩阵来描述，我们这个模型中有六个状态，可以建立一个六乘以六的矩阵（六行六列），矩阵中的每一个元素代表一条转移路径例如X（g,h）代表从g状态转移到h状态的路径，我们图一中总共有12条路径，如果X（g,h）=1，则代表从g状态转移到h状态的路径为第一条。 Mstate包中transMat( )函数可以实现上述转移矩阵，具体代码如下： 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')) 结果： 上图为转移矩阵，第一行第二列为1，代表从TX状态转移至Rec状态所走的路径为第一条路径，NA代表两个状态之间没有转移。 （3） 从R语言的原始数据可以看出，每一条观测都有rec,ae,recae,rel等五条结局时间时间及状态，患者的各种状态路径最多有12种，这在分析中有点混乱，因此R语言在处理前，对数据做了进一步的转换，具体将数据转换成一行为一个状态转移路径，具体格式如下： 上述格式即为，ID号为1的患者的所有状态转移情况。一号患者刚开始状态为1（Tx）可能转移为状态2（Rec）,3(AE),5(Rel),6(Death),患者在第22个月时从状态1转移到了状态2，从第22个月开始，患者可能转移的状态为4（Rec+AE）,5(Rel),6(Death),直到第995个月患者发生删失，这期间患者都在状态2未发生状态转移。数据格式中的id代表患者id，from 指的是起始状态，to指的时转移后的状态，trans 指的是转移矩阵中的第几条路径，Tstart和Tstop指的起始和终止时间，time指的是两个时间的差值，status在当前这个状态中是否发生结局事件，后面的变量proph,year,agecl是指协变量。 （4） 数据集中从当前状态转移到下一个状态的人数频率，占比，Mstate包中events( )函数都可给出： 结果中可以看到每个状态转移的人数，占比，发生结局事件数，占比情况等。 （5） 对于协变量proph,year,agecl，我们假设每个协变量对每个状态的影响是不同的。下节分析时，需要不同状态不同对待，因此对每个状态需要指定相应的协变量。这里我们通过设置哑变量（指示变量）的方式实现。例如对变量Z被分解到每个状态中分析，在i到g的状态转移过程中Zig=Z，其他状态则对应为Zig=0. 例如对协变量year的具体设置如上面的形式。 其他协变量使用同样的方式处理。至此数据整理部分完成。
（1） 第一步分析，先不考虑协变量的情况下，分析状态之间的转移强度和累计转移风险。 考虑状态g转移到状态h过程中，某时刻的瞬间转移强度，公式为：
这条公式假设多状态模型是个马氏过程，即当前的转移强度只与当前状态有关，与它前面的状态没有关系。 那么累计转移风险即为对转移强度求积分（类似累加）：
这个累计转移风险我们先用用survival包中的coxph( )函数实现分层分析:
这里，只要用到上面整理好的数据中的生存时间和生存状态，对转移状态变量trans做分层分析即可，不加任何协变量。 在mstate中msift( )函数以coxph( )函数的输出结果为输入变量，进行二次处理，可以计算出转移风险和协方差。
上面的结果是msfit输出的结果中的每个状态上发生结局事件的时间点的累计转移风险值。注意这里我们的时间time是用年做单位，对时间是先处理后分析，具体处理：
可以用plot( )函数将所有的转移状态的累计风险函数用不同颜色画出来。 Probtrans（）函数可以预测出每个状态之间的累计转移概率
用plot()函数画出堆积的转移概率图 两条相邻曲线之间的距离表示处于对应转移状态的概率。这里的图的状态顺序是可以人为设定的。上图是以tx为基准。 下面两张图可用AE的结果。
（2） 下一步考虑将协变量纳入分析,考虑模型：
此模型是状态g到h的比例风险模型，Zgh是状态g到h的模型的协变量。 首先第一步需要检验：是否协变量对所有状态的影响是相同的？这个过程可以通过likelihood ratio test 实现， 我们先建立全模型，即考虑所有状态，即状态下的协变量的cox模型，使用survival 包里的coxph函数建模： 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') 下表是所有状态下，各模型协变量的cox回归系数，其中黑体加粗是指P值小于0.05.通过表格可以发现，同一个协变量对不同状态下的系数是不同的。
表中数字是回归系数，括号里数字是回归系数的标准误。对每一个协变量，估计的系数有正有负，代表状态转移时的风险高低。通过这张表，可以很好的看出，一个状态向另外一个状态转移风险的高低情况。 （未完待续） 本公众号精彩历史文章： 04：如何在R软件中求一致性指数( Harrell'concordance index:C-index)？ 05：Nomogram 绘制原理及R&SAS实现. 06 : Lasso方法简要介绍及其在回归分析中的应用 07 : 最优模型选择中的交叉验证（Cross validation）方法 08 : 用R语言进行分位数回归(Quantile Regression) 09 : 样本数据中异常值（Outliers）检测方法及SPSS & R实现 10 : 原始数据中几类缺失值（Missing Data）的SPSS及R处理方法 11 : [Survival analysis] Kaplan-Meier法之SPSS实现 12 : [Survival analysis] COX比例风险回归模型在SPSS中的实现 13 : 用R绘制地图：以疾病流行趋势为例 14 : 数据挖掘方法：聚类分析简要介绍 及SPSS&R实现 15 : 医学研究中的Logistic回归分析及R实现 16 : 常用的非参数检验(Nonparametric Tests)总结 17 : 高中生都能看懂的最小二乘法原理 18 : R语言中可实现的常用统计假设检验总结（侧重时间序列） 19 : 如何根据样本例数、均数、标准差进行T-Test和ANOVA 20 : 统计学中自由度的理解和应用 21 : ROC和AUC介绍以及如何计算AUC 22 : 支持向量机SVM介绍及R实现 23 : SPSS如何做主成分分析? 24 : Bootstrap再抽样方法简介 25 : 定量测量结果的一致性评价及 Bland-Altman 法的应用 26 : 使用R绘制热图及网络图 27 : 几种常用的双坐标轴图形绘制 28 : 遗失的艺术—诺谟图（Nomogram） 29 : Nomogram 绘制原理及R&SAS实现（二） 30 : WOE:信用评分卡模型中的变量离散化方法 32 : 重复测量的多因素方差分析SPSS实现操作过程 |
|
|
