精华内容
下载资源
问答
  • 多变量COX回归 R语言

    2021-10-08 15:53:57
    library("survival") library("survminer") res <- coxph(Surv(Survival_time, Death)~cluster_0 + cluster_1 + cluster_2 + cluster_3 + cluster_4 + cluster_5 + cluster_6 + cluster_7 + cluster_8 + cluster_...
    library("survival")
    library("survminer")
    res <- coxph(Surv(Survival_time, Death)~cluster_0 + cluster_1 + cluster_2 + cluster_3 + cluster_4 + cluster_5 + cluster_6 + cluster_7 + cluster_8 + cluster_9 + cluster_10 + cluster_11 + cluster_12 + cluster_13 + cluster_14 + cluster_15 + cluster_16 + cluster_17 + cluster_18 + cluster_19 + cluster_20 + cluster_21+ cluster_22 + cluster_23 + cluster_24 + cluster_25 + cluster_26 + cluster_27 + cluster_28 + cluster_29, data=datas)
    
    summary(res)
    ggforest(res, datas)
    res <- coxph(Surv(Survival_time, Death)~cluster_pro_0 + cluster_pro_1 + cluster_pro_2 + cluster_pro_3 + cluster_pro_4 + cluster_pro_5 + cluster_pro_6 + cluster_pro_7 + cluster_pro_8 + cluster_pro_9 + cluster_pro_10 + cluster_pro_11 + cluster_pro_12 + cluster_pro_13 + cluster_pro_14 + cluster_pro_15 + cluster_pro_16 + cluster_pro_17 + cluster_pro_18 + cluster_pro_19 + cluster_pro_20 + cluster_pro_21+ cluster_pro_22 + cluster_pro_23 + cluster_pro_24 + cluster_pro_25 + cluster_pro_26 + cluster_pro_27 + cluster_pro_28 + cluster_pro_29, data=datas)
    
    summary(res)
    ggforest(res, datas)
    
    展开全文
  • 生存分析之R包survival的单变量和多变量Cox回归续前文生存分析(Survival Analysis)。在前文初步简介了生存分析的概念,以及展示了一种生存分析模型Kaplan-Meier的使用。Kaplan-Meier是一种非参数的单变量分析方法,...
    696d9db39af9033044b0155592ac33b1.gif 生存分析之R包survival的单变量和多变量Cox回归 2188b8f8458fcc06e1e522ff33372006.gif 续前文生存分析 (Survival Analysis)。

    在前文初步简介了生存分析的概念,以及展示了一种生存分析模型Kaplan-Meier的使用。Kaplan-Meier是一种非参数的单变量分析方法,通过估计生存率以及中位生存时间,以生存曲线方式展示生存特征。通常,中位生存时间越长,或者曲线变化幅度平缓则代表着预后较好。Kaplan-Meier分析和对数秩检验(Log-rank test)的结合使用可用于比较两组或多组的生存曲线是否存在显著差异,评估关键的预后因素。

    另一种方法是Cox回归模型(或称Cox比例风险回归模型,Cox proportional hazards regression model),它既适用于定量预测变量也适用于类别变量。此外,Cox回归模型扩展了生存分析方法:尽管Kaplan-Meier分析非常流行,但其一种局限性是只能用于单变量分析,相比之下,Cox回归模型可分析多种风险因素对生存的影响。

    本篇继续简介Cox回归模型,及其在R语言中的计算方法。

        

    Cox回归模型简介

    Cox回归是医学研究中常用的统计回归模型,其原理是将生存时间与协变量联系起来(Cox, 1972)。例如在医学领域,期望寻求找出哪个(或哪些)协变量对患者的生存时间具有最重要的影响,Cox回归模型将建立患者的生存率和几个解释变量之间的关系,评估这些因素对生存时间的效应。

    Cox回归模型无需生存时间有任何特定的分布,其假设不同变量对生存的影响随时间变化是恒定的,并且在特定范围内是累加的。

    Cox回归执行式

    Cox回归的执行式如下:

    45ef50c53c0037223f9624e9b6a2316c.png

    式中:t表示生存时间;h(t)是由一组p个协变量(x1,x2,…,xp)确定的风险函数(hazard function);系数(b1,b2,…,bp)量化了协变量的影响(例如,效应的大小);h0称为基线风险(baseline hazard),如果所有xi=0(即exp(xi)=1),则它对应于风险值;h(t)中的“t”提醒我们,风险可能随时间而变化。

    将上述公式进行对数转换后,可以变换为形似线性方程的结构:

    34845a60e6d9e586bd5068da2029f1f9.png

    因此,Cox回归模型可以视为风险与变量xi的对数的多元线性回归,其中基线风险是一个随时间变化的“截距”项。

    偏回归系数与风险的关系

    Cox回归模型涉及检查每个自变量的系数,如b1、b2等,称之为偏回归系数。正回归系数意味着该变量值较高的患者存在高风险,即预后较差;负回归系数意味着该变量值较高的患者预后较好。

    将偏回归系数作为自然对数的幂,即eb,写作exp(b),称为风险比(hazard ratio,HR)。因此,当HR>1(即b>0),意味着该变量值较高的患者存在高风险,即预后较差;当HR<1(即b<0),意味着该变量值较高的患者预后较好。

    在临床上,也常将HR>1的自变量称之为坏的预后因子,将HR<1的自变量称之为好的预后因子。

        

    R包survival的Cox回归模型

    类似前文Kaplan-Meier分析,R包survival中也提供了执行Cox回归模型的方法,本篇继续以survival包的方法为例展示。

    示例数据

    survival包的内置数据集lung,为NCCTG(North Central Cancer Treatment Group)的肺癌患者临床资料数据。

    library(survival)
    #示例数据,详情 ?lung
    data(lung)
    head(lung)

    2b3651360faa82d6697fc47aa03e068a.png

    该数据集中记录了228例晚期肺癌患者的生存时间(time,单位为天)、生存状态(status,1为在世,2为去世)以及其它状态特征,例如年龄(age)、性别(sex,1为男性,2为女性)等,详情?lung查看帮助文档中的简介即可。

    单变量Cox回归(对于分类变量)

    首先期望比较男性与女性肺癌患者的生存率是否具有显著不同,此时可通过单变量Cox回归来解决。(单变量生存分析中的另一种常用方法Kaplan-Meier分析,可参考前文)

    指定指定数据集中的生存时间(time)、患者生存状态(status)以及分组(sex,性别)等执行单变量Cox回归。

    注:尽管分组是分类变量,但仍需要以1、2等连续的离散数值表示。

    #单变量 Cox 回归,详情 ?coxph
    #生存时间与性别(分类变量)的关系
    cox_sex cox_sex
    summary(cox_sex)
    #结果提取,例如
    #names(summary(cox_sex))
    summary(cox_sex)$coefficients  #变量的回归系数、p 值等
    summary(cox_sex)$logtest  #Likelihood ratio test 的统计量、p 值等
    summary(cox_sex)$waldtest  #Wald test 的统计量、p 值等
    summary(cox_sex)$sctest  #Score (logrank) test 的统计量、p 值等

    5928582eee24696d68a7990e2b758071.png

    结果的下方,给出了三种类型的检验p值(Likelihood ratio test、Wald test、Score (logrank) test),它们所示的是整个回归方程的显著性,均显示p=0.001,结果非常显著,表明两组(男性和女性)患者之间的生存率具有显著差异。

    评估变量的效应,查看上方的自变量的偏回归系数等指标。在单变量Cox回归中,单个自变量的显著性p值和整个回归方程p值是一致的。coef为偏回归系数,exp(coef)即代表了HR,这里HR<1,代表了随数值增大而风险降低。由于1代表了男性,2代表了女性,所以男性患者(组1)比女性患者(组2)的生存期更短,有理由认为男性肺癌患者比女性肺癌患者具有更差的预后,猜测原因可能与男性更多存在吸烟史有关(数据集中没体现)。

    构建好的Cox回归模型可进一步用于预测患者风险。

    如下示例,显示男性患者的风险明显大于女性患者。

    #通过构建好的 Cox 回归模型执行风险预测
    lung$risk_sex head(lung$risk_sex)
    library(ggplot2)
    ggplot(lung, aes(factor(sex), risk_sex, color = factor(sex), label = risk_sex)) +
    geom_point() +
    scale_color_manual(values = c('blue', 'red'), limits = c('1', '2'), labels = c('male', 'female')) +
    theme(panel.background = element_rect(fill = 'transparent', color = 'black'), 
        panel.grid = element_blank()) +
    geom_text(vjust = -0.5, show.legend = FALSE)
    4ae88c937b660cff45b20919437dc91a.png

    对于生存曲线,可通过survminer包绘制,比较男性和女性患者的生存率。

    注:以下过程实际上执行了一个Kaplan-Meier分析,借助Kaplan-Meier的结果绘制生存曲线。

    library(survminer)
    #生存曲线和累积风险曲线
    KM  
    ggsurvplot(KM, conf.int = TRUE, palette = c('blue', 'red'), risk.table = TRUE, pval = TRUE)
    ggsurvplot(KM, conf.int = TRUE, palette = c('blue', 'red'), fun = 'cumhaz')

    a6785797de5dabcf07dfce93d8821062.png

    曲线中同时显示了置信区间以及分位数时间段的生存患者数量等,便于我们比较。但需要注意的是,这里显示的p值实际上是Kaplan-Meier分析的p值,同样是非常显著的。

    通过生存曲线以及风险曲线,可以直观看出男性肺癌患者比女性肺癌患者具有更差的预后,体现在生存曲线下降更快以及风险曲线上升更快。

    单变量Cox回归(对于连续变量)

    由于通常认为癌症与年龄密切相关,即老年人口存在更高的风险,因此年龄因素也是不可忽略的。接下来试图通过单变量Cox回归来探索患者生存时间与年龄的关系。

    指定指定数据集中的生存时间(time)、患者生存状态(status)以及多个自变量(sex,性别以及age,年龄)等执行多变量Cox回归。

    #生存时间与年龄(连续变量)的关系
    cox_age cox_age
    summary(cox_age)

    a7c1156becd43f6bf183ba603442c6c2.png

    结果显示p<0.05,表明患者的生存时间确实与年龄存在关联。并由于exp(coef)>1,即代表了HR>1,表明风险随患者年龄增加而增加,高龄人群对应了更差的预后。

    因为这种连续变量的情况下难以绘制生存曲线,我们不妨绘制散点图查看年龄和生存时间的关系。(以下表示方式可能有些不妥,但只是用于探索趋势,所以散点图+简单的线性回归线似乎仍是可行的)

    #生存时间与年龄的散点图
    ggplot(lung, aes(age, time, color = factor(status), shape = factor(status))) +
    geom_point() +
    scale_shape_manual(values = c(1, 16), limits = c('1', '2')) +
    theme(panel.background = element_rect(fill = 'transparent', color = 'black'),
    panel.grid = element_blank()) +
    geom_smooth(method = lm, se = FALSE, show.legend = FALSE)

    4e67b95b7f312da71b0733a8ae52bd43.png

    通过观察,患者生存时间表现出随年龄增长而下降的趋势(关注status=2的趋势,这些是已经去世的患者;status=1是仍在世的患者,暂且忽略),同样透露着高龄人群存在更高风险。

    类似地,通过已构建好的Cox回归模型执行风险预测,同样可以清楚地看到,患者风险随年龄增加而增加。

    #通过构建好的 Cox 回归模型执行风险预测
    lung$risk_age head(lung$risk_age)
    ggplot(lung, aes(age, risk_age, color = age, shape = factor(status))) +
    geom_point() +
    scale_color_gradientn(colors = c('green', 'yellow', 'red')) +
    scale_shape_manual(values = c(1, 16), limits = c('1', '2')) +
    theme(panel.background = element_rect(fill = 'transparent', color = 'black'),
    panel.grid = element_blank())
    24901338dcafa35ba052f4f9dd3391ea.png

    多变量Cox回归

    再来看一种考虑多组变量的情况。现在我们同时考虑性别和年龄与患者生存时间的关系,通过多变量Cox回归来解决。

    指定指定数据集中的生存时间(time)、患者生存状态(status)以及将多个自变量(sex,性别以及age,年龄)等同时作为协变量执行多变量Cox回归。

    #多变量 Cox 回归,详情 ?coxph
    #多变量时,直接在执行式中以“+”相连即可
    cox2 cox2
    summary(cox2)

    155f443f5fccd5b38117b59f25b5c166.png

    类似的,结果中给出了三种类型的检验p值(Likelihood ratio test、Wald test、Score (logrank) test),它们所示的是整个回归方程的显著性。根据p≤0.001可认为回归方程是统计学显著的,说明在多个自变量中包含了对生存时间具有影响的因素。

    在多变量情况下,若评估各个自变量的相对效应,同样可关注上方偏回归系数那一栏,多变量情况下,比较各个自变量的回归p值。性别(sex)更为显著的,表明性别是患者之间生存期不同的主要因素,患者年龄次之,但作为协变量的情况下表现不明显(p>0.05)。

    随后再审查偏回归系数等指标,关注exp(coef)信息,和上文中分别使用性别(sex)和年龄(age)的单变量Cox回归具有一致的趋势。性别(sex)的exp(coef)<1(即HR<1),年龄的exp(coef)>1(即HR>1),表明男性肺癌患者比女性肺癌患者具有更差的预后,高龄人群可能存在更高的风险。

    通过已构建好的Cox回归模型执行风险预测,观察两种自变量协同作用对患者风险的效应。可以看到,男性患者(sex=1)比女性患者(sex=2)具有更高的风险,并且均表现出随年龄增加而增加的趋势。

    #通过构建好的 Cox 回归模型执行风险预测
    lung$risk_sex_age head(lung$risk_sex_age)
    ggplot(lung, aes(age, risk_sex_age, color = factor(sex), shape = factor(status))) +
    geom_point() +
    scale_color_manual(values = c('blue', 'red'), limits = c('1', '2')) +
    scale_shape_manual(values = c(1, 16), limits = c('1', '2')) +
    theme(panel.background = element_rect(fill = 'transparent', color = 'black'),
    panel.grid = element_blank())
    af2775d6e696ede183dd7f60020c4631.png

    此外,还可通过survminer包绘制森林图,比较多变量的效应。

    如下示例,对本次多变量Cox回归中的性别(sex)和年龄(age)的HR值以及p值等进行可视化。

    #森林图,详情 ?ggforest
    ggforest(cox2, main = 'hazard ratio', refLabel = 'reference', noDigits = 2)
    224d521caf5f56cac4c6b0aa1652df5d.png

    关于多变量Cox回归中自变量的选择

    对于多变量的情况,通常会考虑这样一个问题,就是选择怎样的变量组合是最合适的?以下提供几点参考。

    首先可以根据各个自变量的p值进行评估,如果p值特别地大(非常不显著的情况),则表明这个变量作为协变量时效应不明显,可以考虑剔除它。

    对比多个Cox模型的AIC值(Akaike信息准则)等信息。AIC会在给定模型的复杂性与其拟合优度之间进行权衡,可以将AIC值视为对应了模型的准确性,AIC值越小的模型表明越有可能准确地预测新数据。

    #AIC 值获取,以上述结果为例
    extractAIC(cox2)

    简约性原则也是必需考虑的,也就是避免使用太多的自变量,尽可能使模型简约更利于解释。

        

    参考资料

    https://www.xlstat.com/en/solutions/features/cox-proportional-hazards-models https://www.r-bloggers.com/cox-proportional-hazards-model/ http://dwoll.de/rexrepos/posts/survivalCoxPH.html Cox DR (1972). Regression models and life tables (with discussion). J R Statist Soc B 34: 187–220 e447183c1ec20fff05983e1eacbf2185.png 我就知道你在看
    展开全文
  • 从TCGA上下载数据库和临床数据之后,往往需要进行COX分析,一般的分析思路是先进行单变量,在进行多变量分析。然而,当关注的基因比较是,手动输入就会比较麻烦。接下来介绍一种利用循环的方法,快速的对个...

    从TCGA上下载数据库和临床数据之后,往往需要进行COX分析,一般的分析思路是先进行单变量,在进行多变量的分析。然而,当关注的基因比较多是,手动输入就会比较麻烦。接下来介绍一种利用循环的方法,快速的对多个变量进行分析。

    首先是导入数据,包括基因表达counts数据和临床数据sur,autophage是我下的一个自噬基因集,可根据需要替换为其他需要分析的基因列表,以及要用到的包:

    setwd("D:/A1/Rdata/Autophage/胰腺癌")
    library("survival")
    library("survminer")
    library("clusterProfiler")
    library("DO.db")
    
    #输入数据,data为数据矩阵,sur为生存数据,Autophage为自噬相关基因表
    data <- read.csv("data.csv",header = T)
    sur <- read.csv("sur.csv",header = T)
    autophage <- read.csv("Autophage.csv",header = T)
    
    #去掉重复值
    index <- duplicated(sur$Sample) 
    sur <- sur[!index,]

    data格式如下:

    sur数据格式如下:

    对data进行标准化:

    #log2标准化
    row <- row.names(data)
    data <- as.matrix(data)
    data <- apply(data,2, as.numeric)
    qx <- as.numeric(quantile(data, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
    LogC <- (qx[5] > 100) ||
      (qx[6]-qx[1] > 50 && qx[2] > 0)
    if (LogC) {
      data[which(data <= 0)] <- NA
      data <- log2(data+1)
    }
    data <- as.data.frame(data)
    row.names(data) <- row

    因为data中的基因格式是EMSEMBL,所以将autophage中的基因SYMBOL转换为EMSEMBL:

    autogn <- autophage$Symbol
    autogn <- bitr(autogn, fromType = "SYMBOL",toType = c("ENSEMBL"),OrgDb = "org.Hs.eg.db")

    随后进行单变量cox分析:

    基本思路是,将autophage中的基因EMSEMBL一个个取出,与data中的基因对应,并将该基因按中位数分为高表达组(1)和低表达组(0),与生存数据存入同一个数据框中(sur_r),随后进行单因素cox分析,并将结果存入unicox中。

    a <- colnames(data)
    a <- gsub("\\.","-",a) #替换字符串
    i <- length(a)
    while (i > 0) {a[i] <- substr(a[i],1,12);i <- i-1}
    r <- length(row.names(autogn))
    k <- length(data)
    data$name <- row.names(data)
    unicox <- data.frame(geneID = NA,HR = NA,downCI = NA,upCI = NA,P.val = NA)
    s = 1
    
    #自动将data中的自噬基因匹配并进行单因素cox分析
    while(r > 0){
      gene <- autogn[r,2]
      if (gene %in% row.names(data)){
        b <- data[data$name==gene,1:k]
        b <- as.numeric(b)
        ab <- data.frame(Sample = a,group = b)
        sur_r <- merge(sur,ab,by = "Sample")
        avg <- mean(b)
        sur_r$group[sur_r$group < avg] <- 0
        sur_r$group[sur_r$group >= avg] <- 1
        n <- sur_r$days_to_death
        n <- as.numeric(n)
        sur_r$days_to_death <- n
        n <- sur_r$status
        n <- as.numeric(n)
        sur_r$status <- n
        n <- sur_r$group
        n <- as.numeric(n)
        sur_r$group <- n
        res.cox <- coxph(Surv(days_to_death, status) ~ group, data = sur_r)
        res <- summary(res.cox)
        unicox[s,1] <- gene
        unicox[s,2] <- res$conf.int[1]
        unicox[s,3] <- res$conf.int[3]
        unicox[s,4] <- res$conf.int[4]
        unicox[s,5] <- res$waldtest[3]
        s = s + 1
      }
      ;r <- r-1
    }

    unicox格式如下:

    存储了基因名,HR,95%CI和P值,下面需要选出P值<0.05的基因进入多因素cox分析,最开始纳入了16个变量,经过筛选后仅剩下6个:

    #筛选有统计学意义的变量
    unicox_sub <- unicox[unicox$P.val < 0.05,]
    write.csv(unicox,"unicox_new.csv",row.names = F)
    write.csv(unicox_sub,"unicox_select_new.csv",row.names = F)
    
    #多元cox回归分析
    p_max <- 1
    muticox_list <- unicox_sub$geneID
    i <- length(muticox_list)
    l <- 19
    c <- colnames(sur_r)
    while (i > 0) {
      gene <- muticox_list[i]
      b <- data[data$name==gene,1:k]
      b <- as.numeric(b)
      avg <- mean(b)
      b[b < avg] <- 0
      b[b >= avg] <- 1
      ab <- data.frame(Sample = a,group = b)
      sur_me <- merge(sur,ab,by = "Sample")
      sur_r[,l] <- sur_me$group
      c[l] <- gene
      l <- l + 1
      ;i <- i - 1
    }
    colnames(sur_r) <- c
    l <- length(muticox_list)
    muticox_list[l+1] <- "days_to_death"
    muticox_list[l+2] <- "status"
    sur_muti <- subset(sur_r,select = muticox_list)
    res.muticox <- coxph(Surv(days_to_death, status) ~ ., data = sur_muti)
    sum_muticox <- summary(res.muticox)
    muticox_result <- cbind(sum_muticox[["conf.int"]],sum_muticox[["coefficients"]])
    muticox_result <- as.data.frame(muticox_result)
    colnames(muticox_result)[9] <- "p.val"
    p_max <- max(muticox_result$p.val)
    while (p_max > 0.15){
      muticox_result <- muticox_result[order(muticox_result$p.val,decreasing = T),]
      muticox_result <- muticox_result[-1,]
      muticox_list <- row.names(muticox_result)
      i <- length(muticox_list)
      l <- 19
      c <- colnames(sur_r)
      while (i > 0) {
        gene <- muticox_list[i]
        b <- data[data$name==gene,1:k]
        b <- as.numeric(b)
        avg <- median(b)
        b[b < avg] <- 0
        b[b >= avg] <- 1
        ab <- data.frame(Sample = a,group = b)
        sur_me <- merge(sur,ab,by = "Sample")
        sur_r[,l] <- sur_me$group
        c[l] <- muticox_list[i]
        l <- l + 1
        ;i <- i - 1
      }
      colnames(sur_r) <- c
      l <- length(muticox_list)
      muticox_list[l+1] <- "days_to_death"
      muticox_list[l+2] <- "status"
      sur_muti <- subset(sur_r,select = muticox_list)
      res.muticox <- coxph(Surv(days_to_death, status) ~ ., data = sur_muti)
      sum_muticox <- summary(res.muticox)
      muticox_result <- cbind(sum_muticox[["conf.int"]],sum_muticox[["coefficients"]])
      muticox_result <- as.data.frame(muticox_result)
      colnames(muticox_result)[9] <- "p.val"
      p_max <- max(muticox_result$p.val)
    }
    
    gene_df <- row.names(muticox_result)
    gene_df <- bitr(gene_df, fromType = "ENSEMBL",toType = c("SYMBOL"),OrgDb = "org.Hs.eg.db")
    row.names(muticox_result) <- gene_df$SYMBOL
    write.csv(muticox_result,"muticox_result_new.csv")

    代码比较长,基本思路是先进行一次多因素cox分析,得到模型中P值最大的变量,随后将其存在p_max中,并将改变了剔除,随后利用循环不断剔除p_max最大的变量,直到p_max<0.15,也就是多变量分析的向后选择法,剔除标准为0.15,随后得到结果:

    最终剩余6个变量,对这些变量绘制森林图:

    muticox_list[1:l] <- gene_df$SYMBOL
    colnames(sur_muti) <- muticox_list
    res.muticox <- coxph(Surv(days_to_death, status) ~ ., data = sur_muti)
    
    #绘制森林图
    tiff(file="muticox_result_new.tiff",
         width = 20,          
         height = 20,            
         units ="cm",
         compression="lzw",
         bg="white",
         res=500)
    ggforest(res.muticox)
    dev.off()

     以上便是全过程了,巧用循环语句,可以减少生信分析中很多重复的工作。

    最后还有一点疑问,如何根据上述的结果构建预测癌症预后的风险模型呢?欢迎大家参与讨论!

    展开全文
  • 此下载内容应包括此代码,用于测试目的的随机生成的excel数据文件以及示例输出。 该程序需要四个 excel 列 1. 跟进时间和 2. 失败指标(1 = 发生了失败,0 = 没有发生失败。3. Predictor1 和 4....
  • 转自个人微信公众号【Memo_Cleon】的统计学习笔记:生存分析Cox回归。随访资料的生存分析是一个很大的题目。从分析的因素上看,有单因素分析因素分析。正如“连续资料的单因素分析常用t检验、方差分析,对应的...

    转自个人微信公众号【Memo_Cleon】的统计学习笔记:生存分析之Cox回归

    随访资料的生存分析是一个很大的题目。

    从分析的因素上看,有单因素分析和多因素分析。正如“连续资料的单因素分析常用t检验、方差分析,对应的多因素分析是多重线性回归”、“分类资料的单因素分析方法卡方分析,对应的多因素分析有logistic回归”一样,生存分析的常用单因素(或少数因素)的分析有Life Tables法、Kaplan-Meier法,对应的多因素模型则常用Cox回归模型(Cox风险比例模型)。从采取的分析方法上看,生存分析有非参数法(如Wilcoxon法、Log-rank法)、参数法(如Weibull回归、lognormal回归等)和半参数分析(Cox回归)。

    Cox回归要求满足比例风险假定(proportional-hazards assumption)的前提条件。所谓比例风险假定,就是假定风险比(HR,Hazard Ratio)不随时间t变化而变化。

    在进行生存分析前,你最好对以下的一些概念及其意义有所了解:起始事件、失效事件(Failure Event)/终点事件(Endpoint Event)、生存时间(Survival Time)/失效时间(Failure Time)、中位生存时间(Median Survival Time)、平均生存时间(Mean Survival Time)、删失值/截尾值(censored values)、生存概率(Survival Probability)、生存率(Survival Rate)/积累生存概率/生存函数/积累生存函数(Cumulative Survival Function)、风险函数(Hazard Function)、累积风险函数……风险函数h(t)=概率密度函数f(t)/生存函数S(t),概率密度函数f(t)为累积分布函数F(t)的导数,而F(t)=1-S(t)。可参见《生存分析》。

    模型结构与参数释义可参见颜虹等主编的《医学统计学》,如下。对此不感兴趣而只关心操作和结果解读的,可直接越过。

    eca0e4eea600c67dcba268d7341a52a4.png

    当前笔记用STATA演示Cox回归操作。STATA在进行Cox回归分析前首先需要声明生存时间变量,另外比例风险假定是进行Cox回归的前提条件,需要进行考察和检验。

    示例(陈启光等.医学统计学第3版):探讨某肿瘤的预后,某研究机构收集了41例患者的生存时间(月份)、生存结构及影响因素。影响因素包括性别、年龄、病理分级、是否复发、PD-L1分子。变量赋值与资料如下:

    50619e5ad25703d9a3d3316da76a517a.png

    60baa5fec72bbfbfa671b9fc3e88b75b.png

    【1】数据录入:

    【2】声明时间变量:统计>>生存分析>>模型设定和实用工具>>声明数据集为生存时间数据,将在[时间变量]和[失效变量]中选择相应的变量即可。具体步骤如下:

    ab75022ca0e21dbdc4406b2f83969c45.png

    相应命令如下:

    stset time, failure(status==1)

    也可以不指定失效值,默认不等于0的为失效值。另外scale选项可以重新定义生存时间,如本例生存时间单位是月,可以使用scale(12)后生存时间代表的就是年啦,命令为:stset time, failure(status) scale(12)

    9cad9013ee1f2d231db0f90d3eb14192.png

    【3】Cox回归统计>>生存分析>>回归模型>>Cox比例风险模型,选入相应的自变量即可。主要添加自变量时应选择合适的变量类型即参照水平(默认低水平为参照),如果未经步骤【2】的生存时间变量声明,也可以在此对话框中的[生存设置]中声明。需要说明的是,本例有序变量按连续变量处理,性别、病例分级、复发即PD-L1虽然是分类变量但都是二分类,直接按连续变量分析也不影响结果。命令为stcox i.gender age i.grade i.recur i.PD_L1 或者stcox gender age grade recur PD_L1

    abef9e948adf6d64512472502bacbda3.png

    结果显示相比空白模型,纳入五个变量的整体模型是有统计学意义的(LR Chi2=22.3,P=0.0005)。主要结果中默认显示的是风险比(HR),如果想显示模型系数β,需要在Cox比例风险模型对话框中的[报告]选项卡中选中复选框[报告系数,而不是风险比],HR=exp(β)。根据系数可得出模型为:h(t)=h0(t)exp(-0.113·gender+0.365·age-5.44·grade+1.985·recur+0.190·PD-L1)结果显示校正其他因素,女性相比男性、病例IV期相比III期死亡风险更低,但没有统计学意义;而年龄每增加一个等级、复发相对不复发、PD-L1阳性相比阴性死亡风险更高(分别是1.44、7.28、1.21倍),但只有复发与否具有统计学意义。

    cf8fdc0139ba61c438057d5a0a6f98fc.png

    08725399b25602ce347d465a1e5d1427.png

    很明显,模型给出是强制纳入了所有的变量的结果,但实际上有很多因素并不具有统计学意义,为了精简模型,我们可能需要对模型的变量进行筛选,可以采用逐步回归,SPSS里面可直接在[method]中进行选择相应的方法。STATA里面则可借助stepwise命令,Cox逐步回归菜单操作:

    统计>>其他>>逐步估计,具体操作如下。本例采用向前逐步回归,默认wald检验,变量剔除标准P值为0.1,纳入标准P值为0.05。如果直接进行逐步回归,分析前不要忘记声明生存时间变量。

    b7cb4e42b5f898bac7003269a394c318.png

    相应命令:stepwise, pr(0.1) pe(0.05) forward: stcox gender age grade recur PD_L1

    说明

    stepwise, pr(0.1) pe(0.05) forward : stcox (gender) (age) (grade) (recur) (PD_L1), nohr

    pr(剔除标准),pe(纳入标准)。逐步回归方法:forward、backward,逐步回归方法:lr,默认是wald法。括号内为同进同出的变量为一个回归项,类似于SPSS里面的Block,同一个回归项中变量同进同出,比如同一个分类变量设置为哑变量后可以放在同一个回归项中。nohr表示报告系数而不是风险比。

    逐步回归结果如下:最终引入的变量是recur和age,模型:h(t)=h0(t)exp(1.756·recur+0.451·age)recur和age的HR分别为5.791、1.569,即是否复发与年龄是该肿瘤的死亡风险因素。固定其他因素的影响,患者肿瘤复发的死亡风险是不复发的5.791倍;固定其他因素的影响,患者年龄每增加一个等级,死亡风险增加1.569倍。

    a365ff7297f8e1976cf13b948cfda504.png

    有时候我们还想得到生存曲线。以指定协变量recur,绘制生存率曲线为例,操作如下:

    [图形>>生存分析图>>生存,风险,累积风险或累积发生函数]或者[统计>>后验估计],在打开的[后验估计选择器]中依次选择:设定,诊断和拟合优度分析>>生存函数,风险函数,累积风险等图形,点击[开始]进入[曲线-绘制生存函数,风险函数,累积风险函数或累积发病率函数]对话框,界面操作如下:

    441d9fb4a6cb01f83fa8a80f9bd7ad10.png

    相应命令为:stcurve, survival at1( recur=0 ) at2( recur=1 )

    结果如下(你要是觉得图片丑可以启用图形编辑器进行编辑):

    20eb7d72dc9f6ba82074db97b65b32f4.png

    你也可以绘制Kaplan–Meier生存函数图、Kaplan–Meier失效函数图、Nelson-Aalen累积风险函数图即平滑危险估计图等,可在[统计>>生存分析>>回归模型>>图形] 或 [图形>>生存分析图] 中进行。

    【4】比例风险(PH)假定检验:比例风险假定是Cox回归的前提条件,可以通过计算检验,也可以通过图示法。

    4.1计算检验命令:整体检验estat phtest;单独检验每个协变量estat phtest, detail

    菜单可通过以下途径进入

    统计>>生存分析>>回归模型>>比例风险假设检验

    图形>>生存分析图>>stcox后比例风险假设检验

    统计>>后验估计,在打开的[后验估计选择器]中选择[比例风险假定的检验];

    点击[开始]进入对话框[estat-后验估计统计量],选择[基于Schoenfeld残差的比例风险假设检验(phtest)]。如想进一步检验每个协变量的比例风险假定,可在继续点击对话框[estat-后验估计统计量]中的[选项]按钮,在打开的对话框中选择复选框[单独检验每个协变量的比例风险假设]。

    结果显示P>0.05,满足比例风险假设。

    376affbbb25a82123ecf76151e50a1c2.png

    4.2 图示法可以通过生存函数的双对数图(Log-log plot of survival)Kaplan–MeierCox预测的生存曲线图(Kaplan–Meier and predicted survival plot)

    stphplot plots sln{−ln(survival)} curves for each category of a nominal or ordinal covariate versus ln(analysis time). These are often referred to as “log-log” plots. Optionally, these estimates can be adjusted for covariates. The proportional-hazards assumption is not violated when the curves are parallel.

    stcoxkm plots Kaplan–Meier observed survival curves and compares them with the Cox predicted curves for the same variable. The closer the observed values are to the predicted, the less likely it is that the proportional-hazards assumption has been violated.

    4.2.1 Log-log plot of survival可以拟合命令单独或者分层的Cox模型,并可根据调整变量进行调整。

    命令1:stphplot, by(recur)

    根据分类变量recur拟合单独的Cox模型,结果是下图左;

    命令2:stphplot, strata(recur) adjust(age) zero

    根据分类变量recur拟合分层的Cox模型,并根据变量age进行调整,调增方法是将age值调整为0(默认是均值),结果为下图右。

    菜单可通过以下途径进入stphplot-生存函数的双对数图(Log-log plot of survival)

    统计>>生存分析>>回归模型>>比例风险假设的图形评估

    图形>>生存分析图>>比例风险假设检验

    选入相应的自变量和分层变量即可。

    结果显示变量recur两个水平基本“平行”,满足比例风险的假定。

    29842a6f8403f614977c65d96f3f6019.png

    4.2.2 Kaplan–Meier and predicted survival plot命令:stcoxkm, by(recur)

    当然你可以单独绘制recur两个水平的观测预测曲线:stcoxkm, by(recur) separate

    菜单操作:

    统计>>生存分析>>回归模型>>Kaplan–Meier生存曲线和Cox预测曲线比较;

    图形>>生存分析图>>Kaplan–Meier和Cox生存曲线比较

    选入相应的自变量即可。

    结果显示Kaplan–Meier观测曲线与Cox回归预测曲线重合性较好,满足比例风险假设。

    7f32e014b6e60e4e36d8ca58b5626189.png

    如果运气不好,你的数据违背了比例风险的假定,可以考虑采用含时依协变量的Cox回归,而且含时依协变量的Cox回归也可以作为验证比例风险的假设的一种手段。个人理解,这个所谓时依协变量,其实就是在Cox回归里构建的一个交互作用项。这个关于含时依协变量的Cox回归我们放在以后分享。

    转自个人微信公众号【Memo_Cleon】的统计学习笔记:生存分析之Cox回归

    END

    展开全文
  • Logistic 回归常用于分析二分类因变量个自变量的关 系, 本文通过案例解析分类变量的 Logistic 回归, 借助于 SPSS 软件 实现 Logistic 回归过程, 并对分类因变量的 Logistic 回归做简单 介绍。
  • 以 “survival”包自带的晚期肺癌数据集 lung 为例,演示cox回归及可视化、模型评价与选择及变量筛选。 library(DescTools) library(survival) library(survminer) library(mice) #载入生存资料以survival包自带...
  • 变量的赋值和部分原始数据见表1和表2。 表1. 肺癌患者生存的影响因素与赋值 表2. 两组患者的生存情况 对数据结构的分析 该研究以死亡为结局,治疗方式为主要研究因素,每个研究对象都有生存时间(随访开始到...
  • Cox回归模型又称为比例风险回归模型,该模型以生存结局和生存时间作为因变量,进而分析众多因素对生存期的影响,是一个典型的因素分析方法。 SPSS中就带有Cox回归模型方法,本节将带大家进行深入的了解与探索,话...
  • 史上最详细Cox回归列线图制作教程

    万次阅读 2017-12-20 22:43:31
    列线图是一种直观有效地展示Cox回归结果的一种方法。最有价值的是进行结局的预测,同时可以通过直线的长度来表示不同变量对结局的影响,以及变量的不同取值对结局的影响。正如下图中所示。 举例来说,一个...
  • Cox比例风险回归模型单因素因素生存分析

    万次阅读 多人点赞 2020-03-13 12:00:20
    Cox比例风险回归模型单因素因素生存分析 欢迎使用Markdown编辑器 ...Cox比例风险回归模型,简称Cox回归模型。该模型又英国统计学家D.R.Cox于1972年提出,主要用于肿瘤和其他慢性病的预后分析,也可用于...
  • 除了这些传统的单因素分析方法,我们在阅读文献时也能常常遇到“单因素回归分析”这样的说法,例如我们在之前推送的《如何理解回归模型中的“调整”和“独立作用”》一文中所引用的研究实例。表1. 单因素(Unvariat.....
  • 上一篇文章写了数据分析系列:...这就是归因分析中的渠道归因。 对于渠道归因,有一些启发式的归因方法,比如“最终点击”(将订单归属于最后一个渠道)、“非最终点击”(归属于倒数第二个渠道)、“首次点击...
  • 用R语言进行Cox回归生存分析

    万次阅读 2019-06-29 20:13:00
    欢迎关注”生信修炼手册”!在生存分析中,探究生存时间的影响因素是一个重要的研究内容,通过KM和log-rank test检验的方法,只能够处理单个二分类因素的生存数据。当想探究个因素或...
  • 生存分析(Survival Analysis)、Cox风险比例回归模型(Cox proportional hazards model)及C-index 1. 生存分析 生存分析指的是一系列用来探究所感兴趣的事件的发生的时间的统计方法。常见的有1)癌症患者生存时间...
  • Cox回归+Cox比例风险回归模型

    千次阅读 2021-05-16 10:34:21
    Cox回归+Cox比例风险回归模型 COX回归模型,又称“比例风险回归模型(proportional hazards model,简称Cox模型)”,是由英国统计学家D.R.Cox(1972)年提出的一种半参数回归模型。该模型以生存结局和生存时间为应...
  • Cox风险回归分析

    千次阅读 2017-04-24 20:32:08
    3.log-rank只是对单因素进行了检验,建立回归模型,考虑因素的情况 (中途换了一组数据--恶性肿瘤的预后) 从结果中可以看出,X--(4,5,6)有显著性,RR值即为exp(coef),ceof即为β i值 4,library...
  • 2) Cox回归因素 3)随机森林因子(根据cox回归结果) 4)nomogram ( 根据cox回归结果,建立了中位生存时间,1年5年生存率的概率计算) -------------------????????‍♀️本文只有干货,非常干!????----------...
  • #glmnet包含有线性回归,逻辑回归,泊松计数模型,cox回归模型,分类逻辑回归响应线性回归 #阿法系数=0是岭回归,阿法系数=1,是lasso回归 ############################################################ data...
  • 实现图中描述的cox分析和lasso回归, 附代码和注释![图片说明](https://img-ask.csdn.net/upload/201909/07/1567833588_655859.jpg)
  • variable.list, in.variable = "NULL", data, sle = 0.15, sls = 0.15, vif.threshold = 999) Arguments Time The 'Time' (time to an event) for the sepcified Cox's proportional hazards model as in coxph()....
  • https://zhuanlan.zhihu.com/p/81424108?utm_source=wechat_session
  • 关于cox单因素与因素分析

    千次阅读 2021-07-14 09:04:23
    #写在前面 说来惭愧,跟师姐聊天的时候,我...当被问到为什么要做因素时,我想了一下说是为了去除共线性的问题,也就是a再做单因素分析时可能为显著,但是是其他变量带来的协同效应,所以再进入因素cox回归后,如果
  • https://blog.csdn.net/RicheyLee/article/details/51067339?utm_medium=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-3.channel_param&depth_1-utm_source=distribute.pc_relevant....
  • 在基因数据的分析中,经常会用到...接下来以线性回归为例介绍其在R语言中的实现,当然在logistic回归、cox回归也是可用lasso的。实例数据data("longley")R包(glmnet)library(glmnet)我们用交叉验证来确定lamda的值...
  • SPSS详细操作:生存资料的Cox回归分析 一、问题与数据 某研究者拟观察某新药的抗肿瘤效果,将70名肺癌患者随机分为两组,分别采用该新药和常规药物进行治疗,观察两组肺癌患者的生存情况,共随访2年。研究以死亡为...
  • SPSS教程之生存分析Cox回归模型(比例风险模型)

    万次阅读 多人点赞 2016-04-05 19:22:27
    最近有同学问师兄,“最近我要做生存分析,可是我不太会,也不太懂,师兄能不能教教我”,好吧,今天开一贴,讲讲这个。有同样的问题的同学可以一起来看看,毕竟在临床、科研上,这方面知识还是很受用的。有什么想跟...
  • ![图片说明](https://img-ask.csdn.net/upload/201909/07/1567833184_85101.jpg) 用python如何实现下图中要求的cox生存分析? 最好附代码和注释
  • 【临床研究】---多元回归分析中的变量筛选问题方法选择的思考路径:1、变量筛选方法的归纳1)变量筛选的一般流程:①逐个变量:单因素回归分析②分析P值:依据样本量大小情况调整P值选择范围③纳入规则:将单因素...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 2,064
精华内容 825
关键字:

多变量cox回归分析