使用R语言的LongituRF软件包实现纵向随机森林。

12
我有一些高维重复测量数据,并且我有兴趣使用随机森林模型来研究这种模型的适用性和预测效用。具体而言,我正在尝试实现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(样本标识符的向量)和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 [矩阵]:包含随机效应预测变量的q个的Nxq矩阵。" 在使用实际数据时,如何指定这个矩阵呢?

我理解传统上Z只是组变量的单热(二进制)编码(例如如此描述),因此DataLongGenerator()中的Z应该是nxG(471x6)的稀疏矩阵,对吗?

希望能够明确如何使用实际数据来指定Z参数。

编辑

我的具体示例如下,我有一个响应变量(Y)。样本(用id标识)随机分配到干预(I,干预或不干预)。高维度的特征集合 (X)。特征和响应在两个时间点(Time,基线和终点)进行了测量。我想用XI来预测Y,并且还想提取哪些特征对预测Y最重要(就像Capitaine等人在他们的论文中对HIV所做的那样)。

我将按以下方式调用REEMforest()

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

我应该使用什么来代替Z

1个回答

2
当函数 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²。因此,如果您有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(多效果随机森林)函数没有任何问题。即使在研究论文中,该方法也很可靠。只是觉得值得一提。


在你提到癌症治疗的例子中,Z是什么? - JustGettinStarted
在我的特定示例中,lmm模型看起来像是Y ~ treatment + biomarker + treatment:time + biomarker:time + (1|subject)。我可以轻松地指定Y、X、时间和主题,但我不知道如何指定Z。 - JustGettinStarted
此方法的重点是在高维空间中工作。该论文模拟了HIV情境下的基因表达,因此系统地建模随机效应及其相互作用是不可行的吗? - JustGettinStarted
Z是随机效应,因此Z可以是“白细胞计数*体重”而不是“白细胞计数+体重”。 - Kat
1
那么字面上来说,这是一个包含逐元素乘积/求和列1(截距)、“白细胞计数*体重”和“白细胞计数+体重”的矩阵?在我的特定示例中,我只想为每个样本(即受试者ID)分配随机效应。 - JustGettinStarted
Kat - 你的编辑与其他人关于REEMforest()出现错误的经验一致,但不适用于其他函数。请参见此处答案1的评论 https://stats.stackexchange.com/questions/103730/how-can-i-include-random-effects-or-repeated-measures-into-a-randomforest/495055?noredirect=1#comment960438_495055 - JustGettinStarted

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接