（生物信息学）R语言与统计学入门（九）—— 单因素cox回归分析

## 近期将推出 ​​​​​​​苦逼医学生R语言入门;R语言与高级统计学 付费专栏，感谢大家支持 ##​​​​​​​

COX回归模型，又称“比例风险回归模型(proportional hazards model，简称Cox模型)”，是由英国统计学家D.R.Cox(1972)年提出的一种半参数回归模型。该模型以生存结局和生存时间为因变量，可同时分析众多因素对生存期的影响，能分析带有截尾生存时间的资料，且不要求估计资料的生存分布类型。由于上述优良性质，该模型自问世以来，在医学随访研究中得到广泛的应用，是迄今生存分析中应用最多的多因素分析方法。

下面我们直接测试代码：

 数据准备如下存成CSV格式，方便R语言读取：

setwd("D:\\")
dir()
data <- read.csv("Cox.csv",header = T,sep = ",")
head(data)


# > head(data)
#         sampleID RFS.time RFS    ALY
#1 TCGA-2A-A8VL-01      621   0 2.0320
#2 TCGA-2A-A8VO-01     1701   0 1.8710
#3 TCGA-2A-A8VT-01     1373   0 0.7952
#4 TCGA-2A-A8VV-01      671   0 2.6270
#5 TCGA-2A-A8VX-01     1378   0 2.3310
#6 TCGA-2A-A8W1-01      112   0 1.9390

下面我们下载coin软件包：

install.packages("coin")
library(coin)
cox <- coxph(Surv(RFS.time,RFS)~data$ALY,data)
cox



# > cox
#Call:
#coxph(formula = Surv(RFS.time, RFS) ~ data$ALY, data = data)
#
#           coef exp(coef) se(coef)     z      p
#data$ALY 0.2845    1.3291   0.1143 2.489 0.0128
#
#Likelihood ratio test=5.34  on 1 df, p=0.02089
#n= 421, number of events= 52 

注意我们这里研究的RFS，无病生存期。

结果显示P值是有意义的。

summary(cox)


# > summary(cox)
#Call:
#coxph(formula = Surv(RFS.time, RFS) ~ data$ALY, data = data)
#
#  n= 421, number of events= 52 
#
#           coef exp(coef) se(coef)     z Pr(>|z|)  
#data$ALY 0.2845    1.3291   0.1143 2.489   0.0128 *
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#         exp(coef) exp(-coef) lower .95 upper .95
#data$ALY     1.329     0.7524     1.062     1.663
#
#Concordance= 0.607  (se = 0.049 )
#Likelihood ratio test= 5.34  on 1 df,   p=0.02
#Wald test            = 6.2  on 1 df,   p=0.01
#Score (logrank) test = 6.08  on 1 df,   p=0.01

这里的exp(coef)就是风险比HR的意思。Lower95和upper95是HR置信区间的上下限。

这是针对连续性变量，我们研究的基因表达就是连续性变量。

cox回归也可以进行分类变量的分析，下面我们将基因的表达分为高低表达组：

data$ALY <- ifelse(data$ALY>median(data$ALY),"High expression","Low expression")
head(data)



# > head(data)
#         sampleID RFS.time RFS             ALY
#1 TCGA-2A-A8VL-01      621   0 High expression
#2 TCGA-2A-A8VO-01     1701   0  Low expression
#3 TCGA-2A-A8VT-01     1373   0  Low expression
#4 TCGA-2A-A8VV-01      671   0 High expression
#5 TCGA-2A-A8VX-01     1378   0 High expression
#6 TCGA-2A-A8W1-01      112   0  Low expression

此时是二分类变量：

cox <- coxph(Surv(RFS.time,RFS)~data$ALY,data)
summary(cox)


# > summary(cox)
#Call:
#coxph(formula = Surv(RFS.time, RFS) ~ data$ALY, data = data)
#
#  n= 421, number of events= 52 
#
#                          coef exp(coef) se(coef)      z Pr(>|z|)  
#data$ALYLow expression -0.6992    0.4970   0.2922 -2.393   0.0167 *
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#                       exp(coef) exp(-coef) lower .95 upper .95
#data$ALYLow expression     0.497      2.012    0.2803    0.8812
#
#Concordance= 0.56  (se = 0.039 )
#Likelihood ratio test= 6.04  on 1 df,   p=0.01
#Wald test            = 5.73  on 1 df,   p=0.02
#Score (logrank) test = 5.96  on 1 df,   p=0.01

其中High expression组被设置成了reference，这里的exp(coef)是针对Low expression组的，意思是这个ALY基因低表达的患者死亡风险是ALY基因高表达患者的0.497倍。

除了二分类变量以外，还可以分析多分类，无论是哪种分类变量，都会设置某一个类为一个参照。