Jak wykonać test hipotezy MCMC na modelu regresji z efektem mieszanym z losowymi nachyleniami?

Biblioteka languageR zapewnia metodę (pvals.fnc) do testowania istotności MCMC ustalonych efektów w modelu dopasowania regresji mieszanej przy użyciu lmera. Jednak pvals.fnc podaje błąd, gdy model Lmer zawiera losowe zbocza.

Czy istnieje sposób przeprowadzenia testu hipotetycznego MCMC dla takich modeli?
Jeśli tak to jak? (Aby zostać zaakceptowanym, odpowiedź powinna zawierać sprawdzony przykład w języku R). Jeśli nie, to czy istnieje powód koncepcyjny / obliczeniowy, dlaczego nie ma sposobu?

To pytanie może być związane z tym pytaniem, ale nie zrozumiałem tam wystarczająco dobrze, aby się upewnić.

Edycja 1 : Dowód koncepcji pokazujący, że pvals.fnc () nadal robi coś z modelami lme4, ale nie robi nic z przypadkowymi modelami nachylenia.

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)

Mówi: Błąd w pliku pvals.fnc (primingHeid.lmer.rs): Próbkowanie MCMC nie jest jeszcze zaimplementowane w lme4_0.999375 dla modeli z losowymi parametrami korelacji

Dodatkowe pytanie: Czy plik pvals.fnc działa zgodnie z oczekiwaniami dla modelu przechwytywania losowego? Czy należy ufać wyjściom?

Wygląda na to, że Twój komunikat o błędzie nie dotyczy różnych nachyleń, dotyczy skorelowanych efektów losowych. Możesz również dopasować do nieskorelowanych; to jest model mieszanych efektów z niezależnymi efektami losowymi:

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

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

Oto (przynajmniej większość) rozwiązanie z MCMCglmm.

Najpierw dopasuj równoważny model tylko wariantu przechwytującego z MCMCglmm:

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

Porównywanie pasowań między MCMCglmmi lmer, najpierw odzyskiwanie mojej zhakowanej wersji 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)

Teraz spróbuj z losowymi zboczami:

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

Daje to pewnego rodzaju „wartości p MCMC” ... musisz sam odkryć i zobaczyć, czy to wszystko ma sens ...