R语言与统计-6：生存分析和COX回归

---

R语言与统计-1：t检验与秩和检验
R语言与统计-2：方差分析
R语言与统计-3：卡方检验
R语言与统计-4：线性回归分析与模型诊断
R语言与统计-5：Logistic回归

---

更详细的可以参考之前分享的：
Kaplan-Meier生存分析的结果解读和绘制方法
TCGA生存模型的构建以及模型预测和评估

1. 生存分析

导入数据

library(coin)
data(glioma)
head(glioma)
#   no. age    sex histology group event time
# 1   1  41 Female    Grade3   RIT  TRUE   53
# 2   2  45 Female    Grade3   RIT FALSE   28
# 3   3  48   Male    Grade3   RIT FALSE   69
# 4   4  54   Male    Grade3   RIT FALSE   58
# 5   5  40 Female    Grade3   RIT FALSE   54
# 6   6  31   Male    Grade3   RIT  TRUE   25

生存分析

library(survival)
g3 <- subset(glioma,histology=='Grade3')
fit <- survfit(Surv(time,event)~group,data=g3)
plot(fit,lty = c(2,1),col = c(2,1))
legend('bottomright',legend = c('Control','Treatment'),lty=c(2,1),col=c(2,1))

进行logrank检验

survdiff(Surv(time,event)~group,data=g3)
# Call:
# survdiff(formula = Surv(time, event) ~ group, data = g3)

#                N Observed Expected (O-E)^2/E (O-E)^2/V
# group=Control  6        4     1.49      4.23      6.06
# group=RIT     11        2     4.51      1.40      6.06

#  Chisq= 6.1  on 1 degrees of freedom, p= 0.01 

也可以使用logrank_test()函数

logrank_test(Surv(time,event)~group,data=g3,distribution='exact')

#   Exact Two-Sample Logrank Test

# data:  Surv(time, event) by group (Control, RIT)
# Z = -2.1711, p-value = 0.02877
# alternative hypothesis: true theta is not equal to 1

考虑不同histology的情况

logrank_test(Surv(time,event)~group|histology,data=glioma,distribution=approximate(nresample=1000))

    Approximative Two-Sample Logrank Test

data:  Surv(time, event) by
     group (Control, RIT) 
     stratified by histology
Z = -3.6704, p-value = 0.001
alternative hypothesis: true theta is not equal to 1

2. Cox回归

COX回归模型，又称“比例风险回归模型 (proportional hazards model，简称Cox模型)”，是由英国统计学家D.R.Cox (1972)年提出的一种半参数回归模型。该模型以生存结局和生存时间为因变量，可同时分析众多因素对生存期的影响，能分析带有截尾生存时间的资料，且不要求估计资料的生存分布类型。

鉴于临床数据的特殊性，COX回归比起一般的多重线性回归和Logistic回归在临床研究中具有更为广泛的应用。

COX回归与另外两者的共同点：
(1) 分析目的：均可用于研究自变量影响程度，校正混杂因素以及作预测分析。
(2) 自变量类型：均可为连续型数值变量，或离散型分类变量，顺序变量。
(3) 自变量筛选：均可采用逐步回归方法筛选变量。

COX回归与另外两者的相异点：
(1) 因变量类型：COX回归因变量为生存资料（包含二分类结局变量和连续型生存时间变量），而多重线性回归因变量为数值型变量，Logistic回归因变量为分类或顺序型变量。
(2) 因变量数据分布：COX回归对于数据分布不作要求，而多重线性回归，Logistic回归则要求数据分布分别近似正态分布和二项分布。
(3) 模型参数解释：当自变量增大一个单位时，对于COX回归，改变了风险比HR的自然对数值。对于Logistic回归，改变了优势比OR的自然对数值。而对于多重线性回归，改变的即是Y值本身。
(4) 是否允许数据删失：COX回归允许删失值，而多重线性回归，Logistic回归不允许。

library(survival)
library(tableone)
data('GBSG2',package='TH.data')
head(GBSG2)
#   horTh age menostat tsize tgrade pnodes progrec estrec time cens
# 1    no  70     Post    21     II      3      48     66 1814    1
# 2   yes  56     Post    12     II      7      61     77 2018    1
# 3   yes  58     Post    35     II      9      52    271  712    1
# 4   yes  59     Post    17     II      4      60     29 1807    1
# 5    no  73     Post    35     II      1      26     65  772    1
# 6    no  32      Pre    57    III     24       0     13  448    1

plot(survfit(Surv(time,cens)~horTh,data = GBSG2),lty = c(2,1),col = c(2,1),mark.time = T)
legend('bottomright',legend = c('Control','Treatment'),lty=c(2,1),col=c(2,1))

coxreg <- coxph(Surv(time,cens)~.,data=GBSG2)
summary(coxreg)

#Tips：其中coef是估计得到的两个组之间的风险比的对数-ln(RR)，因此真正的风险比是exp(coef)，其实得到的就是RR（相对危险度）。最后三行是三种检验结果，在大样本的时候，这三种检验结果是一样，但是对于小样本也许会出现差别。

ShowRegTable(coxreg,digits= 2) #计算95%CI和置信区间 (上表的简化)
#              exp(coef) [confint] p     
# horThyes     0.71 [0.55, 0.91]    0.007
# age          0.99 [0.97, 1.01]    0.309
# menostatPost 1.29 [0.90, 1.86]    0.159
# tsize        1.01 [1.00, 1.02]    0.048
# tgrade.L     1.74 [1.20, 2.52]    0.004
# tgrade.Q     0.82 [0.64, 1.04]    0.099
# pnodes       1.05 [1.03, 1.07]   <0.001
# progrec      1.00 [1.00, 1.00]   <0.001
# estrec       1.00 [1.00, 1.00]    0.661

COX回归模型系数意义：
在COX回归模型中，取某一自变量系数为e的幂数，得到的值即为HR值。考虑HR值在临床研究中的实际意义，则当系数大于0（HR>1)时，该自变量为危险因素；当系数小于0（HR<1)时，该自变量为保护因素。

library(party)
tree <- ctree(Surv(time,cens)~.,data=GBSG2)
plot(tree)