前面我們講過一個R函數(shù)搞定風(fēng)險評估散點(diǎn)圖，熱圖，其中LASSO模型的輸入就是單因素cox分析得到的顯著與生存相關(guān)的基因浸间。今天我們就來探討一下如何使用R來做單因素和多因素cox回歸分析。
我們用R的survival包自帶的一套肺癌的數(shù)據(jù)來舉例
#安裝下面兩個R包
install.packages(c("survival", "survminer"))
?
#加載這兩個R包
library("survival")
library("survminer")
?
#加載肺癌這套數(shù)據(jù)
data("lung")
#顯示前6行
head(lung)
這里每一行是一個樣本驶臊，從第三列開始每一列是一個特征
1.單因素cox回歸分析
對單個特征進(jìn)行cox回歸分析籽暇，看它是否與樣本的生存顯著相關(guān)
#單因素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悠菜，風(fēng)險率），lower .95和upper .95為95%的置信區(qū)間
2.批量單因素cox回歸分析
一般我們的關(guān)注的特征都比較多均蜜，用上面的代碼一個一個來做單因素cox回歸分析效率太低了李剖，下面我們來看看如何批量做單因素cox回歸分析芒率。
#假設(shè)我們要對如下5個特征做單因素cox回歸分析
covariates <- c("age", "sex", "ph.karno", "ph.ecog", "wt.loss")
#分別對每一個變量囤耳，構(gòu)建生存分析的公式
univ_formulas <- sapply(covariates,
function(x) as.formula(paste('Surv(time, status)~', x)))
#循環(huán)對每一個特征做cox回歸分析
univ_models <- lapply( univ_formulas, function(x){coxph(x, data = lung)})
#提取HR，95%置信區(qū)間和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%置信區(qū)間
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)
})
#轉(zhuǎn)換成數(shù)據(jù)框偶芍，并轉(zhuǎn)置
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")
得到的結(jié)果如下充择，你會發(fā)現(xiàn)對于sex這個特征來說，結(jié)果跟前面單獨(dú)做得到的結(jié)果是一樣的匪蟀。
3.多因素cox回歸分析
前面是單獨(dú)看每一個特征是否跟生存相關(guān)椎麦，而多因素cox回歸是同時檢測多個特征是否與生存相關(guān)。一般先通過單因素cox回歸分析找出與生存顯著相關(guān)的特征材彪，然后基于這些特征再去做多因素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")
得到的結(jié)果如下
參考資料
1.一個R函數(shù)搞定風(fēng)險評估散點(diǎn)圖段化，熱圖
2.單因素和多因素cox回歸分析
