单因素和多因素cox回归分析

生信交流平台

发布于 2022-09-21 04:58:16

3.5K00

代码可运行

文章被收录于专栏：用户7627119的专栏用户7627119的专栏

运行总次数：0

代码可运行

前面我们讲过一个R函数搞定风险评估散点图，热图，其中LASSO模型的输入就是单因素cox分析得到的显著与生存相关的基因。今天我们就来探讨一下如何使用R来做单因素和多因素cox回归分析。

我们用R的survival包自带的一套肺癌的数据来举例

#安装下面两个R包
install.packages(c("survival", "survminer"))

#加载这两个R包
library("survival")
library("survminer")

#加载肺癌这套数据
data("lung")
#显示前6行
head(lung)

这里每一行是一个样本，从第三列开始每一列是一个特征

1.单因素cox回归分析

对单个特征进行cox回归分析，看它是否与样本的生存显著相关

#单因素cox回归分析，这里看性别sex这个特征
res.cox <- coxph(Surv(time, status) ~ sex, data = lung)
res.cox

可以看到这里算出来的p值是0.00149，是显著的

我们在来看一下summary

summary(res.cox)

这里的exp(coef)就是HR（hazard ratio，风险率），lower .95和upper .95为95%的置信区间

2.批量单因素cox回归分析

一般我们的关注的特征都比较多，用上面的代码一个一个来做单因素cox回归分析效率太低了，下面我们来看看如何批量做单因素cox回归分析。

#假设我们要对如下5个特征做单因素cox回归分析
covariates <- c("age", "sex",  "ph.karno", "ph.ecog", "wt.loss")
#分别对每一个变量，构建生存分析的公式
univ_formulas <- sapply(covariates,
                        function(x) as.formula(paste('Surv(time, status)~', x)))

#循环对每一个特征做cox回归分析
univ_models <- lapply( univ_formulas, function(x){coxph(x, data = lung)})

#提取HR，95%置信区间和p值
univ_results <- lapply(univ_models,
                       function(x){ 
                         x <- summary(x)
                         #获取p值
                         p.value<-signif(x$wald["pvalue"], digits=2)
                         #获取HR
                         HR <-signif(x$coef[2], digits=2);
                         #获取95%置信区间
                         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(p.value,HR)
                         names(res)<-c("p.value","HR (95% CI for HR)")
                         return(res)
                       })
#转换成数据框，并转置
res <- t(as.data.frame(univ_results, check.names = FALSE))
as.data.frame(res)
write.table(file="univariate_cox_result.txt",as.data.frame(res),quote=F,sep="\t")

得到的结果如下，你会发现对于sex这个特征来说，结果跟前面单独做得到的结果是一样的。

3.多因素cox回归分析

前面是单独看每一个特征是否跟生存相关，而多因素cox回归是同时检测多个特征是否与生存相关。一般先通过单因素cox回归分析找出与生存显著相关的特征，然后基于这些特征再去做多因素cox回归分析，或者做LASSO分析。

res.cox <- coxph(Surv(time, status) ~ age + sex + ph.karno + ph.ecog, data =  lung)
x <- summary(res.cox)
pvalue=signif(as.matrix(x$coefficients)[,5],2)
HR=signif(as.matrix(x$coefficients)[,2],2)
low=signif(x$conf.int[,3],2)
high=signif(x$conf.int[,4],2)
multi_res=data.frame(p.value=pvalue,
                     HR=paste(HR," (",low,"-",high,")",sep=""),
                     stringsAsFactors = F
)
multi_res
write.table(file="multivariate_cox_result.txt",multi_res,quote=F,sep="\t")