survival包学习笔记——cox回归（一）

survival包最早1985年开始撰写并应用在生存分析，包中有许多数据集供我们学习使用（肺癌、膀胱癌，急性白血病、糖尿病等），功能多样，包括对生存数据的描述、假设检验、cox风险比例回归模型的构建及竞争风险模型的构建等

同样是cox模型的构建，下面两个函数应该怎么选择？
survival::coxph()
rms::cph()

survival包是专门处理生存数据的包，所以使用survival::coxph()一定没有错
但既然存在rms::cph()就一定有他存在的意义，rms全称regression model strategies即回归模型策略，这个包可以做列线图，而且只能依靠自己构建的模型来生成列线图。

也就是说，survival包和rms包都生成原材料比如都养猪，同时rms还做香肠，rms还只用自己家的猪。害，所以没办法，你家猪再优质，想吃香肠就只能去找rms

一、cox模型的构建

这一部分主要解决以下三个问题：
1、用什么函数构建cox回归模型
2、如何展示模型
3、如何提取模型中的参数

问题1、使用什么函数
coxph(formula, data=, weights, subset, na.action, init, control, ties=c("efron","breslow","exact"), singular.ok=TRUE, robust, model=FALSE, x=FALSE, y=TRUE, tt, method=ties, id, cluster, istate, statedata, nocenter=c(-1, 0, 1), ...)
formula: 构建回归模型的公式，结局变量必须以Surv(time, status) 的形式书写，~后面是因变量，可以是一个（单因素），也可以是多个变量（多因素）
data: 数据集

我们使用survival::lung这个数据集，首先展示以下数据

> library(survival)
> lung <- survival::lung
> str(lung)
'data.frame':   228 obs. of  10 variables:
 $ inst     : num  3 3 3 5 1 12 7 11 1 7 ...
 $ time     : num  306 455 1010 210 883 ...
 $ status   : num  2 2 1 2 2 1 2 2 2 2 ...
 $ age      : num  74 68 56 57 60 74 68 71 53 61 ...
 $ sex      : num  1 1 1 1 1 1 2 2 1 1 ...
 $ ph.ecog  : num  1 0 0 1 0 1 2 2 1 2 ...
 $ ph.karno : num  90 90 90 90 100 50 70 60 70 70 ...
 $ pat.karno: num  100 90 90 60 90 80 60 80 80 70 ...
 $ meal.cal : num  1175 1225 NA 1150 NA ...
 $ wt.loss  : num  NA 15 15 11 0 0 10 1 16 34 ...

接下来使用单因cox回归筛选变量，简单先随机挑选一个变量

> cfit1 <- coxph(Surv(time, status) ~ age, data=lung)
> print(cfit1, digits=3)

image.png

> summary(cfit1, digits=3)

image.png

可以看到，直接展示模型和summary模型的输出存在差异，总的来说：

问题2、如何展示模型，如上直接输出或summary

summary（fit）输出的结果更全面，包括模型中参数：系数、HR即P值，还包括对模型的评价C index、似然比检验、Wald检验、Score值等

可见模型构建并不难，但如何解读结果才是关键。

问题3、如何提取想要的参数

coef或coefficient：提取模型中变量的系数
exp(coef())：可以求得HR值
confint：提取模型中变量的系数的95%置信区间
exp(confint())：HR值的95%置信区间
concordance：模型的一致性指数（C index）

image.png

---

二、单因素和多因素cox模型

这一部分重点关注如下问题：
1、如何构建单因素、多因素模型？
2、实战：如何批量筛选单因素回归有意义的变量？

问题1、单因素和多因素cox回归即纳入单个变量还是多个变量

单因素
如上个部分对年龄做单因素cox回归，这里不再赘述

多因素
还使用lung数据集，将全部变量纳入cox回归分析

cfit2 <- coxph(Surv(time, status) ~ ., data=lung)
print(cfit2, digits=3)
summary(cfit2, digits=3)

image.png
image.png

可以看到在多因素回归的结果中，sex, ecog, karno, wt.loss的p值有意义，在根据他们的临床意义进一步评价model ,女性、wt.loss是保护因素，ecog及karro是肺癌患者预后的危险因素。
但模型的C指数只有0.648，说明模型拟合并不是特别好的。

问题2、如何批量做单因素cox回归？

批量处理数据，这是R的优势，其他统计软件如SPSS所无法做到的
在这里要弄清楚，我们使用什么结果的筛选单因素，使用HR和P值
之前我写过一篇简书，使用循环挑选变量2万个基因cox回归分析 - 简书 (jianshu.com)

临床指标与其类似，我们来操作以下

####实战，循环筛选单因素
##为了方便，将time和event放在数据的前两列
library(dplyr)
precox <- lung %>% 
  select("time","status",everything())
precox <- precox[,-3]
precox$status <- ifelse(precox$status==1,0,1)
precox$sex <- factor(precox$sex)
precox$ph.ecog <- factor( ifelse(precox$ph.ecog<2,0,1))
precox$ph.karno
str(precox)

unicox <- data.frame();a <- c()
for(i in colnames(precox)[3:ncol(precox)]){
  fit <- coxph(Surv(time, status)~get(i),data=precox)
  fit1 <- summary(fit)
  coef <- coef(fit);HR <- exp(coef(fit));
  HRl_5 <- exp(confint(fit))[1];HRh_95 <- exp(confint(fit))[2];p_value <- fit1$coefficients[,5]
  unicox <- rbind(unicox,
                  cbind(i,coef,HR,HRl_5,HRh_95,p_value))
}
unicox[,2:6] <- apply(unicox[,2:6],2,as.numeric)
write.csv(unicox,file = "unicox.csv")

这里我对数据集lung里面的变量进一步处理，因子变量数值型变量分别处理，保留了需要处理的变量。如果你看的比较仔细，会发现status是数值型变量，请问可否设置为因子型变量？

进行cox回归时，status只能是数值型，而且必须为1代表死亡，0代表删失
这个我在一篇简书中探讨过生存曲线绘制出现Error in data.frame(..., check.names = FALSE) : arguments imply differing number of r... - 简书 (jianshu.com)

数据展示： image.png

这里还有一个问题，如何筛选变量？
一般文章中比较常见的是根据单因素、多因素的结果以及临床意义筛选变量

结果展示： image.png

C指数0.651，模型结果很一般

先暂告一段，本想接着谈谈模型效能评估及结果可视化，但内容还有很多，所以分开来写。