本內(nèi)容為【科研私家菜】R語言機(jī)器學(xué)習(xí)與臨床預(yù)測模型系列課程

你想要的R語言學(xué)習(xí)資料都在這里霹肝， 快來收藏關(guān)注【科研私家菜】

01 生存分析

生存分析是指對某給定事件發(fā)生的時間進(jìn)行分析和推斷，研究生存時間和結(jié)局與預(yù)后因子間的關(guān)系及其程度大小的方法，是一種處理刪失數(shù)據(jù)的數(shù)據(jù)分析方法泡一，也稱生存率分析或存活率分析创夜。

患者數(shù)據(jù)進(jìn)行生存分析從而實現(xiàn)對預(yù)后和預(yù)后因子的分析，需要通過隨訪提前對臨床數(shù)據(jù)和預(yù)后生存信息進(jìn)行收集通惫。隨訪通常需要數(shù)年時間對每一個隨訪對象進(jìn)行數(shù)據(jù)采集茂翔。在隨訪中，關(guān)注的事件指研究中規(guī)定的生存研究的終點履腋，在研究開始之前就已確定珊燎。根據(jù)研究目的的不同，關(guān)注的事件可以是患者死亡遵湖、疾病復(fù)發(fā)和疾病轉(zhuǎn)移等悔政。生存時間是指從某一起點到事件發(fā)生所經(jīng)過的時間。生存是一個廣義的概念延旧，不僅僅指醫(yī)學(xué)中的存活谋国，也可以是疾病復(fù)發(fā)前的無病生存時間等。刪失指由于關(guān)注的事件沒有被觀測到或者無法觀測到迁沫，從而使真實生存時間無法獲得的情況烹卒。刪失通常由兩種原因?qū)е拢海?）失訪；（2）研究結(jié)束時弯洗，關(guān)注的事件尚未發(fā)生旅急。因此，生存分析的因變量需要用生存時間和結(jié)局狀態(tài)兩個變量來刻畫牡整，將終點事件是否發(fā)生以及發(fā)生終點事件所經(jīng)歷的時間相結(jié)合藐吮。

  # 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回歸模型，又稱“比例風(fēng)險回歸模型(proportional hazards model，簡稱Cox模型)”谣辞。該模型以生存結(jié)局和生存時間為因變量迫摔，可同時分析眾多因素對生存期的影響，能分析帶有截尾生存時間的資料泥从，且不要求估計資料的生存分布類型句占。

2）原理：

假如要研究一個人從出生開始，到t時刻時死亡的概率為多大躯嫉。那么它會受什么影響呢纱烘？直觀的來看：

一方面，它會受時間推移的影響祈餐。一個健康的人擂啥，隨著年紀(jì)的增大，他死亡的概率也會不可避免的越來越大帆阳；

另一方面哺壶，它會受一些客觀因素影響，比如,一個吸煙的人在某一時刻t死亡的概率蜒谤，比一個不抽煙的同齡人概率會更大山宾；再比如，一個富豪鳍徽，每年都花大價錢為自己養(yǎng)生塌碌、雇傭營養(yǎng)師為自己控制飲食起居，那么他可能就比我這個窮屌絲死亡的概率更小旬盯。

綜上所述，我們抽象出了兩部分的因素翎猛，一部分受時間的影響胖翰，你可以理解為是理想情況下、不受任何外界影響下的死亡的概率切厘、是一個基準(zhǔn)萨咳；另一部分受客觀因素的影響，這些因素會影響整體的概率疫稿，使得它在基準(zhǔn)上增加或減少培他。

生存分析的主要目的在于研究變量X與觀察結(jié)果即生存函數(shù)(累積生存率)S(t,X)之間的關(guān)系。當(dāng)S受很多因素影響遗座，即X為向量時舀凛，傳統(tǒng)的方法是考慮回歸方程——即諸變量X 對S 的影響。但由于生存分析研究中的數(shù)據(jù)包含刪失數(shù)據(jù)途蒋。且時間變量t通常不滿足正態(tài)分布和方差齊性的要求拔鹰，這就造成了用一般的回歸方法研究上述關(guān)系的困難浴捆。

3）計算：D.R.Cox提出了Cox比例風(fēng)險回歸模型能曾，它不是直接考察S 與X的關(guān)系服赎，而是用h（t，X） 作為因變量雹食，模型的基本形式為：

4）RR：接著，就能求出相對危險度RR：

5）Cox回歸要求滿足比例風(fēng)險假定(proportional-hazards assumption)的前提條件。所謂比例風(fēng)險假定茸习，就是假定風(fēng)險比(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檢驗

需要比較兩組生存曲線之間有無差別涕刚，比如某種新藥組的患者生存率是否比常規(guī)藥物組要高。

2×2列聯(lián)表的χ2檢驗乙帮。這個方法有兩個明顯的缺點：一是由于刪失的存在杜漠，難以準(zhǔn)確計算生存率；二是時間長度可以隨意指定察净，帶來分析結(jié)果的偏差驾茴。Log-rank就是一種最常用的方法。它在每個時間點分別構(gòu)建2×2列聯(lián)表氢卡，然后把所有時間點結(jié)合起來分析锈至，克服了單個時間點的缺點。

Log-rank檢驗的零假設(shè)是兩組生存曲線一樣的译秦。如果零假設(shè)成立峡捡，那么兩組內(nèi)的事件發(fā)生個體數(shù)之比應(yīng)該等于兩組樣本數(shù)之比，由此計算出事件發(fā)生的期望數(shù)筑悴。Log-rank方法就是分別將兩組所有時間點的期望數(shù)加起來们拙，與所有觀察數(shù)進(jìn)行比較。

為了詳細(xì)說明阁吝，下面為每個時間點構(gòu)建一個四格表：

如果零假設(shè)成立砚婆，那么總的事件發(fā)生數(shù)d應(yīng)該按樣本量等比例的分布在兩組之間。那么每組的期望數(shù)（即期望的事件發(fā)生數(shù)）可以通過下式計算：

QA: 組A的所有觀察數(shù)之和突勇；

QB: 組B的所有觀察數(shù)之和装盯；

EA：組A所有的期望數(shù)之和；

EB：組B的所有期望數(shù)之和甲馋，可以簡單的計算： [公式]

計算log-rank統(tǒng)計量：

根據(jù)P值埂奈，可以判斷是否拒絕零假設(shè)，判斷兩組生存曲線之間有統(tǒng)計學(xué)差異定躏。從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党窜，F(xiàn)ine-gray） - 嗶哩嗶哩 (bilibili.com)

生存分析入門之二（生存曲線的假設(shè)檢驗Log-rank） - 知乎 (zhihu.com)

關(guān)注R小鹽，關(guān)注科研私家菜(VX_GZH: SciPrivate)借宵，有問題請聯(lián)系R小鹽幌衣。讓我們一起來學(xué)習(xí) R語言機(jī)器學(xué)習(xí)與臨床預(yù)測模型