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

先暂告一段,本想接着谈谈模型效能评估及结果可视化,但内容还有很多,所以分开来写。
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 220,063评论 6 510
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 93,805评论 3 396
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 166,403评论 0 357
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 59,110评论 1 295
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 68,130评论 6 395
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,877评论 1 308
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,533评论 3 420
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 39,429评论 0 276
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,947评论 1 319
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 38,078评论 3 340
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 40,204评论 1 352
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,894评论 5 347
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 41,546评论 3 331
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 32,086评论 0 23
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 33,195评论 1 272
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 48,519评论 3 375
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 45,198评论 2 357

推荐阅读更多精彩内容