R|生存分析-结果整理
参考:https://mp.weixin.qq.com/s?__biz=MzIyNDI1MzgzOQ==&mid=2650393867&idx=1&sn=53cc6a3d57faf0ba30dfc3e1da6481b2&chksm=f01cafebc76b26fd5d3fe809d7868c448918336837c84d31d02627e0cefeb758551b5e21a4ad&scene=21#wechat_redirect
#载入所需的R包
library("survival")
library("survminer")
#载入并查看数据集
data("lung")head(lung)
inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss1 3 306 2 74 1 1 90 100 1175 NA2 3 455 2 68 1 0 90 90 1225 153 3 1010 1 56 1 0 90 90 NA 154 5 210 2 57 1 1 90 60 1150 115 1 883 2 60 1 0 100 90 NA 06 12 1022 1 74 1 1 50 80 513 0
#cox 回归分析
res.cox <- coxph(Surv(time, status) ~ sex, data = lung)res.coxsummary(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 .95sex 0.588 1.701 0.4237 0.816Concordance= 0.579 (se = 0.022 )Rsquare= 0.046 (max possible= 0.999 )Likelihood ratio test= 10.63 on 1 df, p=0.001111Wald test = 10.09 on 1 df, p=0.001491Score (logrank) test = 10.33 on 1 df, p=0.001312
COX回归的结果中需要提取HR,HR的置信区间,wald.test和 p.value的信息,最简单的是在summary结果中进行复制粘贴,当然效率很低。假设当变量成百上前后,会发生什么呢?
还可以构建自定义函数,数据框的形式一次输出所有变量的COX回归结果
#查看待分析的变量
covariates <- names(lung[,4:10])covariates[1] "age" "sex" "ph.ecog" "ph.karno" "pat.karno" "meal.cal" "wt.loss"
#构建自定义函数,以数据框形式输出结果
univ_formulas <- sapply(covariates, function(x) as.formula(paste('Surv(time, status)~', x)))
#设定函数输出的信息
univ_models <- lapply( univ_formulas, function(x){coxph(x, data = lung)})# Extract data 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)))) })
#输出所有变量的COX结果
res <- t(as.data.frame(univ_results, check.names = FALSE))as.data.frame(res)
beta HR (95% CI for HR) wald.test p.valueage 0.019 1 (1-1) 4.1 0.042sex -0.53 0.59 (0.42-0.82) 10 0.0015ph.ecog 0.48 1.6 (1.3-2) 18 2.7e-05ph.karno -0.016 0.98 (0.97-1) 7.9 0.005pat.karno -0.02 0.98 (0.97-0.99) 13 0.00028meal.cal -0.00012 1 (1-1) 0.29 0.59wt.loss 0.0013 1 (0.99-1) 0.05 0.83
OK!可以write了,至于csv还是txt ,啦意随。。。
请关注“恒诺新知”微信公众号,感谢“R语言“,”数据那些事儿“,”老俊俊的生信笔记“,”冷🈚️思“,“珞珈R”,“生信星球”的支持!