临床预测模型概述6-统计模型实操-单/多因素Cox回归
原创
基础知识回顾：
https://mp.weixin.qq.com/s/pXRZ1rYUr3lwH5OlDeB0_Q
https://mp.weixin.qq.com/s/UVR6ZHCwhWqTfFBmPYPV9Q
接下来我们进行cox回归模型的实际操练。
简单回顾一下cox回归，在各种临床/基础数据分析中，经常需要分析各种影响/危险因素对疾病/状态随着时间变化而产生的影响作用，如研究肝癌患者的生存或死亡风险如何受到不同治疗方式、年龄、饮食习惯、饮酒和抽烟等因素的影响，并探索这些因素随时间变化的作用机制。再简单的说就是，不同影响因素，对肝癌患者发生死亡事件在一段时间上发生的概率的影响作用。
Cox回归
1、首先输入数据，数据内容形式可见logistic回归推文
rm(list = ls())
load("TCGA-LIHC_sur_model.Rdata")
head(exprSet)[1:4,1:4]
head(meta)2、check数据
meta$OS <- factor(meta$OS,levels = c("0","1"))
meta$race <- factor(meta$race,levels = c("ASIAN","AMERICAN INDIAN OR ALASKA NATIVE","BLACK OR AFRICAN AMERICAN","WHITE"))
meta$stage <- factor(meta$stage,levels = c("I","II","III","IV"))
meta$T <- factor(meta$T,levels = c("1","2","3","4"))
meta$age <- ifelse(meta$age > 60,">60","<=60")
meta$age <- factor(meta$age,levels = c("<=60",">60"))
str(meta)同样的这里对于有一些参数并没有进行因子化，以及存在NA值，这些情况会在下边的探索中解释。
3、构建模型
# 使用survival包中的coxph函数，这是cox回归中最常用的函数之一
library(survival)
formula <- as.formula(Surv(OS.time,OS==1)~race+age+gender+stage+T+N+M)
# Warning message:
# In coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
# Loglik converged before variable 1,6,8,9,13 ; coefficient may be infinite.
fit <- coxph(formula,data = meta)
# Warning message:
# In coxph.fit(X, Y, istrat, offset, init, control, weights = weights, :
# Loglik converged before variable 1,6,8,9,13 ; coefficient may be infinite. 
这个警告 "Loglik converged before variable 1,6,8,9,13; coefficient may be infinite." 在 Cox 比例风险模型的拟合过程中出现，它表明在模型收敛之前，某些变量（编号为 1, 6, 8, 9, 13，可以对照summary结果看一下）的系数估计可能趋于无穷大。这通常是由于以下几种原因引起的：
- 完全分离（Complete Separation） 如果数据集中某些变量可以完美预测事件的发生(这种情况当然是不可能的啦)，例如某个预测变量在所有事件发生的情况下都是一个值，在没有事件发生的情况下都是另一个值，这将导致该变量的系数估计趋向无穷大。这种情况在生存分析中尤为常见，尤其是在处理较小的数据集或有强烈预测因子的情况下。
- 多重共线性（Multicollinearity） 如果模型中包括的变量间存在高度共线性，即一些变量彼此非常相关，这也可能导致某些变量的系数估计非常不稳定，甚至趋向无限大。共线性问题会使得模型参数的估计变得困难，因为模型难以区分这些相关变量的独立效应。
- 数据稀疏（Sparse Data Issue） 在某些分类变量中，如果某些类别的观测非常少，也可能导致系数估计的不稳定性。当包含的分类变量的某些水平在数据中几乎或完全没有观测时，这种情况可能会发生。
解决这些问题的策略包括：
● 数据检查和预处理：检查数据中是否存在数据分离现象，确认各变量之间的关系，可能需要重新编码或合并一些变量的类别。
● 模型简化：尝试减少模型中的变量，特别是那些引起警告的变量。检查这些变量的必要性和影响，考虑从模型中移除或替换它们。
● 使用正则化技术：考虑使用像 Lasso 或 Ridge 这样的正则化方法，这些方法通过在估计过程中添加惩罚项，可以帮助缓解共线性问题并提高模型的稳定性。
个人经验，一般出现这种情况最主要是数据的问题，需要进行数据检查和更好的预处理！
summary(fit)
# Call:
# coxph(formula = formula, data = meta)
#
# n= 227, number of events= 75
# (138 observations deleted due to missingness)
#
# coef exp(coef) se(coef) z Pr(>|z|)
# raceAMERICAN INDIAN OR ALASKA NATIVE -1.538e+01 2.084e-07 4.528e+03 -0.003 0.997
# raceBLACK OR AFRICAN AMERICAN 9.079e-01 2.479e+00 7.497e-01 1.211 0.226
# raceWHITE -3.791e-01 6.845e-01 2.972e-01 -1.276 0.202
# age>60 2.243e-01 1.251e+00 2.523e-01 0.889 0.374
# genderMALE -2.288e-01 7.955e-01 2.739e-01 -0.835 0.404
# stageII 1.692e+01 2.228e+07 3.898e+03 0.004 0.997
# stageIII 1.176e-01 1.125e+00 1.493e+00 0.079 0.937
# stageIV -1.613e+01 9.929e-08 1.876e+04 -0.001 0.999
# T2 -1.654e+01 6.580e-08 3.898e+03 -0.004 0.997
# T3 1.039e+00 2.827e+00 1.461e+00 0.711 0.477
# T4 1.586e+00 4.885e+00 1.540e+00 1.030 0.303
# N 1.327e+00 3.772e+00 1.035e+00 1.283 0.199
# M 1.649e+01 1.448e+07 1.876e+04 0.001 0.999
#
# exp(coef) exp(-coef) lower .95 upper .95
# raceAMERICAN INDIAN OR ALASKA NATIVE 2.084e-07 4.800e+06 0.00000 Inf
# raceBLACK OR AFRICAN AMERICAN 2.479e+00 4.034e-01 0.57032 10.777
# raceWHITE 6.845e-01 1.461e+00 0.38230 1.226
# age>60 1.251e+00 7.991e-01 0.76326 2.052
# genderMALE 7.955e-01 1.257e+00 0.46509 1.361
# stageII 2.228e+07 4.488e-08 0.00000 Inf
# stageIII 1.125e+00 8.890e-01 0.06026 20.996
# stageIV 9.929e-08 1.007e+07 0.00000 Inf
# T2 6.580e-08 1.520e+07 0.00000 Inf
# T3 2.827e+00 3.538e-01 0.16122 49.561
# T4 4.885e+00 2.047e-01 0.23889 99.912
# N 3.772e+00 2.651e-01 0.49652 28.649
# M 1.448e+07 6.907e-08 0.00000 Inf
#
# Concordance= 0.655 (se = 0.033 )
# Likelihood ratio test= 29.29 on 13 df, p=0.006
# Wald test = 29.84 on 13 df, p=0.005
# Score (logrank) test = 36.03 on 13 df, p=6e-04接下来详细解析一下summary(fit)得到的结果
1）常规参数
● 系数（coef）: 回归系数反映了每个协变量对生存时间风险比的影响。负系数表明该变量与生存时间正相关（降低风险），正系数表明与生存时间负相关（增加风险）。
● 指数形式的系数（exp(coef): 这就是风险比（hazard ratio），表明了相应变量对生存风险的倍数影响。例如，raceBLACK OR AFRICAN AMERICAN 的风险比为2.479，意味着与参考类别相比，该族群的死亡风险高2.479倍。
● 标准误差（se(coef): 系数的标准误差，衡量估计的精确度。
● z 值和 p 值: 用于测试每个系数的统计显著性。Z 值是通过将系数估计值（β）除以其标准误差（SE）来计算的。一旦计算出 Z 值，接下来就可以计算 P 值。P 值是通过查找标准正态分布表或使用统计软件中的相关函数来得出的。
● lower .95 upper .95：exp(coef)的95%置信区间，可信区间越窄，可信度越高。
2）重点关注
● 极大的系数和标准误差: 有些变量如 raceAMERICAN INDIAN OR ALASKA NATIVE, stageII, stageIV, T2, 和 M 显示出极大的系数和/或标准误差，这通常是模型不稳定的迹象，可能是因为样本中这些类别的数据非常稀少或者完全分离的情况。这种情况下是必须要重新检查数据的。
● 风险比的极端值: 如 stageII 的风险比达到了约2200万，这样的极端值通常不是实际情况，反映了数据问题或模型拟合问题。这种情况下是必须要重新检查数据的。
● 模型评估指标: Concordance = 0.655这个就是C-index值: 表示模型对生存数据排序的能力为65.5%，高于随机模型（50%）。似然比检验、Wald检验和Score(log rank)检验: 提供了模型整体拟合优度的统计显著性检验，p 值都显示模型是统计显著的。
4、提取模型的关键信息
# β/Estimate值
coef(fit) #或者 coefficients(fit)
# β值的95%可信区间
confint(fit)
# HR值
exp(coef(fit))
# OR值的95%可信区间
exp(confint(fit))[1:4,1:2]
# P值
summary(fit)$coefficients[,5]具体不展示了
5、模型数据提取
# 单因素cox数据手动提取
# 构建自定义函数
uni_cox<-function(x){
fml<-as.formula(paste0("Surv(OS.time,OS==1)~",x))
fit1<-coxph(fml,data = meta)
fit2<-summary(fit1)
#计算我们所要的指标
HR<-round(fit2$conf.int[,1],2)
CI<-paste0(round(fit2$conf.int[,3:4],2),collapse = '-')
P<-round(fit2$coefficients[,5],3)
#将计算出来的指标制作为数据框
uni_cox<-data.frame('characteristics'=x,'HR'=HR,'CI'=CI,'P'=P)
return(uni_cox)
}
#指定参与分析的自变量
variable<-colnames(meta)[c(4:10)]
variable
#运行上面自定义批量执行函数
uni_cox_res<-lapply(variable,uni_cox)
uni_cox_res
library(plyr)
#生成单变量分析的结果
uni_cox_res<-ldply(uni_cox_res,data.frame)
uni_cox_res
# characteristics HR CI P
# 1 race 0.00 0-0.81-0.89-Inf-4-1.9 0.995
# 2 race 1.80 0-0.81-0.89-Inf-4-1.9 0.150
# 3 race 1.30 0-0.81-0.89-Inf-4-1.9 0.170
# 4 age 1.27 0.89-1.79 0.183
# 5 gender 0.82 0.58-1.17 0.284
# 6 stage 1.39 0.85-1.71-1.65-2.26-3.97-17.32 0.187
# 7 stage 2.61 0.85-1.71-1.65-2.26-3.97-17.32 0.000
# 8 stage 5.35 0.85-1.71-1.65-2.26-3.97-17.32 0.005
# 9 T 1.40 0.89-1.69-2.57-2.22-3.87-10.29 0.147
# 10 T 2.56 0.89-1.69-2.57-2.22-3.87-10.29 0.000
# 11 T 5.15 0.89-1.69-2.57-2.22-3.87-10.29 0.000
# 12 N 1.98 0.48-8.07 0.342
# 13 M 3.91 1.23-12.44 0.021
# broom方式
library(broom)
formula <- as.formula(Surv(OS.time,OS==1)~race+age+gender+stage+T+N+M)
fit <- coxph(formula,data=meta)
tidy(fit,exponentiate = T, conf.int = T) # 将模型结果转化为整洁的数据框
# # A tibble: 13 × 7
# term estimate std.error statistic p.value conf.low conf.high
# <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 raceAMERICAN INDIAN OR ALASKA NATIVE 2.08e-7 4528. -0.00340 0.997 0 Inf
# 2 raceBLACK OR AFRICAN AMERICAN 2.48e+0 0.750 1.21 0.226 0.570 10.8
# 3 raceWHITE 6.84e-1 0.297 -1.28 0.202 0.382 1.23
# 4 age>60 1.25e+0 0.252 0.889 0.374 0.763 2.05
# 5 genderMALE 7.96e-1 0.274 -0.835 0.404 0.465 1.36
# 6 stageII 2.23e+7 3898. 0.00434 0.997 0 Inf
# 7 stageIII 1.12e+0 1.49 0.0788 0.937 0.0603 21.0
# 8 stageIV 9.93e-8 18761. -0.000860 0.999 0 Inf
# 9 T2 6.58e-8 3898. -0.00424 0.997 0 Inf
# 10 T3 2.83e+0 1.46 0.711 0.477 0.161 49.6
# 11 T4 4.89e+0 1.54 1.03 0.303 0.239 99.9
# 12 N 3.77e+0 1.03 1.28 0.199 0.497 28.6
# 13 M 1.45e+7 18761. 0.000879 0.999 0 Inf 手动提取和broom函数两种方法
6、多因素cox回归筛选自变量
# 在用forward/backward/both方法的时候需要去除NA值！
meta <- na.omit(meta)
dim(meta) #去除了蛮多数据的
# [1] 227 10
#1.多因素enter回归
formula <- as.formula(Surv(OS.time,OS==1)~race+age+gender+stage+T+N+M)
model_enter <- coxph(formula,data = meta)
model_enter
summary(model_enter)
# 构建最小模型
start_model <- coxph(Surv(OS.time,OS==1) ~ 1, data = meta) # 仅包含截距
library(MASS)
# 多因素-backward
model_backward <- stepAIC(model_enter, direction = "backward")
summary(model_backward)
# Call:
# coxph(formula = Surv(OS.time, OS == 1) ~ T, data = meta)
#
# n= 227, number of events= 75
#
# coef exp(coef) se(coef) z Pr(>|z|)
# T2 0.3538 1.4245 0.3273 1.081 0.279714
# T3 1.0555 2.8734 0.2766 3.815 0.000136 ***
# T4 1.6914 5.4270 0.4344 3.894 9.87e-05 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# exp(coef) exp(-coef) lower .95 upper .95
# T2 1.424 0.7020 0.750 2.706
# T3 2.873 0.3480 1.671 4.942
# T4 5.427 0.1843 2.316 12.715
#
# Concordance= 0.632 (se = 0.033 )
# Likelihood ratio test= 21.53 on 3 df, p=8e-05
# Wald test = 23.42 on 3 df, p=3e-05
# Score (logrank) test = 26.5 on 3 df, p=7e-06
# 多因素-forward
model_foward <- stepAIC(start_model, scope = list(lower = start_model, upper = model_enter), direction = "forward")
summary(model_forward)
# Call:
# coxph(formula = Surv(OS.time, OS == 1) ~ race + age + gender +
# stage + T + N + M, data = meta)
#
# n= 227, number of events= 75
#
# coef exp(coef) se(coef) z Pr(>|z|)
# raceAMERICAN INDIAN OR ALASKA NATIVE -1.538e+01 2.084e-07 4.528e+03 -0.003 0.997
# raceBLACK OR AFRICAN AMERICAN 9.079e-01 2.479e+00 7.497e-01 1.211 0.226
# raceWHITE -3.791e-01 6.845e-01 2.972e-01 -1.276 0.202
# age>60 2.243e-01 1.251e+00 2.523e-01 0.889 0.374
# genderMALE -2.288e-01 7.955e-01 2.739e-01 -0.835 0.404
# stageII 1.692e+01 2.228e+07 3.898e+03 0.004 0.997
# stageIII 1.176e-01 1.125e+00 1.493e+00 0.079 0.937
# stageIV -1.613e+01 9.929e-08 1.876e+04 -0.001 0.999
# T2 -1.654e+01 6.580e-08 3.898e+03 -0.004 0.997
# T3 1.039e+00 2.827e+00 1.461e+00 0.711 0.477
# T4 1.586e+00 4.885e+00 1.540e+00 1.030 0.303
# N 1.327e+00 3.772e+00 1.035e+00 1.283 0.199
# M 1.649e+01 1.448e+07 1.876e+04 0.001 0.999
#
# exp(coef) exp(-coef) lower .95 upper .95
# raceAMERICAN INDIAN OR ALASKA NATIVE 2.084e-07 4.800e+06 0.00000 Inf
# raceBLACK OR AFRICAN AMERICAN 2.479e+00 4.034e-01 0.57032 10.777
# raceWHITE 6.845e-01 1.461e+00 0.38230 1.226
# age>60 1.251e+00 7.991e-01 0.76326 2.052
# genderMALE 7.955e-01 1.257e+00 0.46509 1.361
# stageII 2.228e+07 4.488e-08 0.00000 Inf
# stageIII 1.125e+00 8.890e-01 0.06026 20.996
# stageIV 9.929e-08 1.007e+07 0.00000 Inf
# T2 6.580e-08 1.520e+07 0.00000 Inf
# T3 2.827e+00 3.538e-01 0.16122 49.561
# T4 4.885e+00 2.047e-01 0.23889 99.912
# N 3.772e+00 2.651e-01 0.49652 28.649
# M 1.448e+07 6.907e-08 0.00000 Inf
#
# Concordance= 0.655 (se = 0.033 )
# Likelihood ratio test= 29.29 on 13 df, p=0.006
# Wald test = 29.84 on 13 df, p=0.005
# Score (logrank) test = 36.03 on 13 df, p=6e-04
# 多因素-both
model_both <- stepAIC(start_model, scope = list(lower = start_model, upper = model_enter), direction = "both")
summary(model_both)
# Start: AIC=719.94
# Surv(OS.time, OS == 1) ~ 1
#
# Df AIC
# + T 3 704.40
# + stage 3 704.87
# + M 1 718.31
# <none> 719.94
# + gender 1 720.49
# + N 1 721.11
# + age 1 721.18
# + race 3 724.24
#
# Step: AIC=704.4
# Surv(OS.time, OS == 1) ~ T
#
# Df AIC
# <none> 704.40
# + N 1 705.36
# + age 1 705.69
# + gender 1 706.09
# + M 1 706.23
# + race 3 707.74
# + stage 3 708.30
# - T 3 719.94
# 计算多个模型的 AIC 值
aic_values <- AIC(model_enter, model_forward, model_backward, model_both)
# 打印 AIC 值
print(aic_values)
# 找到 AIC 值最小的模型
best_model <- aic_values[which.min(aic_values$AIC), ]
print(best_model)
# 模型AIC比较
# 选择AIC值最小的模型
anova(model_enter,model_forward,test = "Chisq")
anova(model_enter,model_backward,test = "Chisq")
anova(model_enter,model_both,test = "Chisq")
anova(model_forward,model_backward,test = "Chisq")
anova(model_forward,model_both,test = "Chisq")
anova(model_backward,model_both,test = "Chisq")
# 使用MuMIn包
library(MuMIn)
# 将所有模型放入一个列表
model_list <- list(Enter = model_enter, Forward = model_forward, Backward = model_backward, Both = model_both)
# 使用 MuMIn 包的 model.sel() 函数进行比较
model_comparison <- model.sel(model_list)
# 打印模型比较结果
print(model_comparison)
# 选择 AICc 值最小的模型
best_model <- get.models(model_comparison, 1)[[1]]
summary(best_model)参考资料：
1、coxph：https://rdrr.io/cran/survival/man/coxph.html
2、MASS： https://cran.r-project.org/web/packages/MASS/MASS.pdf
3、broom：https://github.com/tidymodels/broom
4、MuMIn：https://mumin-dataset.github.io/
注：若对内容有疑惑或者有发现明确错误的朋友，请联系后台(欢迎交流)。更多内容可关注公众号：生信方舟
- END -
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
评论
登录后参与评论
推荐阅读
目录
腾讯云开发者
