玩转生存分析,这一篇就够了
1. 导入数据
library(tidyverse)
library("survival")
library("survminer")test_data=read.table("survival.txt",header = T,sep = "\t",stringsAsFactors = F)
test_data$over55=ifelse(test_data$age >= 55,1,0)
test_data$over55=as.factor(test_data$over55)
test_data$drug=as.factor(test_data$drug)head(test_data)
# studytime died drug age over55
# 1 1 1 0 61 1
# 2 1 1 0 65 1
# 3 2 1 0 59 1
# 4 3 1 0 52 0
# 5 4 1 0 56 1
# 6 4 1 0 67 1summary(test_data)
# studytime died drug age over55
# Min. : 1.00 Min. :0.0000 0:20 Min. :47.00 0:19
# 1st Qu.: 7.75 1st Qu.:0.0000 1:28 1st Qu.:50.75 1:29
# Median :12.50 Median :1.0000 Median :56.00
# Mean :15.50 Mean :0.6458 Mean :55.88
# 3rd Qu.:23.00 3rd Qu.:1.0000 3rd Qu.:60.00
# Max. :39.00 Max. :1.0000 Max. :67.00attach(test_data)
注意,因为我为了自动补全变量名,使用了attach()。如果没有用这个,像survfit() coxph()这些函数还需要加上data参数
2. KM生存曲线
drug_KM=survfit(Surv(studytime,died) ~ drug,data=test_data)
ggsurvplot(drug_KM,data=test_data,legend.title = "drug",pval = TRUE,pval.method=T,palette=c("red","blue"))
ggsave("drug.png",width = 5,height = 5)
over55_KM=survfit(Surv(studytime,died) ~ over55,data=test_data)
ggsurvplot(over55_KM,data=test_data,legend.title = "over55",pval = TRUE,pval.method=T,palette=c("red","blue"))
ggsave("over55.png",width = 5,height = 5)
3. COX比例风险模型
用两个变量进行Cox-PH model,都是分类变量,当然数值变量也可以做
cox.mod=coxph(Surv(studytime,died) ~ drug + over55)
#描述一下
summary(cox.mod)
# Call:
# coxph(formula = Surv(studytime, died) ~ drug + over55)
#
# n= 48, number of events= 31
#
# coef exp(coef) se(coef) z Pr(>|z|)
# drug1 -2.2632 0.1040 0.4582 -4.940 7.82e-07 ***
# over551 0.9274 2.5278 0.3935 2.356 0.0184 *
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# exp(coef) exp(-coef) lower .95 upper .95
# drug1 0.104 9.6135 0.04238 0.2553
# over551 2.528 0.3956 1.16888 5.4668
#
# Concordance= 0.784 (se = 0.039 )
# Likelihood ratio test= 31.03 on 2 df, p=2e-07
# Wald test = 26.92 on 2 df, p=1e-06
# Score (logrank) test = 34.49 on 2 df, p=3e-08
- exp(coef)就是HR,0.104表示:在控制年龄的情况下,用药的人死亡可能性是没有用药的人的0.104倍
- exp(-coef)为9.6135的解释刚好反过来,表示:在控制年龄的情况下,没有用药的人死亡可能性是用药的人的9.6135倍
- Concordance表示预测一致性
- 后面有三个假设:零假设是b1=b2=b3=…=b_k=0,表征的是模型整体的显著性
4. 森林图
一般文献里面比较关心HR,我们可以用森林图来表示
ggforest(cox.mod, main = '',data = test_data)
ggsave(filename = "cox.mod.png",device = "png",width = 20,height = 10,units = c("cm"))
有时候也能看到KM曲线里面加HR的图,这需要结合KM的结果和COX的结果。
到这里,对于paper的画图已经足够了,如果想深入了解Cox比例风险模型,可以看看下面的内容
5. 采用某个变量前后的模型比较
用LRT (likelihood ratio tests) 比较嵌套模型(一个模型含有另一个模型的全部变量),零假设是两个模型没有差异
cox.mod=coxph(Surv(studytime,died) ~ drug + over55)
cox.mod2=coxph(Surv(studytime,died) ~ over55)
anova(cox.mod2,cox.mod,test = "LRT")
# Analysis of Deviance Table
# Cox model: response is Surv(studytime, died)
# Model 1: ~ over55
# Model 2: ~ drug + over55
# loglik Chisq Df P(>|Chi|)
# 1 -98.062
# 2 -83.994 28.136 1 1.131e-07 ***
p值很小,说明有差异,所以我们不能去掉drug变量
6. cox模型也能用数值型变量
cox.mod3=coxph(Surv(studytime,died) ~ age)
summary(cox.mod3)
exp(coef)等于1.08表示:每增加1岁,风险增加0.08
7. 检查COX PH model的假设
- check linearity
模型中有数值型变量
线性回归里面也有这一步,方法类似 - check the proportional hazard’s assumptions
可以简单理解为某个变量对结局事件的影响(回归得到的系数)不随时间而改变
cox.mod4=coxph(Surv(studytime,died) ~ age)
png("Martingale_residuals.png")
plot(predict(cox.mod4),residuals(cox.mod4,type = "martingale"),xlab = "fitted values",ylab = "Martingale residuals",main = "Residual Plot",las=1)
#加水平线
abline(h=0)
#画残差的拟合线
lines(smooth.spline(predict(cox.mod4),residuals(cox.mod4,type = "martingale")),col="red")
dev.off()
换一种残差检查线性,deviance residuals
png("deviance_residuals.png")
plot(predict(cox.mod4),residuals(cox.mod4,type = "deviance"),xlab = "fitted values",ylab = "Deviance residuals",main = "Residual Plot",las=1)
abline(h=0)
lines(smooth.spline(predict(cox.mod4),residuals(cox.mod4,type = "deviance")),col="red")
dev.off()
check proportional hazards assumption
using Schoenfeld test for PH, Ho: HAZARDS are prop, Ha: HAZARDS are NOT prop
结果会返回每个变量,以及整体的p值
cox.zph(cox.mod4)
# chisq df p
# age 0.662 1 0.42
# GLOBAL 0.662 1 0.42
p值较大可以接受原假设。
也可以看某个变量的系数(β)是不是随着时间而改变,如果是(也就是说HR随时间而改变),则说明为non-prop hazard
par(mfrow=c(1,1))
png("cox.zph.png")
plot(cox.zph(cox.mod4)[1])
dev.off()
detach(test_data)
玩转生存分析,这一篇就够了相关推荐
- 量化免费行情源最强对比分析--看这篇就够了
序言 很多想做量化的用户一直苦于没有稳定的行情源,我也是一个,但是其实市面上有很多免费好用的行情源,在这边给大家推荐几个我用过的,给大家做个参考 先做一下对比: INSIGHT Tushare 聚宽 ...
- python数据分析实战:生存分析与电信用户流失预测
文章目录 1.背景 1.1 生存分析.KM曲线及Cox回归 1.2 案例背景 2.AIC向前逐步回归法进行特征选择 3.Cox模型搭建 3.1 特征重要性分析 3.2 模型校准 3.3 对个体进行预测 ...
- 一篇项目走进生存分析(Survival Analysis)的世界【Python版
转载自AI Studio 项目链接https://aistudio.baidu.com/aistudio/projectdetail/3410026 开篇语 生存分析在医学研究中占有很大的比例,而且进 ...
- 数据库生存曲线_WGCNA、生存分析、ROC共同筛选biomarker
Biomarker有助于疾病诊断.判断疾病分期或者用来评价新药或新疗法在目标人群中的安全性及有效性.今天将介绍一篇新的文献:针对多形性胶质母细胞瘤,利用WGCNA筛选关键模块中的hub基因,同时结合生 ...
- 绘制pr曲线图_生存分析如何绘制事件发生累计概率曲线图?
公众号前段时间发了篇推文<ggsurvplot()函数绘制Kaplan-Meier生存曲线>用来介绍生存曲线的绘制,下面的推文内容跟这篇文章结合着看. 在生存分析中我们通常关注个体在时间t ...
- 生存分析简介:Kaplan-Meier估计器
In my previous article, I described the potential use-cases of survival analysis and introduced all ...
- 如何对一个变量数据进行正则判定_生存分析数据中的BuckleyJamesMultipleRegression Model...
一.模型简介 目前,生存分析领域,最常用的是Cox比例风险回归模型,该模型具有良好的特性,不仅可以分析各种自变量对生存时间的影响,而且对基准风险分布不作任何要求(半参数模型).Cox模型使用时要满足一 ...
- r roc曲线 语言_R语言系列6:生存分析中多重时间依赖性ROC曲线绘制 timeROC
上一篇文章,我们讲到R语言实现Cox回归生存预测模型构建,以及如何将Logistic回归中,多条ROC曲线绘制在一个图里 今天主要围绕生存分析中,预测模型验证部分,如何将多条time-depend e ...
- 快速搭建ELK7.5版本的日志分析系统--搭建篇
一.ELK安装部署 官网地址:https://www.elastic.co/cn/ 官网权威指南:https://www.elastic.co/guide/cn/elasticsearch/guide ...
最新文章
- Microbiome:所谓的“富集培养”获得的微生物真的都是被“富集”出来的吗?
- java日期时间的转化
- 衡阳市2017计算机考试,2017湖南衡阳中考各科目满分及分值公告
- ASP.NET控件开发基础5
- 【英语学习】【WOTD】 putsch 释义/词源/示例
- c 语言随机验证码原理,用C生成随机中文汉字验证码的基本原理及代码.doc
- 不懂电脑如何买电脑_买电脑交智商税?5分钟看懂笔记本电脑配置
- 50套可视化报表模板直接用,做报告不用愁了!快收藏
- 计算机启用远程桌面连接失败,开启局域网远程桌面连接不上怎么办
- android 屏幕坐标系,android 屏幕坐标总结
- w10的计算机服务在哪,w10电脑服务界面在哪里
- 音乐艺术与科技有何相关?Erkki Kurenniemi的音讯是如此
- 网吧电影服务器解决方案完全指南(二)
- 京东“百亿补贴”提前20小时上线,电商价格战开打; iPhone 15 Pro玻璃面板泄露;凹语言 0.5.0发布|极客头条
- 正方教务系统换数据库服务器,自己山寨正方教务系统数据库连接解密程序
- 使用随机森林做特征选择
- 西电杨宗凯调研计算机学院,杨宗凯调研指导研究生工作:深化研究生教育改革...
- 四. 常用EMC防护器件选型学习笔记
- 更新Windows\BIOS之后WiFi消失的解决方法
- 同事都说我卷,趁着午休我 —— 彻底熟练使用《Vue3的选项APi》
热门文章
- 带HA的VCSA版本升级
- ssm整合redis和mysql_ssm整合redis - 开源中国首席碉堡了的个人空间 - OSCHINA - 中文开源技术交流社区...
- 香港公司注册详细介绍
- Office Word中横线的增加和删除方法
- 计算机一级考试操作题在线练习,计算机一级考试Photoshop练习题(带答案)
- ARCGIS之设施农用地用地报备坐标txt格式批量导出工具(定制开发版)
- python3.7 安装pip3_Centos7 安装python3、pip3、ipython3
- Git下载到idea代码被删除后如何进行还原
- 倒计时广告关闭案例实现 js js小技巧
- 用Outlook收发gmail邮件