如何对具有随机斜率的混合效应回归模型进行 MCMC 假设检验？

机器算法验证 r 混合模式 统计学意义 蒙特卡洛
2022-02-27 05:58:15

库语言 R 提供了一种方法 (pvals.fnc) 使用 lmer 对混合效应回归模型拟合中的固定效应进行 MCMC 显着性测试。但是，当 lmer 模型包含随机斜率时，pvals.fnc 会出错。

有没有办法对这些模型进行 MCMC 假设检验？
如果是这样，怎么做？（要被接受的答案应该在 R 中有一个工作示例）如果没有，是否有概念/计算原因导致没有办法？

这个问题可能与这个问题有关，但我对那里的内容不太了解，无法确定。

编辑 1：概念证明表明 pvals.fnc() 仍然对 lme4 模型做“某些事情”，但它对随机斜率模型没有任何作用。

library(lme4)
library(languageR)
#the example from pvals.fnc
data(primingHeid) 
# remove extreme outliers
primingHeid = primingHeid[primingHeid$RT < 7.1,]
# fit mixed-effects model
primingHeid.lmer = lmer(RT ~ RTtoPrime * ResponseToPrime + Condition + (1|Subject) + (1|Word), data = primingHeid)
mcmc = pvals.fnc(primingHeid.lmer, nsim=10000, withMCMC=TRUE)
#Subjects are in both conditions...
table(primingHeid$Subject,primingHeid$Condition)
#So I can fit a model that has a random slope of condition by participant
primingHeid.lmer.rs = lmer(RT ~ RTtoPrime * ResponseToPrime + Condition + (1+Condition|Subject) + (1|Word), data = primingHeid)
#However pvals.fnc fails here...
mcmc.rs = pvals.fnc(primingHeid.lmer.rs)

它说： pvals.fnc(primingHeid.lmer.rs) 中的错误：对于具有随机相关参数的模型，MCMC 采样尚未在 lme4_0.999375 中实现

附加问题：对于随机截距模型，pvals.fnc 是否按预期执行？输出是否值得信任？

2个回答

这是（至少大部分）一个解决方案MCMCglmm。

首先拟合等效的仅截距方差模型MCMCglmm：

library(MCMCglmm)
primingHeid.MCMCglmm = MCMCglmm(fixed=RT ~ RTtoPrime * ResponseToPrime + Condition, 
                                random=~Subject+Word, data = primingHeid)

比较 和 之间的拟合MCMCglmm，lmer首先检索我的破解版本arm::coefplot：

source(url("http://www.math.mcmaster.ca/bolker/R/misc/coefplot_new.R"))
  ## combine estimates of fixed effects and variance components
pp  <- as.mcmc(with(primingHeid.MCMCglmm, cbind(Sol, VCV)))
  ## extract coefficient table
cc1 <- coeftab(primingHeid.MCMCglmm,ptype=c("fixef", "vcov"))
  ## strip fixed/vcov indicators to make names match with lmer output
rownames(cc1) <- gsub("(Sol|VCV).", "", rownames(cc1))
  ## fixed effects -- v. similar
coefplot(list(cc1[1:5,], primingHeid.lmer))
  ## variance components -- quite different.  Worth further exploration?
coefplot(list(cc1[6:8,], coeftab(primingHeid.lmer, ptype="vcov")),
         xlim=c(0,0.16), cex.pts=1.5)

现在尝试使用随机斜率：

primingHeid.rs.MCMCglmm = MCMCglmm(fixed=RT ~ RTtoPrime * ResponseToPrime + Condition,
                                   random=~Subject+Subject:Condition+Word, 
                                   data = primingHeid)        
summary(primingHeid.rs.MCMCglmm)

这确实给出了某种“MCMC p值”......你必须自己探索，看看整个事情是否有意义......

看起来您的错误消息不是关于变化的斜率，而是关于相关的随机效应。您也可以拟合不相关的；即具有独立随机效应的混合效应模型： 

Linear mixed model fit by REML
Formula: Reaction ~ Days + (1 | Subject) + (0 + Days | Subject)
Data: sleepstudy

来自http://www.stat.wisc.edu/~bates/IMPS2008/lme4D.pdf

其它你可能感兴趣的问题
上一篇测试重复测量方差分析的正态假设？（在 R 中） 下一篇估计动态线性模型的参数