随机效应因子的最小推荐组数是多少?

机器算法验证 混合模式 样本量
2022-02-14 18:10:36

我在R( lme4) 中使用混合模型来分析一些重复测量数据。我有一个响应变量(粪便的纤维含量)和 3 个固定效应(体重等)。我的研究只有 6 名参与者,每人有 16 次重复测量(尽管两个只有 12 次重复)。受试者是在不同“治疗”中给予不同食物组合的蜥蜴。

我的问题是:我可以使用主题 ID 作为随机效果吗?

我知道这是纵向混合效应模型中通常的做法,考虑到受试者的随机抽样性质以及受试者内部的观察比受试者之间的观察更密切相关的事实。但是,将受试者 ID 视为随机效应涉及估计该变量的均值和方差。

  • 由于我只有 6 个受试者(该因子的 6 个水平),这足以准确表征均值和方差吗?

  • 在这方面,我对每个主题都有很多重复测量这一事实是否有帮助(我不明白这有什么关系)?

  • 最后,如果我不能将受试者 ID 作为随机效果,将其作为固定效果包含在我是否可以控制我重复测量的事实?

编辑:我想澄清一下,当我说“我可以”使用主题 ID 作为随机效果时,我的意思是“这样做是个好主意”。我知道我可以用一个只有 2 个级别的因子来拟合模型,但这肯定是站不住脚的吗?我在问什么时候考虑将受试者视为随机效应变得明智?似乎文献建议 5-6 级是一个下限。在我看来,在有 15 个以上的因子水平之前,随机效应的均值和方差的估计不会很精确。

4个回答

简短回答:是的,您可以使用 ID 作为 6 个级别的随机效果。

稍长的答案:@BenBolker 的 GLMM 常见问题解答在标题“我应该将因子 xxx 视为固定的还是随机的? ”下说(除其他外)以下内容

与“现代”混合模型估计(而不是“经典”矩量法估计)特别相关的一点是,出于实际目的,必须有合理数量的随机效应水平(例如块)——多于最少5个或6个。

所以你在下限,但在它的右边。

为了弄清楚多级模型的最小组数,我查看了 Gelman 和 Hill(2007)的《使用回归和多级/分层模型的数据分析》一书。

他们似乎在第 11 章第 5 节(第 247 页)中讨论了这个主题,他们写道,当组数小于 5 时,多级模型通常比经典模型添加的少。然而,他们似乎写道,应用多级模型几乎没有风险。

同样的作者似乎在第 12 章第 9 节(第 275-276 页)中回到了这个主题。他们在那里写道,关于多级模型的最少组数的建议是错误的。他们再次说,当组数较少时,多级模型通常比经典模型添加的少。然而,他们还写道,多级模型的表现不应该比无池回归差(其中无池似乎意味着在经典回归中使用组指标)。

在第 275-276 页上,作者有一个特定的小节针对一两个群体的情况(例如,男性与女性)。在这里,他们写道,他们通常以经典形式表达模型。然而,他们指出,即使只有一两个组,多级建模也很有用。他们写道,通过一组或两组多级建模可以简化为经典回归。

我对此的印象是,经典回归是模型连续体的一个末端,即多级模型的一个特例。

基于上述,我的印象是,当只有两组时,经典回归和多级建模将返回几乎相同的估计值,并且使用只有一组、二、三、四、五或六组的多级模型是可以的。

将来我将尝试使用R代码和一个小型数据集来修改这个答案,该数据集比较使用两组时使用两种方法获得的估计值。

对于它的价值,我做了一些模拟研究来查看相对简单的 LMM 的方差估计的稳定性(使用通过提供的sleepstudy数据集lme4)。第一种方法为主题ngroups数量生成所有可能的主题组合,并为每个可能的组合重新拟合模型。第二个需要几个随机的主题子集。

library(lme4)
library(ggplot2)
library(tidyr)

m0 <- lmer(Reaction ~ Days + (1|Subject), data = sleepstudy,
           control = lmerControl(optimizer = "nloptwrap"))
# set the number of factor levels
ngroups <- 3:18 
# generate all possible combinations
combos <- lapply(X = ngroups, 
                 FUN = function(x) combn(unique(sleepstudy$Subject), x)) 

# allocate output (sorry, this code is entirely un-optimized)
out <- list(matrix(NA, ncol(combos[[1]]), 1), matrix(NA, ncol(combos[[2]]), 1),
            matrix(NA, ncol(combos[[3]]), 1), matrix(NA, ncol(combos[[4]]), 1),
            matrix(NA, ncol(combos[[5]]), 1), matrix(NA, ncol(combos[[6]]), 1),
            matrix(NA, ncol(combos[[7]]), 1), matrix(NA, ncol(combos[[8]]), 1),
            matrix(NA, ncol(combos[[9]]), 1), matrix(NA, ncol(combos[[10]]), 1),
            matrix(NA, ncol(combos[[11]]), 1), matrix(NA, ncol(combos[[12]]), 1),
            matrix(NA, ncol(combos[[13]]), 1), matrix(NA, ncol(combos[[14]]), 1),
            matrix(NA, ncol(combos[[15]]), 1), matrix(NA, ncol(combos[[16]]), 1))
# took ~ 2.5 hrs on my laptop, commented out for safety
#system.time(for(ii in 1:length(combos)) {
#    for(jj in 1:ncol(combos[[ii]])) {
#    sls <- sleepstudy[sleepstudy$Subject %in% combos[[ii]][,jj],]
#    out[[ii]][jj] <- attr(VarCorr(update(m0, data = sls))$Subject, 'stddev')
#        }
#    })

# pad with zeros, not all were equal
# from http://stackoverflow.com/questions/11148429/r-convert-asymmetric-list-to-matrix-number-of-elements-in-each-sub-list-diffe
max.len <- max(sapply(out, length))
corrected.list <- lapply(out, function(x) {c(x, rep(NA, max.len - length(x)))})
mat <- do.call(rbind, corrected.list)
mat <- data.frame(t(mat))
names(mat) <- paste0('s',3:18)
mat <- gather(mat, run, value)

ggplot(mat, aes(x = value, fill = run)) + 
    geom_histogram(bins = 60) +
    geom_vline(xintercept = 37.12, linetype =  'longdash', 
               aes(colour = 'original')) +
    facet_wrap(~run, scales = 'free_y') +
    scale_x_continuous(breaks = seq(0, 100, by = 20)) + 
    theme_bw() + 
    guides(fill = FALSE)

黑色虚线是方差的原始点估计,分面代表不同数量的受试者(s3三人一组,s4四人等)。 在此处输入图像描述

另一种方法:

ngroups <- 3:18
reps <- 500
out2<- matrix(NA, length(ngroups), reps)

for (ii in 1:length(ngroups)) {
    for(j in 1:reps) {
        sls <- sleepstudy[sleepstudy$Subject %in% sample(unique(sleepstudy$Subject), ngroups[i], replace = FALSE),]
        out2[i,j] <- attr(VarCorr(update(m0, data = sls))$Subject, 'stddev')
    }
}
out2 <- data.frame(t(out2))
names(out2) <- paste0('s',3:18)
out2 <- gather(out2, run, value)

ggplot(out2, aes(x = value, fill = run)) + 
    geom_histogram(bins = 60) +
    geom_vline(xintercept = 37.12, linetype =  'longdash', 
               aes(colour = 'original')) +
    facet_wrap(~run, scales = 'free_y') +
    scale_x_continuous(breaks = seq(0, 100, by = 20)) + 
    theme_bw() + 
    guides(fill = FALSE)

在此处输入图像描述

看来(对于这个例子,无论如何)方差并没有真正稳定,直到至少有 14 个受试者,如果不是更晚的话。

您还可以使用贝叶斯混合模型 - 在这种情况下,在计算 95% 预测可信区间时会充分考虑随机效应估计中的不确定性。例如,新的 R 包brms和函数brm允许从lme4频率混合模型到贝叶斯模型的非常容易的转换,因为它具有几乎相同的语法。