时间戳

首页 > 解决方案 > lmer中方差矩阵的提取

r - lmer中方差矩阵的提取

问题描述

我试图从 lmer 中提取随机效应的方差分量，但在一个非常具体和复杂的设置中，即 fMRI。

我将首先解释问题，然后提供一个可重现的示例。

我正在处理功能性 fMRI 数据和测试聚合技术。在这个特定的上下文中，我有 3 个主题，每个主题都有一个包含第一级估计的矩阵。这里每一行代表大脑中的一个特定节点i，系数是节点ij的回归参数（暂时它们如何计算或为什么无关紧要。我想做的是使用线性聚合结果混合建模。

问题是，至少从我的经验来看，在 lmer 中使用矩阵几乎是不可能的。所以这就是我所做的：首先，我创建了一些随机矩阵，然后将它们转换为向量

   set.seed(49830)
   S.1 <- matrix(rexp(400, rate = 10),20,20)
   S.2 <- matrix(rexp(400, rate = 10),20,20)
   S.3 <- matrix(rexp(400, rate = 10),20,20)

   S1.vec <- as.vector(S.1)
   S2.vec <- as.vector(S.2)
   S3.vec <- as.vector(S.3)
   S.Vec.all <- c(S1.vec, S2.vec, S3.vec)

其次，我正在创建一个固定效应（假设每个人都平等的平均激活）和分组变量。在这种情况下，我只是创建了一个以数据为中心的向量。center_scale函数的代码可以在以下位置找到： http ://www.gastonsanchez.com/visually-enforced/how-to/2014/01/15/Center-data-in-R/

     constant <- as.vector(center_scale(S.Vec.all)) 
     Group1 <- rep(1, length(S1.vec))
     Group2 <- rep(2, length(S2.vec))
     Group3 <- rep(3, length(S3.vec))
     Group.all <- c(Group1, Group2, Group3)

第三，制作数据集，指定分组变量为因子，运行模型。

            data.comp <- data.frame(Group.all, constant, S.Vec.all) 
            data.comp$Group.all <- as.factor(data.comp$Group.all)

            mod <- lmer(S.Vec.all ~ constant + (constant|Group.all),
            data = data.comp,
            REML = T,
            control = lmerControl(calc.derivs = FALSE))

            summary(mod)
            ranef(mod)$Group.all

现在，问题是，即使我找到了一种方法来提取随机效应并将它们放在矩阵格式中，我仍然需要一种方法来计算这些单独矩阵的方差以进行进一步计算。

有没有办法提取随机系数的方差并将它们再次放入矩阵格式？我找到了一些伪代码来计算它们，使用 Laird and Ware(1982) 的公式，这是我想使用的公式，但它仍然给了我一个令人难以置信的大矩阵，这对我来说是完全无法理解的。这是代码：

            var.d <- crossprod(getME(mod,'Lambdat'))
            Zt <- getME(mod,"Zt")
            vr <- sigma(mod)^2
            var.b <- vr*(t(Zt) %*% var.d %*% Zt)
            sI <- vr * Diagonal(nrow(data.comp))
            var.y <- var.b + sI

我希望解释足够清楚，非常感谢所有可能并愿意帮助我的人！

标签： rmatrixhierarchyhierarchical-datalme4

解决方案

推荐阅读