R语言生存分析03-Cox比例风险模型

作者：白介素2
相关阅读：
[R语言生存分析03-Cox比例风险模型]
R语言生存分析-02-ggforest
R语言生存分析-01
ggpubr-专为学术绘图而生（二）
ggstatsplot-专为学术绘图而生(一)
生存曲线
R语言GEO数据挖掘01-数据下载及提取表达矩阵
R语言GEO数据挖掘02-解决GEO数据中的多个探针对应一个基因
R语言GEO数据挖掘03-limma分析差异基因
R语言GEO数据挖掘04-功能富集分析

如果没有时间精力学习代码，推荐了解：零代码数据挖掘课程

Cox比例风险模型

  • The Cox proportional-hazards model (Cox, 1972) Cox比例风险模型是常用与医学研究中的回归模型，通常用于研究预测变量与生存时间的关系
  • Kaplan-Meier curves 与 logrank test tests属于单因素分析的例子，他们研究的是单一变量与生存的关系,并且Kaplan-Meier 与log-rank检验只适用于分类变量，却并不适用于数值型变量，比如我们常见的基因表达
  • Cox比例风险回归分析的优势在于可以分析分类变量与数值变量，并且将生存分析的范围由单因素拓展到多因素的分析.

Cox比例风险模型基础知识

  • 在临床研究中，存在许多情况，其中几个已知量（称为协变量）可能影响患者预后。
  • 例如，假设比较两组患者：具有和没有特定基因型的患者。如果其中一组也包含较老的个体，则存活率的任何差异可归因于基因型或年龄或两者。因此，在研究与任何一个因素相关的生存时，通常需要调整其他因素的影响。
  • 统计模型是一种常用工具，可以同时分析几个因素的生存。另外，统计模型提供每个因素的效果大小。
  • cox比例风险模型是用于对生存分析数据建模的最重要方法之一。
  • 该模型的目的是同时评估几个因素对生存的影响。换句话说，它允许我们检查特定因素如何影响特定时间点发生的特定事件（例如，感染，死亡）的发生率。这个比率通常被称为危险率。 预测变量（或因子）通常在生存分析文献中称为协变量。

h(t)=h0(t)×exp(b1x1+b2x2+...+bpxp)公式可以描述为t代表生存时间，h(t)是风险函数,由协变量x1,x2,x3等决定,b1,b2,b3等代表的是协变量的系数,即该变量对生存时间影响的大小,h0(t)表示的是基线风险
exp(bi)指的是Hazard Ratio,HR,

  • HR = 1: 无效应
  • HR < 1: 风险减少
  • HR > 1: 风险增加

实例代码演示

Sys.setlocale('LC_ALL','C')
library("survival")
library("survminer")

计算风险模型的函数 coxph()

survival::coxph(formula, data, method)
formula是带有生存对象的线性模型，生存对象用Surv(time,event)
data是带有协变量的数据框,method是方法，一般默认efron

加载实例数据集

data("lung")
head(lung)
##   inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
## 1    3  306      2  74   1       1       90       100     1175      NA
## 2    3  455      2  68   1       0       90        90     1225      15
## 3    3 1010      1  56   1       0       90        90       NA      15
## 4    5  210      2  57   1       1       90        60     1150      11
## 5    1  883      2  60   1       0      100        90       NA       0
## 6   12 1022      1  74   1       1       50        80      513       0

单因素Cox回归分析

res.cox <- coxph(Surv(time, status) ~ sex, data = lung)
res.cox

## Call:
## coxph(formula = Surv(time, status) ~ sex, data = lung)
## 
##        coef exp(coef) se(coef)      z       p
## sex -0.5310    0.5880   0.1672 -3.176 0.00149
## 
## Likelihood ratio test=10.63  on 1 df, p=0.001111
## n= 228, number of events= 165

生成完整的结果报告

summary(res.cox)
## Call:
## coxph(formula = Surv(time, status) ~ sex, data = lung)
## 
##   n= 228, number of events= 165 
## 
##        coef exp(coef) se(coef)      z Pr(>|z|)   
## sex -0.5310    0.5880   0.1672 -3.176  0.00149 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##     exp(coef) exp(-coef) lower .95 upper .95
## sex     0.588      1.701    0.4237     0.816
## 
## Concordance= 0.579  (se = 0.021 )
## Likelihood ratio test= 10.63  on 1 df,   p=0.001
## Wald test            = 10.09  on 1 df,   p=0.001
## Score (logrank) test = 10.33  on 1 df,   p=0.001
str(summary(res.cox))##是一个list
## List of 14
##  $ call        : language coxph(formula = Surv(time, status) ~ sex, data = lung)
##  $ fail        : NULL
##  $ na.action   : NULL
##  $ n           : int 228
##  $ loglik      : num [1:2] -750 -745
##  $ nevent      : num 165
##  $ coefficients: num [1, 1:5] -0.53102 0.588 0.16718 -3.17638 0.00149
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr "sex"
##   .. ..$ : chr [1:5] "coef" "exp(coef)" "se(coef)" "z" ...
##  $ conf.int    : num [1, 1:4] 0.588 1.701 0.424 0.816
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr "sex"
##   .. ..$ : chr [1:4] "exp(coef)" "exp(-coef)" "lower .95" "upper .95"
##  $ logtest     : Named num [1:3] 10.6336 1 0.00111
##   ..- attr(*, "names")= chr [1:3] "test" "df" "pvalue"
##  $ sctest      : Named num [1:3] 10.32515 1 0.00131
##   ..- attr(*, "names")= chr [1:3] "test" "df" "pvalue"
##  $ rsq         : Named num [1:2] 0.0456 0.9986
##   ..- attr(*, "names")= chr [1:2] "rsq" "maxrsq"
##  $ waldtest    : Named num [1:3] 10.09 1 0.00149
##   ..- attr(*, "names")= chr [1:3] "test" "df" "pvalue"
##  $ used.robust : logi FALSE
##  $ concordance : Named num [1:2] 0.5786 0.0207
##   ..- attr(*, "names")= chr [1:2] "C" "se(C)"
##  - attr(*, "class")= chr "summary.coxph"

结果解释

  • Statistical significance:z给出的是waldtest的统计量，它对应于每个回归系数与其标准误差的比率(z = coef/se(coef)。wald统计量评估给定变量的β（β）系数是否在统计学上显着不同于0从上面的输出，我们可以得出结论，变量性别具有高度统计上显著的系数。
  • coef: regression coefficients，回归系数，coef正值表示高风险,在R中HR值是给出第二组与第一组的比值,比如这里的coef=-0.53意味着female Vs. male 95%CI: lower upper
  • 整体统计学意义： 最后，输出给出了模型总体显着性的三个替代测试的p值：似然比检验，Wald检验和得分数据统计.这三种方法是渐近等价的.对于足够大的N，它们将给出类似的结果。 对于小N，它们可能有所不同。对于小样本量，似然比检验具有更好的行为，因此通常是优选的。

批量进行单因素Cox生存分析

配合lapply函数进行批量分析,写出单因素批量分析的函数

covariates <- c("age", "sex",  "ph.karno", "ph.ecog", "wt.loss")

得到一个列表分别为，对每个变量构建的生存对象公式

univ_formulas <- sapply(covariates,
                        function(x) as.formula(paste('Surv(time, status)~', x)))
univ_formulas
## $age
## Surv(time, status) ~ age
## <environment: 0x00000000136a6138>
## 
## $sex
## Surv(time, status) ~ sex
## <environment: 0x00000000136aa1e0>
## 
## $ph.karno
## Surv(time, status) ~ ph.karno
## <environment: 0x0000000013778520>
## 
## $ph.ecog
## Surv(time, status) ~ ph.ecog
## <environment: 0x000000001377d908>
## 
## $wt.loss
## Surv(time, status) ~ wt.loss
## <environment: 0x000000001377eec0>

巧妙的配合lapply

univ_models <- lapply( univ_formulas, function(x){coxph(x, data = lung)})
univ_models
## $age
## Call:
## coxph(formula = x, data = lung)
## 
##         coef exp(coef) se(coef)     z      p
## age 0.018720  1.018897 0.009199 2.035 0.0419
## 
## Likelihood ratio test=4.24  on 1 df, p=0.03946
## n= 228, number of events= 165 
## 
## $sex
## Call:
## coxph(formula = x, data = lung)
## 
##        coef exp(coef) se(coef)      z       p
## sex -0.5310    0.5880   0.1672 -3.176 0.00149
## 
## Likelihood ratio test=10.63  on 1 df, p=0.001111
## n= 228, number of events= 165 
## 
## $ph.karno
## Call:
## coxph(formula = x, data = lung)
## 
##               coef exp(coef)  se(coef)     z       p
## ph.karno -0.016448  0.983686  0.005854 -2.81 0.00496
## 
## Likelihood ratio test=7.56  on 1 df, p=0.005966
## n= 227, number of events= 164 
##    (1 observation deleted due to missingness)
## 
## $ph.ecog
## Call:
## coxph(formula = x, data = lung)
## 
##           coef exp(coef) se(coef)     z        p
## ph.ecog 0.4759    1.6095   0.1134 4.198 2.69e-05
## 
## Likelihood ratio test=17.57  on 1 df, p=2.773e-05
## n= 227, number of events= 164 
##    (1 observation deleted due to missingness)
## 
## $wt.loss
## Call:
## coxph(formula = x, data = lung)
## 
##             coef exp(coef) se(coef)     z     p
## wt.loss 0.001319  1.001320 0.006079 0.217 0.828
## 
## Likelihood ratio test=0.05  on 1 df, p=0.8289
## n= 214, number of events= 152 
##    (14 observations deleted due to missingness)

提取出分析结果

  • 创建提取批量分析结果的函数
  • 提取分析结果，转换为数据框输出 注意是列表元素的提取,
  • exp(coef)指的是HR，beta值是系数
univ_results <- lapply(univ_models,
                       function(x){ 
                          x <- summary(x)
                          p.value<-signif(x$wald["pvalue"], digits=2)
                          wald.test<-signif(x$wald["test"], digits=2)
                          beta<-signif(x$coef[1], digits=2);#coeficient beta
                          HR <-signif(x$coef[2], digits=2);#exp(beta)
                          HR.confint.lower <- signif(x$conf.int[,"lower .95"],2)
                          HR.confint.upper <- signif(x$conf.int[,"upper .95"],2)
                          HR <- paste0(HR, " (", 
                                       HR.confint.lower, "-", HR.confint.upper, ")")
                          res<-c(beta, HR, wald.test, p.value)
                          names(res)<-c("beta", "HR (95% CI for HR)",           "wald.test", "p.value")
                          return(res)
                          #return(exp(cbind(coef(x),confint(x))))
                         })
class(univ_results)
## [1] "list"
str(univ_results)
## List of 5
##  $ age     : Named chr [1:4] "0.019" "1 (1-1)" "4.1" "0.042"
##   ..- attr(*, "names")= chr [1:4] "beta" "HR (95% CI for HR)" "wald.test" "p.value"
##  $ sex     : Named chr [1:4] "-0.53" "0.59 (0.42-0.82)" "10" "0.0015"
##   ..- attr(*, "names")= chr [1:4] "beta" "HR (95% CI for HR)" "wald.test" "p.value"
##  $ ph.karno: Named chr [1:4] "-0.016" "0.98 (0.97-1)" "7.9" "0.005"
##   ..- attr(*, "names")= chr [1:4] "beta" "HR (95% CI for HR)" "wald.test" "p.value"
##  $ ph.ecog : Named chr [1:4] "0.48" "1.6 (1.3-2)" "18" "2.7e-05"
##   ..- attr(*, "names")= chr [1:4] "beta" "HR (95% CI for HR)" "wald.test" "p.value"
##  $ wt.loss : Named chr [1:4] "0.0013" "1 (0.99-1)" "0.05" "0.83"
##   ..- attr(*, "names")= chr [1:4] "beta" "HR (95% CI for HR)" "wald.test" "p.value"
res <- t(as.data.frame(univ_results, check.names = FALSE))
as.data.frame(res)
##            beta HR (95% CI for HR) wald.test p.value
## age       0.019            1 (1-1)       4.1   0.042
## sex       -0.53   0.59 (0.42-0.82)        10  0.0015
## ph.karno -0.016      0.98 (0.97-1)       7.9   0.005
## ph.ecog    0.48        1.6 (1.3-2)        18 2.7e-05
## wt.loss  0.0013         1 (0.99-1)      0.05    0.83

通过单因素分析我们得出age,sex，ph.ecog,ph.karno具有统计学意义,由于wt.loss变量无统计学意义，我们在接下来的多因素Cox回归分析中跳过改变量，将四个变量纳入进入分析. 多因素分析回答的是，这些因素如何共同影响生存时间

多因素Cox回归

res.cox <- coxph(Surv(time, status) ~ age + sex + ph.ecog + ph.karno, data =  lung)
res.cox
## Call:
## coxph(formula = Surv(time, status) ~ age + sex + ph.ecog + ph.karno, 
##     data = lung)
## 
##               coef exp(coef)  se(coef)      z        p
## age       0.012868  1.012951  0.009404  1.368 0.171226
## sex      -0.572802  0.563943  0.169222 -3.385 0.000712
## ph.ecog   0.633077  1.883397  0.176034  3.596 0.000323
## ph.karno  0.012558  1.012637  0.009514  1.320 0.186842
## 
## Likelihood ratio test=31.27  on 4 df, p=2.695e-06
## n= 226, number of events= 163 
##    (2 observations deleted due to missingness)
summary(res.cox)
## Call:
## coxph(formula = Surv(time, status) ~ age + sex + ph.ecog + ph.karno, 
##     data = lung)
## 
##   n= 226, number of events= 163 
##    (2 observations deleted due to missingness)
## 
##               coef exp(coef)  se(coef)      z Pr(>|z|)    
## age       0.012868  1.012951  0.009404  1.368 0.171226    
## sex      -0.572802  0.563943  0.169222 -3.385 0.000712 ***
## ph.ecog   0.633077  1.883397  0.176034  3.596 0.000323 ***
## ph.karno  0.012558  1.012637  0.009514  1.320 0.186842    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##          exp(coef) exp(-coef) lower .95 upper .95
## age         1.0130     0.9872    0.9945    1.0318
## sex         0.5639     1.7732    0.4048    0.7857
## ph.ecog     1.8834     0.5310    1.3338    2.6594
## ph.karno    1.0126     0.9875    0.9939    1.0317
## 
## Concordance= 0.632  (se = 0.025 )
## Likelihood ratio test= 31.27  on 4 df,   p=3e-06
## Wald test            = 30.61  on 4 df,   p=4e-06
## Score (logrank) test = 31.06  on 4 df,   p=3e-06

可视化生存时间分布

survfit() 用于评估生存比例，默认各协变量的平均值 可视化基线生存时间分布

survfit(res.cox)
## Call: survfit(formula = res.cox)
## 
##       n  events  median 0.95LCL 0.95UCL 
##     226     163     329     286     364
ggsurvplot(survfit(res.cox),data = lung, color = "#2E9FDF",
           ggtheme = theme_minimal())
## Warning: Now, to change color palette, use the argument palette= '#2E9FDF'
## instead of color = '#2E9FDF'
image.png

展示感兴趣的变量如何影响生存分布

原文中的教程较复杂，且并可行，这里我只用自己的简单方法

构建surfut生存对象，与感兴趣的变量建立公式

fit<- survfit(Surv(time, status) ~ sex, data = lung)
fit
##
##Call: survfit(formula = Surv(time, status) ~ sex, data = lung)

##       n events median 0.95LCL 0.95UCL
##sex=1 138    112    270     212     310
##sex=2  90     53    426     348     550
##

绘制简单的曲线非常简单::ggsurvplot函数

ggsurvplot(fit, data = lung)
image.png

深度定制，这个前几期的内容似乎已经涉及了，就不再重复了

ggsurvplot(fit, data = lung,
 surv.median.line = "hv", # Add medians survival

 # Change legends: title & labels
 legend.title = "Sex",
 legend.labs = c("Male", "Female"),
 # Add p-value and tervals
 pval = TRUE,

 conf.int = TRUE,
 # Add risk table
 risk.table = TRUE,
 tables.height = 0.2,
 tables.theme = theme_cleantable(),

 # Color palettes. Use custom color: c("#E7B800", "#2E9FDF"),
 # or brewer color (e.g.: "Dark2"), or ggsci color (e.g.: "jco")
 palette = c("#E7B800", "#2E9FDF"),
 ggtheme = theme_bw() # Change ggplot2 theme
)
image.png

本期内容就到这里，我是白介素2，下期再见
参考资料链接

©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 204,590评论 6赞 478
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 86,808评论 2赞 381
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 151,151评论 0赞 337
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 54,779评论 1赞 277
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 63,773评论 5赞 367
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 48,656评论 1赞 281
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 38,022评论 3赞 398
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 36,678评论 0赞 258
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 41,038评论 1赞 299
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,659评论 2赞 321
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,756评论 1赞 330
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 33,411评论 4赞 321
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 39,005评论 3赞 307
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,973评论 0赞 19
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,203评论 1赞 260
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 45,053评论 2赞 350
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 42,495评论 2赞 343

推荐阅读更多精彩内容

  • 💞甜心小时光~为赋新词强说愁的甜宝
    甜宝感冒了，昨晚撒娇要求我拍拍睡觉，鼻塞，辗转反则睡不着，说：“妈妈，鼻子不通气，睡不着，蓝痩香菇，……” 隔会儿...
    瑞敏谈成长阅读 174评论 0赞 0
  • 小霸王孙策
    梁惠王最近写了一部新书《刺杀孙策》，据看过的人讲，完全颠覆了孙策在正史里面的形象，在这部书里面，孙策奸诈阴...
    深雨点阅读 811评论 2赞 1
  • 谁又会为你的矫情买单
    四月，大家常说人间四月天 其实不然，在我看来，四月是个慵懒，困乏的月份，这样的天不冷不热，总感觉让人提不起神来，这...
    若時光溫柔以待F阅读 641评论 0赞 2