在使用anova()的lmer()模型中测试随机效果时,是否需要设置refit = FALSE？

Question

我目前正在测试是否应该在我的lmer模型中包含某些随机效果.我使用anova函数.我的过程,到目前为止是一个函数调用模型适合lmer()用REML=TRUE(默认选项).然后我打电话anova()给两个模型,其中一个确实包括要测试的随机效果,另一个没有.但是,众所周知,该anova()函数使用ML重新生成模型,但在新版本中,anova()您可以anova()通过设置选项来防止这样做refit=FALSE.为了测试随机效应refit=FALSE,我应该在我的调用中设置anova() or not?(如果我设置refit=FALSEp值往往更低.当我设置时,p值是反保守的refit=FALSE吗？)

方法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.

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

谢谢你的帮助

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)

现在使用完整模型和简化模型重新调整null数据,并保存anova()with和without 生成的p值的分布refit=FALSE.这基本上是零假设的参数自举测试; 我们想看看它是否具有适当的特征(即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 

其他要点:

• 您可以使用包中的PBmodcomp函数执行类似的练习pbkrtest.
• 该RLRsim软件包的设计精确适用于关于随机效应项的零假设的快速随机化(参数引导)测试,但在这种稍微复杂的情况下似乎不起作用
• 请参阅相关的GLMM faq部分以获取类似信息,包括为什么您可能根本不想测试随机效应的重要性的论据...
• 为了额外的功劳,你可以使用偏差(-2对数似然)差异而不是p值作为输出重做参数自举运行,并检查结果是否符合a chi^2_0(点质量为0)和chi^2_n分布(其中)之间的混合n为大概 2,但我不会为这种几何形状可以肯定的)