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

前面我们讲过一个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_reswrite.table(file="multivariate_cox_result.txt",multi_res,quote=F,sep="\t")

                得到的结果如下

                参考资料：

                1. 一个R函数搞定风险评估散点图，热图

                赞 (0)

                相关推荐