Xianxiong Ma
2020年3月16日
某研究人员收集了本市2007年确诊为轻度认知损害(MCI)的518例老年患者资料，包括基本人口学特征、生活方式、体格检查和合并疾病信息等，并于2010~2013年完成6次随访调查，主要观察结局为发生阿尔兹海默病(AD)。随访期间，共发生AD 78例，失访84例，其中28例搬迁、31例退出、25例死亡。试问影响MCI向AD转归的因素都有哪些?
生存分析是预后研究中比较常见的统计分析方法，但是经典的生存分析-般只关心一个终点事件(即研究者感兴趣的结局)，而医学研究中观察的终点往往并不唯一(即出现不感兴趣的结局)。
比如MCI患者在观察期间死于癌症、心血管疾病、车祸等原因而未发生AD，就不能为AD的发病做出贡献，即死亡“竞争”了AD的发生。传统统计方法将发生AD前死亡的个体、失访个体和未发生AD个体均按删失数据(censored data) 处理，可能会导致估计偏差。
上述例子应该选用竞争风险模型。竞争风险模型(Competing Risk Model)是一种处理多种潜在结局生存数据的分析方法，早在1999年Fine和Gray就提出了部分分布的半参数比例风险模型，通常使用的终点指标是累积发生率函数(Cumulative incidence function, CIF)
本例中可以将发生AD前死亡作为AD的竞争风险事件，采用竞争风险模型进行分析。竞争风险的单因素分析常用来估计关心终点事件的发生率，多因素分析常用来探索预后影响因素及效应值。
对于死亡率较高的老年人群，当有竞争风险事件存在时，采用传统生存分析方法(K-M法、Cox比例风险回归模型)会高估所研究疾病的发生风险，产生竞争风险偏倚，有人专门]研究发现约46%的文献可能存在这种偏倚。
一般的统计软件没有进行竞争风险分析的相应模块，可以用R软件的“cmprsk” 程序包进行编程。同时，SAS 9.4版本也新增了竞争风险分析模块，可以用PHREG语句完成，主要是在原Model选项中增加了“eventcode=“ 语句，用来指明哪个取值为感兴趣的结局，同时还增加了直接绘图功能。
研究骨髓移植对比血液移植治疗白血病的疗效，结局定义为“复发”，假定患者移植后不幸因为移植不良反应死亡，那这些发生移植相关死亡的患者就无法观察到“复发”的终点，也就是说“移植相关死亡”与“复发”存在竞争风险。
Sex : 性别； M=男，F=女
D : 疾病； ALL，AML
Phase : 疾病所处阶段； CR1, CR2, CR3, Relapse
Source : 移植类型； BM+PB，PB
Age : 年龄； 年
Ftime : 失败时间（事件发生时长）； 月
Status : 结局状态； 0=删失，1=复发，2=竞争风险事件
library(foreign)
bmt <- read.csv("bmtcrr.csv")
head(bmt, 10)
## Sex D Phase Age Status Source ftime
## 1 M ALL Relapse 48 2 BM+PB 0.67
## 2 F AML CR2 23 1 BM+PB 9.50
## 3 M ALL CR3 7 0 BM+PB 131.77
## 4 F ALL CR2 26 2 BM+PB 24.03
## 5 F ALL CR2 36 2 BM+PB 1.47
## 6 M ALL Relapse 17 2 BM+PB 2.23
## 7 M ALL CR1 7 0 BM+PB 124.83
## 8 F ALL CR1 17 2 BM+PB 10.53
## 9 M ALL CR1 26 0 BM+PB 123.90
## 10 F ALL Relapse 8 1 BM+PB 2.00
bmt$D <- as.factor(bmt$D)
library(survival)
library(cmprsk)
library(splines)
attach(bmt)
crmod <- cuminc(ftime, Status, D) #构建单因素生存函数
crmod
## Tests:
## stat pv df
## 1 2.8623325 0.09067592 1
## 2 0.4481279 0.50322531 1
## Estimates and Variances:
## $est
## 20 40 60 80 100 120
## ALL 1 0.3713851 0.3875571 0.3875571 0.3875571 0.3875571 0.3875571
## AML 1 0.2414530 0.2663827 0.2810390 0.2810390 0.2810390 NA
## ALL 2 0.3698630 0.3860350 0.3860350 0.3860350 0.3860350 0.3860350
## AML 2 0.4439103 0.4551473 0.4551473 0.4551473 0.4551473 NA
##
## $var
## 20 40 60 80 100 120
## ALL 1 0.003307032 0.003405375 0.003405375 0.003405375 0.003405375 0.003405375
## AML 1 0.001801156 0.001995487 0.002130835 0.002130835 0.002130835 NA
## ALL 2 0.003268852 0.003373130 0.003373130 0.003373130 0.003373130 0.003373130
## AML 2 0.002430406 0.002460425 0.002460425 0.002460425 0.002460425 NA
# 上面1代表复发（1）相对于删失，2代表竞争风险事件（2，移植？）相对于删失
plot(crmod, xlab = "Month", ylab = "CIF", lwd = 2, lty = 1, col = c("red", "blue",
"black", "forestgreen"))
# 定义所有可能为竞争风险因素的协方差矩阵
cov1 <- data.frame(age = bmt$Age,
sex_F = ifelse(bmt$Sex=='F',1,0),
dis_AML = ifelse(bmt$D=='AML',1,0),
phase_cr1 = ifelse(bmt$Phase=='CR1',1,0),
phase_cr2 = ifelse(bmt$Phase=='CR2',1,0),
phase_cr3 = ifelse(bmt$Phase=='CR3',1,0), #没设置relapse就是以relapse为参照
source_PB = ifelse(bmt$Source=='PB',1,0)) ## 手动设置哑变量
cov1
## age sex_F dis_AML phase_cr1 phase_cr2 phase_cr3 source_PB
## 1 48 0 0 0 0 0 0
## 2 23 1 1 0 1 0 0
## 3 7 0 0 0 0 1 0
## 4 26 1 0 0 1 0 0
## 5 36 1 0 0 1 0 0
## 6 17 0 0 0 0 0 0
## 7 7 0 0 1 0 0 0
## 8 17 1 0 1 0 0 0
## 9 26 0 0 1 0 0 0
## 10 8 1 0 0 0 0 0
## 11 29 0 0 0 0 0 0
## 12 13 0 1 0 0 0 0
## 13 29 1 1 0 1 0 0
## 14 7 0 0 0 0 0 0
## 15 40 0 0 0 0 0 0
## 16 20 1 0 0 1 0 0
## 17 18 0 0 0 0 1 0
## 18 29 0 1 1 0 0 0
## 19 6 0 0 1 0 0 0
## 20 33 0 1 0 1 0 0
## 21 27 1 0 0 0 1 0
## 22 22 1 0 0 1 0 1
## 23 15 0 0 0 0 1 1
## 24 8 0 0 0 1 0 1
## 25 46 0 1 0 1 0 1
## 26 21 0 0 0 0 0 1
## 27 21 1 0 0 0 0 1
## 28 7 0 1 0 0 0 1
## 29 4 0 0 0 0 0 1
## 30 53 1 1 0 0 0 1
## 31 27 0 1 0 1 0 1
## 32 28 0 1 0 0 0 1
## 33 23 1 0 0 1 0 1
## 34 17 0 0 0 0 0 1
## 35 19 1 0 0 0 0 1
## 36 39 0 1 0 0 0 1
## 37 15 0 0 0 1 0 1
## 38 47 0 1 0 1 0 1
## 39 16 1 1 0 0 0 1
## 40 17 1 0 1 0 0 1
## 41 26 0 1 0 0 1 1
## 42 44 0 1 0 1 0 1
## 43 30 1 1 0 0 0 1
## 44 32 0 1 0 0 0 1
## 45 20 0 1 1 0 0 1
## 46 38 0 1 1 0 0 1
## 47 17 0 0 0 0 0 1
## 48 38 0 1 0 0 0 1
## 49 36 0 1 0 1 0 1
## 50 34 0 1 0 1 0 1
## 51 42 1 1 0 1 0 1
## 52 58 1 1 0 0 0 1
## 53 30 1 1 0 0 0 1
## 54 37 1 0 1 0 0 1
## 55 25 0 0 1 0 0 1
## 56 47 1 1 1 0 0 1
## 57 46 1 0 1 0 0 1
## 58 21 0 1 0 1 0 1
## 59 17 1 1 0 0 0 1
## 60 25 0 0 0 0 0 1
## 61 45 1 1 0 1 0 1
## 62 17 0 0 1 0 0 1
## 63 44 0 1 0 0 0 1
## 64 20 1 1 1 0 0 1
## 65 18 0 1 0 1 0 1
## 66 46 1 1 0 1 0 1
## 67 49 1 1 0 0 0 1
## 68 17 0 0 0 0 0 1
## 69 36 1 0 0 0 0 1
## 70 14 0 0 0 1 0 1
## 71 22 0 0 0 0 0 1
## 72 38 0 1 0 0 1 1
## 73 20 1 0 0 1 0 1
## 74 26 1 0 0 0 0 1
## 75 20 0 0 0 0 0 1
## 76 28 0 1 1 0 0 1
## 77 45 1 1 0 0 0 1
## 78 17 0 1 0 0 0 1
## 79 40 0 1 0 1 0 1
## 80 9 1 0 0 0 0 1
## 81 18 0 1 0 0 0 1
## 82 33 1 1 0 0 0 1
## 83 33 0 0 1 0 0 1
## 84 17 1 1 0 1 0 1
## 85 17 1 0 0 0 0 1
## 86 28 0 1 0 1 0 1
## 87 26 1 1 1 0 0 1
## 88 32 0 0 0 1 0 1
## 89 62 0 1 0 0 0 1
## 90 46 0 1 0 0 0 1
## 91 20 1 1 1 0 0 1
## 92 50 1 1 0 1 0 1
## 93 40 0 1 0 0 0 1
## 94 43 0 1 1 0 0 1
## 95 31 1 1 0 0 1 1
## 96 46 1 0 1 0 0 1
## 97 51 0 1 0 0 0 1
## 98 39 0 0 0 0 0 1
## 99 19 1 0 0 0 0 1
## 100 33 1 0 1 0 0 1
## 101 23 1 1 0 0 0 1
## 102 50 1 0 1 0 0 1
## 103 36 1 1 1 0 0 1
## 104 36 0 1 0 1 0 1
## 105 35 0 0 1 0 0 1
## 106 41 1 1 0 0 0 1
## 107 28 1 1 1 0 0 1
## 108 30 0 1 1 0 0 1
## 109 39 1 1 1 0 0 1
## 110 52 0 1 1 0 0 1
## 111 38 0 1 1 0 0 1
## 112 35 0 1 1 0 0 1
## 113 27 1 1 0 0 1 1
## 114 38 0 1 1 0 0 1
## 115 22 1 1 0 0 0 1
## 116 26 0 0 1 0 0 1
## 117 47 1 1 0 1 0 1
## 118 17 0 1 0 0 0 1
## 119 21 0 0 0 0 0 1
## 120 28 0 0 1 0 0 1
## 121 24 0 0 0 1 0 1
## 122 35 1 1 1 0 0 1
## 123 19 1 0 0 0 1 1
## 124 23 0 0 1 0 0 1
## 125 33 0 0 0 1 0 1
## 126 38 1 0 0 0 1 1
## 127 36 0 1 1 0 0 1
## 128 48 0 1 0 0 0 1
## 129 25 1 1 0 0 0 1
## 130 22 0 0 0 1 0 1
## 131 37 1 1 0 0 0 1
## 132 52 0 1 0 1 0 1
## 133 53 1 1 1 0 0 1
## 134 22 0 1 0 1 0 1
## 135 20 0 0 0 1 0 1
## 136 51 0 1 1 0 0 1
## 137 33 0 1 0 0 0 1
## 138 33 1 1 0 0 1 1
## 139 20 0 0 1 0 0 1
## 140 37 0 1 0 0 0 1
## 141 14 1 1 1 0 0 1
## 142 33 0 1 1 0 0 1
## [ reached 'max' / getOption("max.print") -- omitted 35 rows ]
mod1 <- crr(bmt$ftime, bmt$Status, cov1, failcode=1, cencode=0)
# 构建多因素竞争风险模型，cov1为可能影响结局的协方差矩阵，failcode为感兴趣的结局事件，cencode为删失，其余为竞争风险因素
summary(mod1)
## Competing Risks Regression
##
## Call:
## crr(ftime = bmt$ftime, fstatus = bmt$Status, cov1 = cov1, failcode = 1,
## cencode = 0)
##
## coef exp(coef) se(coef) z p-value
## age -0.0185 0.982 0.0119 -1.554 0.1200
## sex_F -0.0352 0.965 0.2900 -0.122 0.9000
## dis_AML -0.4723 0.624 0.3054 -1.547 0.1200
## phase_cr1 -1.1018 0.332 0.3764 -2.927 0.0034
## phase_cr2 -1.0200 0.361 0.3558 -2.867 0.0041
## phase_cr3 -0.7314 0.481 0.5766 -1.268 0.2000
## source_PB 0.9211 2.512 0.5530 1.666 0.0960
##
## exp(coef) exp(-coef) 2.5% 97.5%
## age 0.982 1.019 0.959 1.005
## sex_F 0.965 1.036 0.547 1.704
## dis_AML 0.624 1.604 0.343 1.134
## phase_cr1 0.332 3.009 0.159 0.695
## phase_cr2 0.361 2.773 0.180 0.724
## phase_cr3 0.481 2.078 0.155 1.490
## source_PB 2.512 0.398 0.850 7.426
##
## Num. cases = 177
## Pseudo Log-likelihood = -267
## Pseudo likelihood ratio test = 24.4 on 7 df,
library(aod)
wald.test(mod1$var,mod1$coef,Terms = 4:6) #比空模型有统计学意义
## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 14.0, df = 3, P(> X2) = 0.0029