在 R 中使用 LongituRF 包实现纵向随机森林

Question

我有一些高维重复测量数据，我有兴趣拟合随机森林模型来研究这些模型的适用性和预测效用。具体来说，我正在尝试实现LongituRF 包中的方法。这个包背后的方法在这里详细介绍：

Capitaine, L., et al. Random forests for high-dimensional longitudinal data. Stat Methods Med Res (2020) doi:10.1177/0962280220946080.

作者方便地提供了一些有用的数据生成函数进行测试。所以我们有

install.packages("LongituRF")
library(LongituRF)

让我们使用DataLongGenerator() 生成一些数据，它以 n=样本大小、p=预测变量的数量和 G=具有时间行为的预测变量的数量作为参数。

my_data <- DataLongGenerator(n=50,p=6,G=6)

my_data 是您期望 Y（响应向量）的列表， X（固定效应预测矩阵），Z（随机效应预测矩阵）， id（样本标识符向量）和 time（时间测量向量）。简单地拟合随机森林模型

model <- REEMforest(X=my_data$X,Y=my_data$Y,Z=my_data$Z,time=my_data$time,
                    id=my_data$id,sto="BM",mtry=2)

这里大约需要 50 秒，请耐心等待

到目前为止一切顺利。现在我很清楚这里的所有参数，除了Z。 什么是 Z 什么时候我要根据我的实际数据来拟合这个模型？

看着my_data$Z。

dim(my_data$Z)
[1] 471   2
head(my_data$Z)
      [,1]      [,2]
 [1,]    1 1.1128914
 [2,]    1 1.0349287
 [3,]    1 0.7308948
 [4,]    1 1.0976203
 [5,]    1 1.3739856
 [6,]    1 0.6840415

每一行看起来像一个截距项（即 1）和取自均匀分布 runif() 的值。

REEMforest() 的文档表明“Z [矩阵]：一个 Nxq 矩阵，包含随机效应的 q 预测器。” 在使用实际数据时如何指定这个矩阵？

我的理解是，传统上 Z 只是组变量的单热（二进制）编码（例如 as described here），所以 DataLongGenerator() 中的 Z 应该是 nxG (471x6) 稀疏矩阵不？

如果能明确说明如何使用实际数据指定 Z 参数，我们将不胜感激。

编辑

我的具体例子如下，我有一个响应变量（Y）。样本（用id标识）被随机分配到干预组（I，干预或不干预）。一组高维特征 (X)。在两个时间点（Time，基线和终点）测量特征和响应。我有兴趣预测Y，使用X 和I。我也有兴趣提取哪些特征对预测 Y 最重要（就像 Capitaine 等人在他们的论文中对 HIV 所做的那样）。

我会拨打REEMforest()如下

REEMforest(X=cbind(X,I), Y=Y, time=Time, id=id)

Z 应该用什么？

Answer 1

当函数DataLongGenerator() 创建Z 时，它是一个矩阵中的随机均匀数据。实际编码是

Z <- as.matrix(cbind(rep(1, length(f)), 2 * runif(length(f))))

其中f 表示代表每个元素的矩阵的长度。在您的示例中，您使用了 6 组 50 名参与者和 6 个固定效果。这导致长度为 472。

据我所知，由于此函数旨在模拟纵向数据，因此这是对该数据的随机效应的模拟。如果您使用的是真实数据，我认为它会更容易理解。

虽然这个示例没有使用 RE-EM 森林，但我认为它非常清楚，因为它使用有形元素作为示例。您可以在第 1.2.2 节固定 v. 随机效应中阅读有关随机效应的信息。 https://ademos.people.uic.edu/Chapter17.html#32_fixed_effects

查看第 3.2 节，了解使用真实数据时可以有意建模的随机效应示例。

另一个例子：您正在进行一项抗癌药物试验。您每周收集患者人口统计数据：体重、体温和 CBC 面板以及不同的给药组：每天 1 个单位、每天 2 个单位和每天 3 个单位。

在传统回归中，您需要对这些变量进行建模以确定模型识别结果的准确程度。固定效应是可解释方差或R²。因此，如果您有 0.86 或 86%，那么 14% 是无法解释的。它可能是一种相互作用导致了噪音，完美与模型确定的结果之间无法解释的差异。

假设白细胞计数非常低和超重的患者对治疗的反应要好得多。或者红头发的患者反应更好；这不在您的数据中。就纵向数据而言，假设关系（交互关系）仅在经过一段时间后才出现。

您可以尝试对不同的关系进行建模，以评估数据中的随机交互。不过，我认为您最好采用多种方法来系统地评估交互，而不是随机尝试识别随机效应。

已编辑我开始用 @JustGettinStarted 在 cmets 中写这个，但它太多了。

没有背景 - 实现此目的的最简单方法是运行 REEMtree::REEMtree() 之类的东西，将随机效果参数设置为 random = ~1 | time / id)。运行后，提取它计算的随机效应。你可以这样做：

data2 <- data %>% mutate(oOrder = row_number()) %>%  # identify original order of the data
     arrange(time, id) %>%
     mutate(zOrder = row_number()) # because the random effects will be in order by time then id

extRE <- data.frame(time = attributes(fit$RandomEffects[2][["id"]])[["row.names"]]) %>% 
     separate(col = time,
              into = c("time", "id"),
              sep = "\\/") %>%
     mutate(Z = fit$RandomEffects[[2]] %>% unlist(),
            id = as.integer(id),
            time = time))        # set data type to match dataset for time

data2 <- data2 %>% left_join(extRE) %>% arrange(oOrder) # return to original order

Z = cbind(rep(1, times = nrows(data2)), data2$Z)

另外，我建议您从随机生成随机效果开始。你开始的随机效果只是一个起点。最后的随机效果会不一样。

无论我尝试将LongituRF::REEMforest() 与真实数据使用多少种方式，我都遇到了错误。我每次都有一个不可逆的矩阵失败。

我注意到DataLongGenerator() 生成的数据按id 排序，然后是时间。我试图以这种方式订购数据（和 Z），但它没有帮助。当我从包LongituRF 中提取所有功能时，我毫无问题地使用了MERF（多效随机森林）功能。即使在研究论文中，这种方法也是可靠的。只是觉得值得一提。