临床预测模型/机器学习-Coxboost算法学习

CoxBoost 是一种用于生存分析的统计和机器学习方法，特别适合处理高维数据（例如基因组数据）中的 Cox 回归模型。它将 Cox 比例风险模型（Cox Proportional Hazards Model）与 Boosting（提升）算法结合，用来在大量特征（变量）和相对较少样本的数据集中进行生存时间预测和变量筛选。

这里需要提一下，高维数据的概念，通常来说特征(变量)远大于观测数量(样本)的就是高维数据(比如：测序数据)。

目前笔者所了解到的高低维度数据的界定没有严格的标准，如果变量数不多比如在10-100之间，样本量又很多远超变量数那就可以看做是低维数据(如果这里的概念有明确错误的话请尽管批评指正)。

低维度数据通常使用常规的算法构建的模型就可具有良好的可解释性，而高维度数据是可以用常规算法去进行分析，但通常建议使用更加“强大”的工具，比如Coxboost。

稍微读一读Coxboost的文献(PMID: 19244389IF: 4.4 Q1 B3)

通常而言如果协变量对发生事件的特定风险有影响的时候，那么忽略协变量的方式其实是不可取的。有一种方法是构建多个COX比例风险模型并结合解释(保护/危害因素)，但在低维度数据中，通过累积发生函数（Cumulative Incidence Function, CIF）对两个特定原因的模型结合进行综合解释已经是有难度了，那么在高维度数据中通过这种方法去解释就更加麻烦(可以通过CIF自行探索)。Fine 和 Gray（1999）提出了一个适用于竞争风险情境的亚分布风险模型（Subdistribution Hazard Model），通过对CIF建模，来描述特定事件在有竞争风险时的累积发生概率。该模型中的亚分布风险与 CIF 直接相关，因此只需为感兴趣的事件拟合一个模型即可。相比之下，传统的特定原因风险模型（Cause-specific Hazard Model）与 CIF 没有直接关系，无法准确反映在竞争风险下的累积发生概率。

累积发生函数（CIF）是什么？

1. 累积发生函数表示某一特定事件（例如癌症复发）在一段时间内发生的累积概率，考虑到了竞争事件（例如非癌症原因死亡）的存在。
2. 由于竞争事件的存在，每个患者只能经历一个事件，因此 CIF 计算的是“在某一时间之前特定事件发生的概率”，而不是单纯地累计特定原因的发生次数。

亚分布风险（Subdistribution Hazard）与 CIF 的关系

1. Fine 和 Gray 的亚分布风险模型直接针对 CIF 建模，描述的是在考虑竞争风险的情况下特定事件的“累积发生率”。
2. 这个模型可以解释为，假如我们从一开始就“坚持观察”这个特定事件，即使有些个体因其他事件（例如非癌症死亡）退出了观察，我们仍然通过调整模型，将这些退出的人纳入总体，使得最终得到的是“在竞争风险条件下该事件的累积发生概率”。
3. 因此，亚分布风险模型能够直接反映特定事件的CIF，它既考虑了感兴趣事件的发生，也考虑了竞争事件的影响。

特定原因风险（Cause-specific Hazard）与 CIF 的关系

1. 特定原因风险模型关注的是特定事件在“没有竞争事件的干扰下”的瞬时发生率，即在每个时间点该事件发生的风险率。
2. 这种模型是“条件风险”，只考虑该事件的发生情况，而一旦发生了竞争事件（例如非癌症死亡），个体就会被视为“删失”数据（即不再纳入分析），而不会继续影响模型。
3. 因此，特定原因风险模型并不能直接描述 CIF，因为它没有考虑到竞争事件的发生情况对特定事件累积概率的影响。开发者使用一种改进的Boosting方法来适应高维数据的比例亚分布风险模型，这种方法基于Fine和Gray醍醐的加权似然，能够更好的识别出具有价值的变量。预测误差曲线（即真实与预测累积发生函数之间的均方差）在单一事件类型的时间-事件分析中已被广泛应用，通过 bootstrap .632+ 技术可以估计预测误差，而无需预留测试集。

CoxBoost中的重要知识/概念

Cox 回归模型

1. Cox 回归模型是一种广泛应用于生存分析的模型，用来评估多种变量对生存时间的影响。它基于比例风险假设，即假设不同个体的风险比随时间保持不变。
2. CoxBoost 使用 Cox 回归模型来处理生存数据，但相比于传统 Cox 回归，它通过 Boosting 提高了对高维数据的适应性。

Boosting 算法

1. Boosting 是一种集成学习方法，通过多次迭代来增强弱学习器的性能。
2. 在 CoxBoost 中，Boosting 的作用是通过逐步添加重要的变量并优化系数，使模型能够在包含大量变量的生存数据中找到对生存时间最有影响的特征。

高维数据适用性

1. CoxBoost 特别适用于高维数据集，尤其是在样本数量少于变量数量的场景下，比如基因表达数据、转录组数据等 -omics 数据。这类数据集常常具有成千上万个特征，而传统 Cox 模型很难有效处理。
2. 通过 Boosting，CoxBoost 可以自动筛选出对生存时间有显著影响的变量，从而减少维度，提高模型的预测能力。

CoxBoost 的主要功能

1. 变量筛选：CoxBoost 能够在高维数据中自动选择与生存时间密切相关的变量，这一功能对基因组数据分析尤其有用，因为它可以从成千上万个基因中识别出最关键的基因。
2. 高维生存分析：它能处理变量数量远大于样本数量的情况，使得 CoxBoost 成为高维生存分析的一种重要工具。
3. 避免过拟合：通过 Boosting 和正则化，CoxBoost 控制了模型的复杂度，帮助避免高维数据中的过拟合问题。

分析步骤

1.导入

rm(list = ls())
library(CoxBoost)
load("data.Rdata")

2.数据预处理

# 把基因数据转置之后跟生存信息整合在一起
# 行为样本，列为生存信息+变量
meta <- meta[,c(1:3)]
head(meta)
#                                ID OS  OS.time
# TCGA-CR-7374-01A TCGA-CR-7374-01A  0 1.000000
# TCGA-CV-A45V-01A TCGA-CV-A45V-01A  1 1.066667
# TCGA-CV-7102-01A TCGA-CV-7102-01A  1 1.866667
# TCGA-MT-A67D-01A TCGA-MT-A67D-01A  0 1.866667
# TCGA-P3-A6T4-01A TCGA-P3-A6T4-01A  1 2.066667
# TCGA-CV-7255-01A TCGA-CV-7255-01A  1 2.133333
identical(rownames(meta),colnames(exprSet))
# [1] TRUE
meta <- cbind(meta,t(exprSet))
head(meta)[1:5,1:5]
#                                ID OS  OS.time    WASH7P AL627309.6
# TCGA-CR-7374-01A TCGA-CR-7374-01A  0 1.000000 0.5808846   3.117962
# TCGA-CV-A45V-01A TCGA-CV-A45V-01A  1 1.066667 1.4177642   6.250413
# TCGA-CV-7102-01A TCGA-CV-7102-01A  1 1.866667 0.6501330   1.219729
# TCGA-MT-A67D-01A TCGA-MT-A67D-01A  0 1.866667 1.2045780   3.038835
# TCGA-P3-A6T4-01A TCGA-P3-A6T4-01A  1 2.066667 1.3470145   3.799571
str(meta)
# 'data.frame': 493 obs. of  18238 variables:
#  $ ID                       : chr  "TCGA-CR-7374-01A" "TCGA-CV-A45V-01A" "TCGA-CV-7102-01A" "TCGA-MT-A67D-01A" ...
#  $ OS                       : int  0 1 1 0 1 1 1 1 1 0 ...
#  $ OS.time                  : num  1 1.07 1.87 1.87 2.07 ...
#  $ WASH7P                   : num  0.581 1.418 0.65 1.205 1.347 ...
#  $ AL627309.6               : num  3.12 6.25 1.22 3.04 3.8 ...
#  $ AL627309.7               : num  3.73 6.38 2.13 2.96 4.43 ...


# 数据分割 7:3,8:2 均可
# 划分是随机的，设置种子数可以让结果复现
set.seed(123)
ind <- sample(1:nrow(meta), size = 0.7*nrow(meta))
train <- meta[ind,]
test <- meta[-ind, ]

3.Coxboost分析

在 CoxBoost 模型中，penalty 和 stepno 是控制模型复杂性、选择变量、避免过拟合的两个重要参数，尤其是在高维数据中（即变量数量远多于样本数量的情况）。以下是它们在模型中的重要性和作用：

1. penalty（惩罚参数）的重要性

penalty 参数决定了模型在每次 Boosting 迭代时对变量系数更新的强度，即对每个变量进行惩罚的力度。惩罚越大，模型会更倾向于将变量的系数缩小到零，从而使模型更加稀疏。

1. 控制模型的稀疏性：较大的 penalty 值会使模型的更新更保守，许多变量的系数会被缩小到零。这在高维数据中尤其重要，因为模型在保留少量关键变量的同时，可以剔除许多不重要的变量，从而减少噪声。
2. 变量选择：通过设置适当的 penalty 值，CoxBoost 能够自动选择与生存时间最相关的变量，而不相关或影响较小的变量会因高惩罚而被排除（即系数变为零）。
3. 避免过拟合：在高维数据中，过拟合是常见问题。penalty 有助于控制模型的复杂度，从而避免模型过拟合训练数据。
4. 提高模型的解释性：高惩罚使得模型稀疏化，只包含少量重要变量。相比包含大量变量的复杂模型，稀疏模型的解释性更强，更容易解读哪些变量对生存时间有显著影响。

2. stepno（Boosting步数）的重要性

stepno 决定了模型迭代的步数，即 Boosting 过程中的更新次数。Boosting 的每一步都会选择一个变量并更新其系数，随着步数的增加，模型会逐渐拟合训练数据中的模式。

1. 控制模型的拟合程度：stepno 决定了模型对数据拟合的精细程度。较小的步数会产生一个较为简单的模型，而较大的步数会使模型更复杂。合适的步数可以让模型很好地捕捉数据的规律，但过多的步数则可能导致过拟合。
2. 防止过拟合：如果 stepno 设置得太高，模型可能会过度拟合训练数据，尤其是在高维数据中，过多的 Boosting 步数会导致模型对训练数据的噪声进行拟合。通过选择合适的 stepno，可以控制模型的复杂性，从而提高模型的泛化能力。
3. 模型复杂性的主控参数：在 CoxBoost 中，stepno 是主要的模型复杂性控制参数。适当的 stepno 设置能够保证模型足够复杂以拟合数据的规律，但不会过度复杂而导致不必要的变量或噪声进入模型。

penalty 和 stepno 的相互作用

1. 稀疏性与拟合之间的平衡：penalty 越大，模型越稀疏；而 stepno 越高，模型对数据的拟合程度越高。适当的 penalty 可以在一定程度上抵消高 stepno 可能引起的过拟合。
2. 变量选择与迭代步数的平衡：通过合理地设置 penalty 和 stepno，可以在变量选择和模型拟合之间取得平衡，使得最终模型既包含了关键变量，也不会对数据的噪声进行过度拟合。

确定最佳penalty值

# 输入矩阵哦~ 
dat <- as.matrix(train[,-c(1:3)])
penalty <- optimCoxBoostPenalty(time = train$OS.time, 
                                status = train$OS,                            
                                x = dat, 
                                trace = T,                               
                                parallel = T)
penalty$penalty
# [1] 2556

确定最佳stepno值

# 输入矩阵哦~ 
# dat <- as.matrix(train[,-c(1:3)])
stepno <- cv.CoxBoost(time = train$OS.time, 
                      status = train$OS,                            
                      x = dat,                     
                      maxstepno = 500,                     
                      K = 10,# 交叉验证的折数                    
                      type = "verweij", #  Verweij 和 van Houwelingen 方法     
                      penalty = penalty$penalty,
                      multicore = 4
                      )
stepno$optimal.step
#[1] 79

构建CoxBoost模型

set.seed(123)
fit <- CoxBoost(time = train$OS.time, 
                status = train$OS, 
                x = dat, # 要确认一下输入的矩阵情况哦
                stepno = stepno$optimal.step, #重采样次数
                #maxstepno = 200, # Boosting最大步数
                #multicore = T, # 多核运行
                #stratnotinfocus = 0, #等于0的样本不会被用于模型拟合
                penalty = penalty$penalty #惩罚参数用于控制模型复杂度
                #criterion = "hscore",
                #unpen.index = NULL # 所有变量都使用惩罚
                )

# summary(fit)提取关键变量
summary(fit)
# 79 boosting steps resulting in 29 non-zero coefficients  
# partial log-likelihood: -662.5563 
# 
# parameter estimates > 0:
#  TMCO1, TBCK, SNX14, RPL12, TPP1, NAT10, AC090589.1, TMEM223, EMC4, CCDC32, RSL24D1, SEC11A, DNAJA2, LINC02128, NOB1, NUDT7 
# parameter estimates < 0:
#  PIK3C2B, CELSR3, IMPG2, SPINK6, CPNE5, CUL9, FGD3, MS4A1, CATSPERB, PITPNM3, ZNF266, NIBAN3, AC118344.4 

# plot
pdf("CoxBoost.pdf",width = 9,height = 7)
plot(fit)
dev.off()

Boosting步数和非零系数数量

1. "79 boosting steps resulting in 29 non-zero coefficients" 表示在进行了 79 次 Boosting 步骤后，模型中有 29 个变量的系数不为零。这意味着在逐步变量选择过程中，这 29 个变量被认为对生存时间（或风险）具有显著影响，其余变量的系数被缩小到零，表明它们对模型贡献不大，被排除在模型之外。
2. 这种稀疏化的结果有助于在众多变量中筛选出关键的特征，从而提高模型的解释性和预测性能。

部分对数似然值（Partial Log-Likelihood）

1. "partial log-likelihood: -662.5563" 表示模型的部分对数似然值。对数似然值是生存模型拟合优度的衡量指标之一，对数似然值越高，说明模型与数据的拟合效果越好。
2. 此值可以用于与其他模型进行比较，选择拟合效果最好的模型。

正参数估计（Parameter Estimates > 0）

1. 这些变量的系数为正值，即随着这些变量值的增加，事件发生的风险也增加。换句话说，这些基因或特征变量对生存期有不利影响。
2. 正系数变量：TMCO1, TBCK, SNX14, RPL12, TPP1, NAT10, AC090589.1, TMEM223, EMC4, CCDC32, RSL24D1, SEC11A, DNAJA2, LINC02128, NOB1, NUDT7

负参数估计（Parameter Estimates < 0）

1. 这些变量的系数为负值，即随着这些变量值的增加，事件发生的风险降低。换句话说，这些基因或特征变量对生存期有保护作用。
2. 负系数变量：PIK3C2B, CELSR3, IMPG2, SPINK6, CPNE5, CUL9, FGD3, MS4A1, CATSPERB, PITPNM3, ZNF266, NIBAN3, AC118344.4这里就可以提取summary(fit)中不为0的变量后续再进行其他算法的联合分析~ 如果不进行后续的分析，也可以直接预测风险分数。

4.验证集分析及预测

# 输入矩阵哦~ 
data <- as.matrix(test[,-c(1:3)])
# predict
verify <- predict(fit, 
                  newdata = data, 
                  newtime = test$OS.time, 
                  newstatus = test$OS,
                  type = "logplik")  # lp
verify
# [1] -306.3869

pred.train <- predict(fit, newdata = dat) # lp
head(pred.train)[,1:5]
# [1] -0.90537999 -0.07241175  0.74893233  0.52619319 -0.36423249

pred.test <- predict(fit, newdata = data) # lp
head(pred.test)[,1:5]
# [1] -0.19091011  0.32812666  0.32118698  0.13696966 -0.03886175

type = "logplik"

1. "logplik" 表示计算 对数似然值（log-likelihood）。
2. 当 type = "logplik" 时，predict 函数会计算模型的部分对数似然值，用于衡量模型在给定数据集上的拟合优度。这个值通常用于评估模型的性能，特别是在不同步数的 Boosting 过程中，观察对数似然值的变化可以帮助选择最佳的 Boosting 步数。

type = "lp"

1. "lp" 表示计算 线性预测值（linear predictor）。
2. 当 type = "lp" 时，predict 函数会返回模型对新数据的线性预测值，即风险评分。这个风险评分是根据 Cox 模型的线性组合（线性预测器）计算得出的，每个样本的风险分数是基于其特征变量和模型系数计算的。

5.提取风险分数

# 这里丢失了样本信息
# 笔者结合了各大up主的内容来看，这里是跟样本顺序一致的。
RS_train = t(pred.train)
head(RS_train)
#             [,1]
# [1,] -0.90537999
# [2,] -0.07241175
# [3,]  0.74893233
# [4,]  0.52619319
# [5,] -0.36423249
# [6,] -0.21650543

RS_test = t(pred.test)
head(RS_test)
#             [,1]
# [1,] -0.19091011
# [2,]  0.32812666
# [3,]  0.32118698
# [4,]  0.13696966
# [5,] -0.03886175
# [6,]  0.03699118

得到了风险分数后续就很好办了~

参考资料：

1. github: https://github.com/binderh/CoxBoost
2. Boosting for high-dimensional time-to-event data with competing risks. Bioinformatics. 2009 Apr 1;25(7):890-6.
3. 恒峰基因：https://mp.weixin.qq.com/s/hkKemN59x-Ja0jPl3BUiLA
4. 生信补给站：https://mp.weixin.qq.com/s/jLxIQQ5CLUbPGp4BVbb3Dg
5. 叉叉滴同学的生信笔记：https://mp.weixin.qq.com/s/LPYKxkBoadpdoAuYXhYAxA

注：若对内容有疑惑或者有发现明确错误的朋友，请联系后台(欢迎交流)。更多内容可关注公众号：生信方舟