转自个人微信公粽号【易学统计】的统计学习笔记：R语言与竞争风险模型
题记：本文是个人的读书笔记，仅用于学习交流使用。
01解决什么问题
研究治疗方案A和肿瘤复发的关系，如果患者在复查的路上出车祸意外死亡，就观察不到肿瘤复发了，也就是说“车祸死亡” 和“复发”存在竞争。这种与复发无关的死亡，失访等，就称为竞争事件。
当存在竞争事件时，传统的做法，通常是将竞争事件的发生视作刪失（censor），来计算结局事件的生存概率，采用cox模型分析，但这会导致结局事件的生存概率估计有偏差，一般会高估累计风险发生率，因此在存在竞争风险时，我们优先使用竞争风险模型，而不是cox等其他风险模型。
竞争风险的单因素分析常用来估计关心终点事件的发生率，多因素分析常用来探索预后影响因素及效应值。
02案例分析
本案例数据来自R自带数据集(survival包的mgus2)，共计1341名单克隆丙种球蛋白病患者，结局事件定义为“发生浆细胞恶性肿瘤”，而某些患者在“发生浆细胞恶性肿瘤”之前因为其他原因死亡，那这些发生其他死亡的患者就无法观察到“发生浆细胞恶性肿瘤”的终点，也就是说“其他死亡”与“发生浆细胞恶性肿瘤”存在竞争风险。故采用竞争风险模型分析。
2.1 加载数据
加载竞争风险模型的包cmprsk，加载数据，并定义结局为因子变量
library(cmprsk)
library(survival)
dt<- mgus2
str(dt)
dt$etime <- ifelse(dt$pstat==0, dt$futime, dt$ptime)
dt$event <- ifelse(dt$pstat==0, 2*dt$death, 1)
dt$event <- factor(dt$event, 0:2)
2.2 数据说明
hgb 血红蛋白；
creat 肌酐；
ptime 直至发展为浆细胞恶性肿瘤（PCM）或最后一次访视的时间（以月为单位）；
pstat 出现PCM(浆细胞恶性肿瘤)：0 =否，1 =是；
futime 直到死亡或最后一次接触的时间，以月为单位；
death 发生死亡：0 =否，1 =是；这里把PCM作为结局事件，death作为竞争事件。
2.3 单因素分析
类比两组生存资料log-rank检验，考虑竞争风险事件，同样可以进行单因素分析与多因素分析。下面我们就可以使用cuminc()函数进行单因素的Fine-Gray检验
fit0 <- cuminc(dt$etime,dt$event,dt$sex)
fit0
plot(fit0,xlab = 'Month', ylab = 'CIF',lwd=2,lty=1,
col = c('red','blue','black','forestgreen'))
结果解读：
Test第一行统计量=1.194508, P=0.274422157,表示在控制了竞争风险事件（即第二行计算的统计量和P值）后，男性和女性的累计发展为PCM的风险无统计学差异P=0.274422157。
est表示估计的各时间点男性和女性的累计发展为PCM发生率与累计竞争风险事件发生率（分别用1和2来区分，与第一行第二行一致）。
var表示估计的各时间点男性和女性的累计发展为PCM发生率与累计竞争风险事件发生率的方差（分别用1和2来区分，与第一行第二行一致）。
绘制曲线：
下面我们画出累计PCM发生率与累计竞争风险事件发生率的生存曲线，直观表示上述数字化结果。
图形解读：
纵坐标表示累计发生率pcm，横坐标是时间轴。我们从F1对应的红色曲线和M1对应的蓝色曲线可以得出，女性组的发生pcm风险较男性组高，但未达到统计学意义，P=0.27。同理，F2对应的黑色曲线在M2对应的草绿色曲线下方，我们可以得出，女性组的竞争风险事件发生率较男性组低，具有统计学差异，P=0.0064。简单来讲，这个图可以用一句话来概括：在控制了竞争风险事件后，女性组和男性组累计PCM发生率无统计学差异P=0.27。
2.4 竞争风险模型（多因素分析）
下面我们进行考虑竞争风险事件的生存资料的多因素分析方法。在cmprsk包中crr()函数可以很方便的实现多因素分析。
###多因素竞争风险
cov <- data.frame(age = dt$age,
sex = ifelse(dt$sex=='M',1,0), ## 设置哑变量
hgb = dt$hgb)
##构建多因素的竞争风险模型。此处需要指定failcode=1, cencode=0, 分别代表结局事件赋值1与截尾赋值0，其他赋值默认为竞争风险事件2。
fit2 <- crr(dt$etime, dt$event, cov, failcode=1, cencode=0)
summary(fit2)
结果解读：
在控制了竞争风险事件后，age，即疾病所处阶段是发生pcm的独立影响因素，P=0.0025。其HR及95% CI分别为0.983(0.971,0.994)，说明年龄是个保护因素，年龄越大，发生PCM的风险越低。
2.5 预测
predict_crr <- predict.crr(fit2, data.frame('age'=c(80),'sex'=c(1),'hgb'=c(15.2)))
plot(predict_crr,col='green')
2.6 与cox模型对比
可以看出，在存在竞争风险时，若不适时使用竞争风险模型，而使用cox等模型，会高估Cumulative event rate。
如果您觉得有用，请点赞，转发哦~
更多统计小知识，请关看 公粽号 易学统计
更多阅读
Nomo图出现这种问题，如何修改！
线性回归异方差问题，如何修正！
R语言：样条回归构建
R语言|广义相加模型(GAM)
R语言：多水平统计模型
R语言：广义估计方程(GEE)
R语言|两因素重复测量方差分析
R语言|基于Cox模型pec包深度验证
R语言|中位生存时间列线图绘制
R语言|Cox模型校准度曲线绘制
R语言|中位生存时间列线图绘制
基于Lasso回归筛选变量构建Cox模型并绘制Nomogram
R语言Logistic回归模型验证及Nomogram绘制
如何进行高维变量筛选和特征选择(一)？Lasso回归
