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

12

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?

russellpierce
źródło
3
(1) Jestem zaskoczony, że pvals.fnc nadal w ogóle działa; Myślałem, że mcmcsamp (na którym opiera się pvals.fnc) od dłuższego czasu nie działa w lme4. Jakiej wersji lme4 używasz? (2) Nie ma pojęciowego powodu, dla którego losowe zbocza powinny złamać wszystko, co się robi, aby uzyskać test istotności (3) Połączenie testu istotności z MCMC jest trochę niespójne statystycznie, chociaż rozumiem potrzebę zrobienia tego (uzyskanie pewności siebie) interwały są bardziej obsługiwane) (4) jedynym związkiem między tym
pytaniem
Wersja lme4, której używam, zależy od komputera, na którym siedzę. Ta konsola ma lme4_0.999375-32, ale rzadko używam tej do analizy. Zauważyłem kilka miesięcy temu, że pvals.fnc () rozrywa wyniki lme4 po analizie - w tym czasie stworzyłem obejście tego problemu, ale nie wydaje się to już problemem. W najbliższej przyszłości będę musiał zadać kolejne pytanie dotyczące twojego trzeciego punktu.
russellpierce

Odpowiedzi:

4

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

Patrick McCann
źródło
8

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

Ben Bolker
źródło
Wielkie dzięki Ben. Wygląda na to, że wskaże mi właściwy kierunek. Muszę tylko poświęcić trochę czasu na przeczytanie pomocy i powiązanych artykułów dla MCMCglmm, aby zobaczyć, czy mogę zawinąć głowę wokół tego, co się dzieje.
russellpierce
1
Czy model losowych nachyleń ma w tym przypadku korelację między losowym nachyleniem a przypadkowym przechwytywaniem, czy w tych ramach taki pomysł jest bezsensowny?
russellpierce
Poprawiłem twój kod tutaj, aby ułatwić czytanie, @Ben; jeśli ci się nie podoba, nie krępuj się przywrócić go z przeprosinami.
Gung - Przywróć Monikę