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，模型结果很一般
