本内容为【科研私家菜】R语言机器学习与临床预测模型系列课程

你想要的R语言学习资料都在这里， 快来收藏关注【科研私家菜】

---

01 生存分析

生存分析是指对某给定事件发生的时间进行分析和推断，研究生存时间和结局与预后因子间的关系及其程度大小的方法，是一种处理删失数据的数据分析方法，也称生存率分析或存活率分析。

患者数据进行生存分析从而实现对预后和预后因子的分析，需要通过随访提前对临床数据和预后生存信息进行收集。随访通常需要数年时间对每一个随访对象进行数据采集。在随访中，关注的事件指研究中规定的生存研究的终点，在研究开始之前就已确定。根据研究目的的不同，关注的事件可以是患者死亡、疾病复发和疾病转移等。生存时间是指从某一起点到事件发生所经过的时间。生存是一个广义的概念，不仅仅指医学中的存活，也可以是疾病复发前的无病生存时间等。删失指由于关注的事件没有被观测到或者无法观测到，从而使真实生存时间无法获得的情况。删失通常由两种原因导致：（1）失访；（2）研究结束时，关注的事件尚未发生。因此，生存分析的因变量需要用生存时间和结局状态两个变量来刻画，将终点事件是否发生以及发生终点事件所经历的时间相结合。

  # log-rank 方法 
        data.survdiff <- survdiff(Surv(futime, fustat) ~ Group, data = sur_TCGA_LUAD)
        p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
        HR = (data.survdiff$obs[2]/data.survdiff$exp[2])/(data.survdiff$obs[1]/data.survdiff$exp[1])
        up95 = exp(log(HR) + qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))
        low95 = exp(log(HR) - qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))
   
        # cox 方法
        coxmodul <- coxph(Surv(futime, fustat) ~ Group, data = sur_TCGA_LUAD)
        # m=coxph(mySurv ~ group, data =  survival_dat)
        beta <- coef(coxmodul)
        se <- sqrt(diag(vcov(coxmodul)))
        HR_c <- exp(beta)
        HRse <- HR_c * se

---

02 Cox 回归方法

1）COX回归模型，又称“比例风险回归模型(proportional hazards model，简称Cox模型)”。该模型以生存结局和生存时间为因变量，可同时分析众多因素对生存期的影响，能分析带有截尾生存时间的资料，且不要求估计资料的生存分布类型。

2）原理：

假如要研究一个人从出生开始，到t时刻时死亡的概率为多大。那么它会受什么影响呢？直观的来看：

一方面，它会受时间推移的影响。一个健康的人，随着年纪的增大，他死亡的概率也会不可避免的越来越大；

另一方面，它会受一些客观因素影响，比如,一个吸烟的人在某一时刻t死亡的概率，比一个不抽烟的同龄人概率会更大；再比如，一个富豪，每年都花大价钱为自己养生、雇佣营养师为自己控制饮食起居，那么他可能就比我这个穷屌丝死亡的概率更小。

综上所述，我们抽象出了两部分的因素，一部分受时间的影响，你可以理解为是理想情况下、不受任何外界影响下的死亡的概率、是一个基准；另一部分受客观因素的影响，这些因素会影响整体的概率，使得它在基准上增加或减少。

生存分析的主要目的在于研究变量X与观察结果即生存函数(累积生存率)S(t,X)之间的关系。当S受很多因素影响，即X为向量时，传统的方法是考虑回归方程——即诸变量X 对S 的影响。但由于生存分析研究中的数据包含删失数据。且时间变量t通常不满足正态分布和方差齐性的要求，这就造成了用一般的回归方法研究上述关系的困难。

3）计算：D.R.Cox提出了Cox比例风险回归模型，它不是直接考察S 与X的关系，而是用h（t，X） 作为因变量，模型的基本形式为：

4）RR：接着，就能求出相对危险度RR：

5）Cox回归要求满足比例风险假定(proportional-hazards assumption)的前提条件。所谓比例风险假定，就是假定风险比(HR，Hazard Ratio)不随时间t变化而变化，就是等比例。另外，样本含量不易过小。

library(survival)
library("survminer")
require("survival")
library(plyr)

load("data/TCGA_LUAD_suri_exp_tumor.Rda")
View(TCGA_LUAD_suri_exp_tumor[1:100,1:100])

UniCox <- function(x){
  Group <- ifelse(coxd[,x]>= median(coxd[,x]),"High","Low")
  coxd$Group <- ifelse(coxd[,x]>= median(coxd[,x]),"High","Low")
  # Group <- factor(coxd$Group,levels = c("Low","High"))
  
  # fit_os <- survfit(Surv(OS..months., OS.censor) ~ CCT3_Group, data = sdata)
  # BaSurv <- Surv(time=sdata$OS..months.,event = sdata$OS.censor)
  BaSurv <- Surv(time=coxd$futime,event = coxd$fustat)
  # BaSurv <- Surv(time=coxd$OS_month,event = coxd$OS_censor)
  FML <- as.formula(paste0('BaSurv ~',coxd$Group))
  # FML <- as.formula(paste0('BaSurv ~',x))
  res.cox <- coxph(FML, data = coxd)
  GSum <- summary(res.cox)
  HR <- round(GSum$coefficients[,2],2)
  PValue <- round(GSum$coefficients[,5],3)
  CI <- paste0(round(GSum$conf.int[,3:4],2),collapse = '-')
  Unicox <- data.frame('Characteristics'=x,
                       'Hazard Ratio'= HR,
                       'CI95' = CI,
                       'H CI95' = paste0(HR," (",CI,")"),
                       'P Value' = PValue)
  return(Unicox)
}

03 log-rank方法

log-rank检验

需要比较两组生存曲线之间有无差别，比如某种新药组的患者生存率是否比常规药物组要高。

2×2列联表的χ2检验。这个方法有两个明显的缺点：一是由于删失的存在，难以准确计算生存率；二是时间长度可以随意指定，带来分析结果的偏差。Log-rank就是一种最常用的方法。它在每个时间点分别构建2×2列联表，然后把所有时间点结合起来分析，克服了单个时间点的缺点。

Log-rank检验的零假设是两组生存曲线一样的。如果零假设成立，那么两组内的事件发生个体数之比应该等于两组样本数之比，由此计算出事件发生的期望数。Log-rank方法就是分别将两组所有时间点的期望数加起来，与所有观察数进行比较。

为了详细说明，下面为每个时间点构建一个四格表：

如果零假设成立，那么总的事件发生数d应该按样本量等比例的分布在两组之间。那么每组的期望数（即期望的事件发生数）可以通过下式计算：

QA: 组A的所有观察数之和；

QB: 组B的所有观察数之和；

EA：组A所有的期望数之和；

EB：组B的所有期望数之和，可以简单的计算： [公式]

计算log-rank统计量：

根据P值，可以判断是否拒绝零假设，判断两组生存曲线之间有统计学差异。从km生存曲线也可以看出，两条曲线有没有差异。

mySurv=with(sur_TCGA_LUAD,Surv(futime, fustat))
View(sur_TCGA_LUAD[1:100,1:100])
exprSet <- as.matrix(t(sur_TCGA_LUAD[,-c(1,2)]))

log_rank_p <- apply(exprSet, 1 , function(gene){
  # gene=exprSet[1,]
  sur_TCGA_LUAD$group=ifelse(gene>=median(gene),'high','low')  
  sur_TCGA_LUAD$group <- factor(sur_TCGA_LUAD$group,levels = c("low","high"))
  if(length(unique(sur_TCGA_LUAD$group))>1){ 
      data.survdiff=survdiff(mySurv~group,data=sur_TCGA_LUAD)
      p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
      # p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
      HR = (data.survdiff$obs[2]/data.survdiff$exp[2])/(data.survdiff$obs[1]/data.survdiff$exp[1])
      up95 = exp(log(HR) + qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))
      low95 = exp(log(HR) - qnorm(0.975)*sqrt(1/data.survdiff$exp[2]+1/data.survdiff$exp[1]))
      
      
      tmp <- round(cbind(p = p.val,
                         HR = HR, 
                         UP95= up95,
                         LOW95= low95), 3)
      
      
      # return(p.val)
      # print(gene)
      return(tmp)
      
  }
  
})

04 参考资料

生存分析（Cox，K-M，log-rank，Fine-gray） - 哔哩哔哩 (bilibili.com)

生存分析入门之二（生存曲线的假设检验Log-rank） - 知乎 (zhihu.com)

---

关注R小盐，关注科研私家菜(VX_GZH: SciPrivate)，有问题请联系R小盐。让我们一起来学习 R语言机器学习与临床预测模型