¿Cómo puede uno hacer una prueba de hipótesis MCMC en un modelo de regresión de efecto mixto con pendientes al azar?

Question

La biblioteca languageR proporciona un método (pvals.fnc) para hacer MCMC pruebas de significación de los efectos fijos en un modelo de regresión de efectos mixtos ajuste utilizando lmer. Sin embargo, pvals.fnc da un error cuando el lmer modelo incluye aleatoria de pistas.

Es allí una manera de hacer una MCMC prueba de hipótesis de estos modelos?
Si es así, ¿cómo? (Para ser aceptado una respuesta debe tener un ejemplo en R) Si no, hay conceptual/cálculo de la razón de por qué no hay ninguna forma?

Esta pregunta podría estar relacionado con este , pero yo no entendía el contenido no lo suficientemente bien como para ser ciertas.

Edit 1: Una prueba de concepto que muestra que pvals.fnc() sigue 'algo' con lme4 modelos, pero que no hace nada con el azar de la pendiente de los modelos.

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)

Dice: Error en pvals.fnc(primingHeid.lmer.el ld) : MCMC de muestreo no está aún implementado en lme4_0.999375 para los modelos con azar de correlación de los parámetros de

Pregunta adicional: Es pvals.fnc funcionando como se espera por azar interceptar modelo? Deben las salidas es de confianza?

Answer 1

Aquí (al menos la mayoría) una solución con MCMCglmm.

Primero ajuste el equivalente de intercepción de la varianza modelo sólo con MCMCglmm:

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

La comparación se ajusta entre MCMCglmm y lmer, la primera de recuperar mi versión hackeada de 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)

Ahora intente con aleatoria de pistas:

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

Esto da una especie de "MCMC p-valores" ... tendrás que explorar por sí mismo y ver si todo tiene sentido ...

Answer 2

Parece que el mensaje de error no se trata de diferentes pistas, se trata de efectos aleatorios correlacionados. Caben sin correlación es decir, un modelo de efectos mixtos con efectos aleatorios independientes:

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

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