Jus*_*ted 9 regression r random-forest mixed-models longitudinal
我有一些高维重复测量数据，我对拟合随机森林模型感兴趣，以研究此类模型的适用性和预测效用。具体来说，我正在尝试实现LongituRF包中的方法。此处详细介绍了此包背后的方法：
Capitaine, L. 等。用于高维纵向数据的随机森林。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（样本标识符向量）和时间（时间测量向量）的列表。简单地拟合随机森林模型
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 [矩阵]：包含随机效应的 q 预测变量的 Nxq 矩阵。” 使用实际数据时如何指定这个矩阵？
我的理解是，传统Z是简单地一热（二进制）组变量（例如编码这里描述），所以Z从DataLongGenerator()应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什么？
当函数DataLongGenerator()创建时Z，它是矩阵中的随机均匀数据。实际的编码是
Z <- as.matrix(cbind(rep(1, length(f)), 2 * runif(length(f))))
其中f表示代表每个元素的矩阵的长度。在您的示例中，您使用了 6 组，每组 50 名参与者，具有 6 个固定效果。结果长度为 472。
据我所知，由于该函数旨在模拟纵向数据，因此这是对该数据的随机效应的模拟。如果您使用真实数据，我认为会更容易理解。
虽然这个例子没有使用 RE-EM 森林，但我认为它非常清楚，因为它使用有形元素作为示例。您可以在 1.2.2 固定效应与随机效应一节中阅读有关随机效应的内容。https://ademos.people.uic.edu/Chapter17.html#32_fixed_effects
请参阅第 3.2 节，了解随机效应的示例，如果您使用真实数据，您可以有意建模这些随机效应。
另一个例子：您正在进行一项癌症药物试验。您每周收集患者人口统计数据：体重、体温和 CBC 组以及不同的给药组：每天 1 个单位、每天 2 个单位和每天 3 个单位。
在传统回归中，您将对这些变量进行建模，以确定模型识别结果的准确程度。固定效应是可解释的方差或R 2。因此，如果您有 0.86 或 86%，那么 14% 是无法解释的。这可能是一种导致噪音的相互作用，即完美与模型确定的结果之间无法解释的差异。
假设白细胞计数非常低且超重的患者对治疗的反应要好得多。或者也许红头发的患者反应更好；这不在你的数据中。就纵向数据而言，假设关系（交互关系）仅在经过一定时间后才出现。
您可以尝试对不同的关系进行建模来评估数据中的随机交互作用。不过，我认为您最好使用多种方法来系统地评估相互作用，而不是随机尝试识别随机效应。
编辑我开始在@JustGettinStarted 的评论中写下这个，但是太多了。
如果没有背景 -实现此目的的最简单方法是运行 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（多重效应随机森林）函数，没有任何问题。即使在研究论文中，这种方法也是可靠的。只是觉得值得一提。