## 近期将推出 苦逼医学生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倍。
除了二分类变量以外，还可以分析多分类，无论是哪种分类变量，都会设置某一个类为一个参照。
