0

我有一个带有 、 和 变量的数据框subject，wd以及group一个value响应变量。每个受试者被分配到一组，每个工作日进行 7 次测量。因为每个主题都完全嵌套在一个组中，所以我想使用嵌套随机效应模型subject和group，以及为 . 添加第三个随机效应wd。目前，我正在使用它来这样做：

model = lmer(value ~ 1+ (1|wd) + (1|group) + (1|subject), 
             data = dframe, REML = 0)

我在本指南的第 40 页找到了我基于此的代码。我都用过REML = TRUE和REML = 0。但是，当我使用 时VarCorr(model)$variances，我得到

Groups   Name        Std.Dev.  
subject  (Intercept) 94.9534363
wd       (Intercept) 42.5931401
group    (Intercept)  0.0015608
Residual              0.9589836

此组方差与我用于生成数据的代码冲突，该代码的组均值分别为 36.9、28.78 和 -15.269。当我查看预测随机效应（使用ranef）与真实随机效应的“残差”时，我得到的残差与它们所在的组具有非常高的相关性（如果我建模residuals ~ group，则 R 平方值将超过 0.9）。

如何在 R 中正确拟合嵌套随机效应模型？我更喜欢使用 lme4，但任何软件包都足够了。

这是我用来生成数据的代码：

library(dplyr)
generate_data <- function(n = 10, g = 3, seed = 1, mean.overall = 300,
                          sigma.g = 50, sigma.wd = 50, 
                          sigma.subject = 100, sigma. = 30) {
    set.seed(seed)
    means.wd = rnorm(7) * sigma.wd
    means.g = rnorm(g) * sigma.g
    means.subject = rnorm(n*g) * sigma.subject
    dframe = data.frame(subject = rep(1:(g*n), each = 7),
                        wd = rep(1:7, g*n), 
                        group = rep(1:g, each = (7*n)))
    dframe = mutate(dframe,
       value = mean.overall + means.wd[wd] +    
           means.subject[subject] + means.g[group] + rnorm(7*g*n),
       subject = factor(subject, levels = 1:(n*g)),
       wd = factor(wd), 
       group = factor(group, levels = 1:g))
    dframe$value = round(pmax(5,dframe$value))
    truefx = list(wd = means.wd, group = means.g, 
                  subject = means.subject)
    list(data = dframe, effects = truefx)
}

dframe = generate_data()$data
4

1 回答 1

0

我猜您希望分组作为固定效果，因为只有几个级别。此外，由于每个主题/工作日组合只有一个观察值，因此不需要在主题内嵌套工作日。如果是这样，您只需要

lmer(value ~ group + (1|subject), data = dframe)

尚不清楚工作日是否真正嵌套在主题中。如果所有科目的工作日都是相同的，那么其他的可能更合适。不过，对于 stats.stackexchange.com 来说，这是一个更好的问题。

如果你真的想嵌套这些，这样的事情可能会奏效。

lmer(value ~ (1|group/subject/wd), data = dframe)
于 2015-05-15T17:11:54.030 回答