动态预测模型最新R包-JMbayes2包教程(9)-超级学习(最终章)

教程官网：https://drizopoulos./JMbayes2

本次教程地址：https://drizopoulos./JMbayes2/articles/Super_Learning.html

超级学习

目的和理论

包括了纵向数据和时间-事件数据的联合模型用于对相应的动态预测结果的计算。它可以随着新的可用信息不断加入而更新。因此在精准医疗取得广泛应用，尤其是癌症和心血管疾病中。以往的联合模型通常应用单一模型来获取动态预测，当然这有难度，特别是当计算多个纵向数据结果时。并且因为是动态预测，不同的模型可能在不同的随访时间点具备不同水平的准确性。所以，我们考虑采用多个联合模型进行动态预测，并将这些模型的结果结合起来，以优化准确性。我们使用超级学习（SL）的概念来实现这一点。SL是一种集成方法，允许研究人员将多种不同的预测算法组合成一个。它使用V折交叉验证构建候选算法库所产生的预测结果的最优加权组合，并且可以由使用者指定的目标函数定义需要优化的特征，例如最小化均方误差或最大化ROC曲线下面积。

V折交叉验证是一种用于评估模型性能和泛化能力的常用方法。在V折交叉验证中，将原始数据集分成V个相似大小的子集，然后进行V次迭代的训练和测试。在每一次迭代中，其中一个子集被保留作为测试集，而其余的V-1个子集被用作训练集来训练模型。这个过程重复进行V次，以确保每个子集都被用作测试集一次。通过多次迭代，可以获得多个性能评估指标的值，然后可以计算它们的平均值来得出对模型性能的综合评估。这有助于减少因随机性引起的性能评估偏差，并提供更可靠的性能估计。V的值通常是一个使用者来定义的参数，可以根据问题的性质和数据集的大小进行调整，常见的选择包括5折交叉验证和10折交叉验证。

超级学习背后的基本思想是推导出能够优化交叉验证预测的模型权重。更具体地说，我们让 表示一个包含 个模型的模型库。其中，这些联合模型在纵向子模型中的时间趋势规定和函数形式以及事件子模型等方面存在差异，但是这个库中包含的模型数量无需任何限制，甚至可以说多多益善。我们将原始数据集 分成 个折叠。 的值将取决于 的事件数和样本量，因为对于其中的每个折叠，我们需要有足够数量的事件来稳健地量化预测性能。我们在 个折叠的组合（训练集）中拟合了 个模型，并计算出来留在外面的第 个折叠（验证集）的预测结果，最后进行交叉验证。由于预测的动态性质，我们希望在不同的随访时间点得出最佳权重。具体来说，可以设置时间点序列 。这些时间点的数量和位置应结合 中可用的事件信息来考量。对于 ，我们定义 来表示在时间 处于风险状态且属于第 个折叠的受试者。对于 中的所有受试者，我们计算交叉验证预测结果，

这些预测结果基于模型库 中的模型 计算，其中，该模型是在训练集，即排除了第 个折叠中的患者的数据集 的基础上拟合的，计算方式基于蒙特卡洛方法。我们用 来表示 个模型的预测结果的凸组合，即：

其中，，对于 ，且 。请注意，权重 是随时间变化的，即在不同的时间点， 个模型的不同组合可以产生更加准确的预测。

1.蒙特卡洛方法（Monte Carlo method）是一种统计学和计算机科学领域常用的数值模拟方法。它的基本思想是通过随机抽样和概率统计来解决复杂的数学问题或进行模拟实验。蒙特卡洛方法的名称来源于摩纳哥蒙特卡洛市，因为这个方法最早是用来模拟赌博游戏的概率统计问题。主要步骤包括：1. 生成随机样本：根据问题的特性，生成符合某种概率分布的随机数样本。2. 进行模拟：使用生成的随机样本进行模拟实验或计算，以得出所需的结果。3. 统计分析：对多次模拟实验的结果进行统计分析，得出问题的估计值或概率分布。

蒙特卡洛方法广泛应用于众多领域，包括物理学、工程学、金融学、统计学、计算机科学等。它特别适用于处理复杂的多维积分、概率分布、随机过程等问题，可以用来解决无法用解析方法求解的问题。蒙特卡洛方法的精确性通常随着抽样次数的增加而提高，因此在计算资源允许的情况下，可以获得足够精确的结果。

2.凸组合（Convex Combination）凸组合在凸优化、凸几何、线性规划等领域有广泛的应用，它是描述和构建凸结构的基础概念。对于给定的向量或点集合{v₁, v₂, ..., vₖ}和对应的非负权重{α₁, α₂, ..., αₖ}，如果这些权重满足以下两个条件：1. 非负性，权重都不小于零。2. 归一性，即权重的总和等于1。那么，通过将这些向量或点按照权重线性组合起来，得到的结果就是一个凸组合。数学表示如下：

C = α₁v₁ α₂v₂ ... αₖvₖ

凸组合的一个重要特性是，结果位于连接原始点集合中所有点的线性组合的凸壳内。这意味着凸组合的结果不会超出原始点的凸包（最小凸多边形或凸多面体）。

对于任何时间 ，我们将使用评分规则来选择权重 ，以优化联合交叉验证预测模型的预测性能。如果真实分布达到了最优期望得分，那么这个评分规则 也被认为是合适的。例：

表示真实模型下的条件风险概率， 是 的估计值。期望是针对真实模型下的生存结果的条件密度 进行的，其中 。评分规则 的含义是更低的得分表示更高的准确性。

Brier分数是一种合适的评分规则，它结合了判别能力和校准来测量整体的预测性能。在随访时间 和医学相关的时间窗口 内，我们定义Brier分数为：

作为替代的评分规则，在区间 中，我们采用了期望预测交叉熵（EPCE）的一种改进方式：

其中期望基于真实模型下的 。在此处， 和 都是使用交叉验证预测 的凸组合来计算的。使用超级学习的过程中，我们得到了按评分规则（ ，或者 都可以）得分最小（对预测最好）的的权重 ，

权重的约束条件：，其中 ；以及 。

chatGPT解读：

该部分说明了统计学中如何使用超级学习（Super Learning）来优化预测模型的性能。主要内容包括：

1. 超级学习的基本思想：超级学习是一种集成方法，旨在将多个不同的预测算法结合成一个最优的组合。它使用交叉验证来建立来自候选算法库的预测的最佳权重组合，以优化性能度量。

2. 交叉验证和模型库：在超级学习中，数据集被分为多个折叠（folds），模型库包含多个不同的模型，这些模型可以在不同的折叠上进行训练和测试。

3. 时间序列预测和评分规则：超级学习用于动态预测时，使用适当的评分规则来衡量模型的性能。这些评分规则包括Brier Score（BS）和Expected Predictive Cross-Entropy（EPCE），它们用于比较模型的预测与真实观测之间的差异。

4. 权重优化：超级学习通过选择一组随着时间变化的权重，以将评分规则的值求最小来优化模型的预测性能。这些权重确保了不同模型在不同时间点的组合，以提高预测准确性。

5. 权重的约束：权重的选择受到一些约束的限制，包括权重必须为正数且总和为1，以确保权重的合法性和归一性。

总的来说，超级学习是一种强大的方法，用于优化动态预测模型的性能，特别是在需要考虑多个模型和不同时间点的情况下。通过选择最优的权重组合，可以提高预测模型的准确性。

举例

使用PBC数据集来展示将超级学习与联合模型进行结合并应用于动态预测。我们首先使用create_folds()函数将pbc2数据库分五折。

CVdats <- create_folds(pbc2, V = 5, id_var = 'id')


此函数的第一个参数是我们想要分成V个折叠的数据框。参数id_var指定了此数据集中代表受试者主体变量的名称id。create_folds()的输出是一个带有两个元素，分别名为training'和testing的列表（也就是训练集和验证集）。每个元素又是一个带有V个数据框的列表。

接下来，定义一个用于拟合动态预测联合模型的函数。

1. 这个函数需要有一个数据框(data.frame)参数，用于拟合联合模型。
2. 为优化计算性能，我们将使用并行计算来将这些模型拟合到不同的训练数据集上。因此，在函数内部，使用library('JMbayes2')来在每个工作进程中加载JMbayes2包。
3. 使用as.numeric（）进行二进制转换，表示复合时间（死亡或者一致）发生与否。！表示逻辑取反，即如果status不为alive，不是取1，而是取0。并且将二进制数据储存在status2中。
4. 创建一个data_id列，储存受试者标识id。用！duplicate（）去掉重复的id，是每个受试者只被考虑1次。
5. 使用lem（）函数，和纵向数据（serBilir）拟合线性混合模型。血清胆红素值经过对数转换，对数值作为响应变量，时间year作为解释变量。还包括随机效应，~ year为随机斜率，|id为随机截距，代表着每个受试者的基线情况不同，随时间变化的程度也不同。以及优化控制参数。
6. 
   

该函数的输出应该是一个拟合的联合模型的列表，其类别为'jmList'。将这个类别分配给结果列表将有助于以后合并预测。

fit_models <- function (data) {
    library('JMbayes2')
    data$status2 <- as.numeric(data$status != 'alive')
    data_id <- data[!duplicated(data$id), ]
    lmeFit <- lme(log(serBilir) ~ year, data = data,
                  random = ~ year | id,
                  control = lmeControl(opt = 'optim'))
    CoxFit <- coxph(Surv(years, status2) ~ 1, data = data_id)
    jmFit1 <- jm(CoxFit, lmeFit, time_var = 'year')
    jmFit2 <- update(jmFit1, 
                     functional_forms = ~ slope(log(serBilir)))
    jmFit3 <- update(jmFit1, 
                     functional_forms = ~ value(log(serBilir))   area(log(serBilir)))
    ###
    lmeFit2 <- lme(log(serBilir) ~ ns(year, 3, B = c(0, 14.4))   sex   age, 
                   data = data, random = ~ ns(year, 3, B = c(0, 14.4)) | id,
                   control = lmeControl(opt = 'optim'))
    CoxFit2 <- coxph(Surv(years, status2) ~ sex   age, data = data_id)
    jmFit4 <- jm(CoxFit2, lmeFit2, time_var = 'year')
    jmFit5 <- update(jmFit4, 
                     functional_forms = ~ slope(log(serBilir)))
    out <- list(M1 = jmFit1, M2 = jmFit2, M3 = jmFit3, M4 = jmFit4, M5 = jmFit5)
    class(out) <- 'jmList'
    out
}

具体而言，我们考虑了一个关于纵向结果serBilir（血清胆红素）和复合事件（移植或死亡）的单变量联合模型库。前三个模型考虑了血清胆红素的简单线性混合效应模型，其中每个受试者都有随机截距和随机斜率，并且无其他协变量。此外，在复合事件的Cox模型中，我们同样没有指定任何基线协变量；因此，复合事件的风险仅取决于血清胆红素。

最后的模型考虑了线性混合模型的更详细规范，其中包括自然立方样条(natural cubic splines)在固定效应和随机效应中，以允许对对数血清胆红素轨迹中的非线性进行建模，并在混合模型和Cox模型中都考虑了性别和年龄的主效应。这些功能形式仍然是当前值和当前斜率。

我们在训练数据集上使用并行计算来拟合这些模型，这可以通过使用parallel包来实现（请注意：这和后续的计算需要一些时间，具体取决于您的计算机性能；在我的计算机上，配备了Intel(R) Core(TM) i9-10885H CPU @ 2.40GHz和 32.0 GB RAM，大约需要20分钟来运行）：

cl <- parallel::makeCluster(5L)
Models_folds <- parallel::parLapply(cl, CVdats$training, fit_models)
parallel::stopCluster(cl)


我们计算权重来组合这五个模型的预测，以优化集成Brier分数和预期的交叉熵，针对随访时间t=6年和Δt=2年的相关窗口。函数tvBrier()通过提供作为第一个参数的在训练数据集中拟合的联合模型列表，自动执行此任务。集成Brier分数是使用在newdata参数中提供的测试数据集来计算的：

tstr <- 6
thor <- 8

Brier_weights <- tvBrier(Models_folds, newdata = CVdats$testing, 
                         integrated = TRUE, Tstart = tstr, Thoriz = thor)
Brier_weights
#> 
#> Cross-Validated Prediction Error using the Library of Joint Models 'Models_folds'
#> 
#> Super Learning Estimated Integrated Brier score: 0.056
#> In the time interval: [6, 8)
#> For the 166 subjects at risk at time 6
#> Number of subjects with an event in [6, 8): 18
#> Number of subjects with a censored time in [6, 8): 44
#> Accounting for censoring using model-based weights
#> 
#> Integrated Brier score per model: 0.0596 0.0588 0.0613 0.0629 0.065
#> Weights per model: 0.1134 0.565 0.3029 0.0086 0.0101
#> Number of folds: 5

我们观察到前三个模型占据了主导地位，权重较大。我们还注意到，基于组合预测的集成Brier分数小于每个单独模型的集成Brier分数。要使用预期的预测交叉熵来计算模型权重，可以使用tvEPCE()函数，其调用方式与Brier分数几乎相同：

EPCE_weights <- tvEPCE(Models_folds, newdata = CVdats$testing, 
                       Tstart = tstr, Thoriz = thor)
EPCE_weights
#> 
#> Cross-Validated Expected Predictive Cross-Entropy using the Library of Joint Models 'Models_folds'
#> 
#> Super Learning Estimated EPCE: 0.3109
#> In the time interval: [6, 8)
#> For the 166 subjects at risk at time 6
#> Number of subjects with an event in [6, 8): 18
#> Number of subjects with a censored time in [6, 8): 44
#> 
#> EPCE per model: 0.3568 0.3599 0.3639 0.3992 0.4321
#> Weights per model: 0.0028 0.5386 0.4586 0 0
#> Number of folds: 5


EPCE的结果与集成Brier分数的结果一致；然而，现在只有模型M2和M3共享权重。再次观察到，基于组合的交叉验证预测的EPCE小于基于每个单独模型的交叉验证预测的EPCE。

要将这些权重应用到实际问题中，首先需要在原始数据集上重新拟合我们考虑的五个联合模型。

Models <- fit_models(pbc2)

然后，我们构建一个数据集，其中包括在第六年处于风险状态的受试者，并考虑在此随访时间之前收集的纵向测量。此外，在这个数据集中，将观察到的事件时间设置为六，事件变量设置为零，即表示患者在此时间之前没有发生事件：

# 选择第六年处于风险状态的受试者
ND <- pbc2[pbc2$years > tstr & pbc2$year <= tstr, ]

# 删除不需要的列，并确保'id'列为标量而不是向量
ND$id <- ND$id[, drop = TRUE]

# 将'years'列设置为tstr
ND$years <- tstr

# 将'status2'列设置为0，表示在此时间之前没有发生事件
ND$status2 <- 0


作为示例，我们将使用EPCE权重来组合患者8的累积风险预测。我们可以使用'jmList'类对象的predict()方法来实现这一点，还需要提供权重参数。其余的参数与'jm'对象的predict()方法相同（也可以参考Dynamic Predictions文档）：

# 从EPCE权重中提取模型权重
model_weights <- EPCE_weights$weights

# 使用模型权重进行事件预测
predsEvent <- predict(Models, weights = model_weights, newdata = ND[ND$id == 8, ],
                      process = 'event', return_newdata = TRUE)

# 使用模型权重进行纵向预测
predsLong <- predict(Models, weights = model_weights, newdata = ND[ND$id == 8, ],
                     return_newdata = TRUE)

# 绘制纵向和事件预测的图形
plot(predsLong, predsEvent)