在使用anova()函数测试lmer()模型的随机效应时，是否需要将refit=FALSE设置为假？

Question

在使用anova()函数测试lmer()模型的随机效应时，是否需要将refit=FALSE设置为假？

13

我目前在测试是否应该在我的lmer模型中包含某些随机效应。为此，我使用anova函数。到目前为止，我的步骤是使用lmer()函数调用拟合模型，并使用REML=TRUE(默认选项)。然后，在其中一个模型上调用anova()，该模型包括要测试的随机效应，而另一个模型则不包括。然而，众所周知，anova()函数会重新拟合ML模型，但是在新版本的anova()中，您可以通过设置选项refit=FALSE来防止anova()这样做。为了测试随机效应，我应该在对anova()的调用中设置refit=FALSE还是不设置？（如果我设置refit=FALSE，则p值往往会更低。当我设置refit=FALSE时，p值是否反保守？）

方法1：

    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。

哪种方法是正确的，以测试是否应该包括随机效应？

感谢任何帮助。

- lord.garbage

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Ben Bolker · Accepted Answer

一般来说，在这种情况下使用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)

现在，使用完整模型和简化模型重新安装空数据，并使用anova()生成包括refit=FALSE和不包括refit=FALSE的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)

空分布的直方图

当α=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部分以获取类似信息，包括为什么您可能根本不想测试随机效应重要性的论据...
如果要获得额外的学分，您可以使用偏差（-2对数似然）差异而不是p值作为输出重新执行参数自助法运行，并检查结果是否符合介于chi^2_0（0处点质量）和chi^2_n分布（其中n可能是可能 2，但我不能确定针对此几何形状）的混合物。