根据上面的生存分析的介绍可以大概的了解了生存分析的概念和原理以及KM曲线的绘制。但是生存分析中COX回归的结果不容易直接输出，本文简单的介绍一种自定义函数，批量并且规则的输出结果的方式。
#载入所需的R包
library("survival") library("survminer")
#载入并查看数据集
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 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.022 ) Rsquare= 0.046 (max possible= 0.999 ) Likelihood ratio test= 10.63 on 1 df, p=0.001111 Wald test = 10.09 on 1 df, p=0.001491 Score (logrank) test = 10.33 on 1 df, p=0.001312
COX回归的结果中需要提取HR，HR的置信区间，wald.test和 p.value的信息，最简单的是在summary结果中进行复制粘贴，当然效率很低。假设当变量成百上前后，会发生什么呢？
--------------------复制粘贴N*成百上千次！！！
还可以构建自定义函数，数据框的形式一次输出所有变量的COX回归结果
#查看待分析的变量
res.cox <- coxph(Surv(time, status) ~ sex, data = lung) res.cox 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.022 ) Rsquare= 0.046 (max possible= 0.999 ) Likelihood ratio test= 10.63 on 1 df, p=0.001111 Wald test = 10.09 on 1 df, p=0.001491 Score (logrank) test = 10.33 on 1 df, p=0.001312
#构建自定义函数，以数据框形式输出结果
covariates <- names(lung[,4:10]) covariates [1] "age" "sex" "ph.ecog" "ph.karno" "pat.karno" "meal.cal" "wt.loss"
#设定函数输出的信息
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.value age 0.019 1 (1-1) 4.1 0.042 sex -0.53 0.59 (0.42-0.82) 10 0.0015 ph.ecog 0.48 1.6 (1.3-2) 18 2.7e-05 ph.karno -0.016 0.98 (0.97-1) 7.9 0.005 pat.karno -0.02 0.98 (0.97-0.99) 13 0.00028 meal.cal -0.00012 1 (1-1) 0.29 0.59 wt.loss 0.0013 1 (0.99-1) 0.05 0.83
OK！可以write了，至于csv还是txt ，啦意随。。。
