r - 使用 anova() 测试 lmer() 模型中的随机效应时，是否需要设置 refit=FALSE？

Question

我目前正在测试是否应该在我的 lmer 模型中包含某些随机效应。我为此使用方差分析函数。到目前为止，我的过程是通过使用 REML=TRUE(默认选项)对 lmer() 的函数调用来拟合模型。然后，我对两个模型调用 anova()，其中一个模型包含要测试的随机效应，而另一个则不包含。然而，众所周知，anova() 函数通过 ML 重新拟合模型，但在新版本的 anova() 中，您可以阻止 anova() 通过设置选项 refit=FALSE 来避免这样做。为了测试随机效应，我是否应该在对 anova() 的调用中设置 refit=FALSE ？(如果我确实设置了 refit=FALSE code> p 值往往较低。当我设置 refit=FALSE 时，p 值是否是反保守的？)

方法一:

    mod0_reml <- lmer(x ~ y + z + (1 | w), data=dat)
    mod1_reml <- lmer(x ~ y + z + (y | w), data=dat)
    anova(mod0_reml, mod1_reml)

这将导致 anova() 使用 ML 而不是 REML 重新拟合模型。 (较新版本的 anova() 函数也将输出有关此的信息。)

方法2:

    mod0_reml <- lmer(x ~ y + z + (1 | w), data=dat)
    mod1_reml <- lmer(x ~ y + z + (y | w), data=dat)
    anova(mod0_reml, mod1_reml, refit=FALSE)

这将导致 anova() 在原始模型上执行计算，即使用 REML=TRUE。

为了测试我是否应该包含随机效应，这两种方法中哪一种是正确的？

感谢您的帮助

Answer 1

一般来说，我会说在这种情况下使用 refit=FALSE 是合适的，但让我们继续尝试模拟实验。

首先将没有随机斜率的模型拟合到 sleepstudy 数据集，然后模拟来自该模型的数据:

library(lme4)
mod0 <- lmer(Reaction ~ Days + (1|Subject), data=sleepstudy)
## also fit the full model for later use
mod1 <- lmer(Reaction ~ Days + (Days|Subject), data=sleepstudy)
set.seed(101)
simdat <- simulate(mod0,1000)

现在使用完整模型和简化模型重新拟合空数据，并保存使用和不使用 refit=FALSE 时由 anova() 生成的 p 值的分布。这本质上是原假设的参数引导测试；我们想看看它是否具有适当的特征(即 p 值的均匀分布)。

sumfun <- function(x) {
    m0 <- refit(mod0,x)
    m1 <- refit(mod1,x)
    a_refit <- suppressMessages(anova(m0,m1)["m1","Pr(>Chisq)"])
    a_no_refit <- anova(m0,m1,refit=FALSE)["m1","Pr(>Chisq)"]
    c(refit=a_refit,no_refit=a_no_refit)
}

我喜欢 plyr::laply 因为它很方便，尽管您也可以轻松地使用 for 循环或其他 *apply 之一> 方法。

library(plyr)
pdist <- laply(simdat,sumfun,.progress="text")

library(ggplot2); theme_set(theme_bw())
library(reshape2)
ggplot(melt(pdist),aes(x=value,fill=Var2))+
     geom_histogram(aes(y=..density..),
        alpha=0.5,position="identity",binwidth=0.02)+
     geom_hline(yintercept=1,lty=2)
ggsave("nullhist.png",height=4,width=5)

alpha=0.05 的 I 类错误率:

colMeans(pdist<0.05)
##   refit no_refit 
##   0.021    0.026

您可以看到，在这种情况下，两个过程给出了几乎相同的答案，并且两个过程都非常保守，众所周知的原因与假设的零值有关检验是在其可行空间的边界上进行的。对于测试单个简单随​​机效应的具体情况，将 p 值减半给出适当的答案(参见 Pinheiro 和 Bates 2000 等)；这实际上似乎在这里给出了合理的答案，尽管它并不真正合理，因为这里我们删除了两个随机效应参数(斜率的随机效应以及斜率和截距随机效应之间的相关性):

colMeans(pdist/2<0.05)
##   refit no_refit 
##   0.051    0.055 

其他要点:

• 您也许可以使用 pbkrtest 包中的 PBmodcomp 函数进行类似的练习。
• RLRsim 软件包专门针对随机效应项的零假设进行快速随机化(参数引导)测试而设计，但在这种稍微复杂的情况下似乎不起作用
• 查看相关GLMM faq section 类似的信息，包括为什么您可能根本不想测试随机效应的重要性的论点......
• 为了获得额外的积分，您可以使用偏差(-2 对数似然)差异而不是 p 值作为输出来重做参数 Bootstrap 运行，并检查结果是否符合 chi^2_0之间的混合code>(点质量为 0)和 chi^2_n 分布(其中 n可能 2，但我不确定这个几何形状)