Wie kann ein MCMC-Hypothesentest an einem Mixed-Effect-Regressionsmodell mit zufälligen Steigungen durchgeführt werden?

12

Die Bibliothek languageR bietet eine Methode (pvals.fnc) zum Testen der MCMC-Signifikanz der festen Effekte in einem Regressionsmodell für gemischte Effekte unter Verwendung von lmer. Pvals.fnc gibt jedoch einen Fehler aus, wenn das lmer-Modell zufällige Steigungen enthält.

Gibt es eine Möglichkeit, einen MCMC-Hypothesentest für solche Modelle durchzuführen?
Wenn das so ist, wie? (Um akzeptiert zu werden, sollte eine Antwort ein funktionierendes Beispiel in R haben.) Wenn nicht, gibt es einen begrifflichen / rechnerischen Grund, warum es keinen Weg gibt?

Diese Frage könnte mit dieser zusammenhängen, aber ich habe den Inhalt dort nicht gut genug verstanden, um sicherzugehen.

Edit 1 : Ein Proof-of-Concept zeigt, dass pvals.fnc () mit lme4-Modellen immer noch "etwas" macht, mit Modellen mit zufälliger Steigung jedoch nichts.

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)

Es heißt: Fehler in pvals.fnc (primingHeid.lmer.rs): Die MCMC-Abtastung ist in lme4_0.999375 für Modelle mit zufälligen Korrelationsparametern noch nicht implementiert

Zusätzliche Frage: Funktioniert pvals.fnc wie erwartet für ein zufälliges Intercept-Modell? Sollten die Ausgaben vertrauenswürdig sein?

russellpierce
quelle
3
(1) Ich bin überrascht, dass pvals.fnc überhaupt noch funktioniert. Ich dachte, dass mcmcsamp (auf das sich pvals.fnc stützt) in lme4 schon eine ganze Weile nicht mehr funktioniert. Welche Version von lme4 benutzt du? (2) Es gibt keinen begrifflichen Grund, warum zufällige Steigungen brechen sollten, was immer man tut, um einen Signifikanztest zu erhalten intervalle sind besser zu unterstützen) (4) Die einzige Beziehung zwischen diesem und dem anderen ist 'MCMC' (dh praktisch keine)
Ben Bolker
Die Version von lme4, die ich verwende, hängt von dem Computer ab, an dem ich sitze. Diese Konsole hat lme4_0.999375-32, aber ich benutze diese selten zur Analyse. Ich habe vor einigen Monaten bemerkt, dass pvals.fnc () die lme4-Ergebnisse nach der Analyse auseinandergerissen hat - ich habe damals eine Lösung dafür entwickelt, aber es scheint kein Problem mehr zu sein. Ich werde in naher Zukunft eine weitere Frage zu Ihrem dritten Punkt stellen müssen.
Russellpierce

Antworten:

4

Es sieht so aus, als ob es bei Ihrer Fehlermeldung nicht um unterschiedliche Steigungen geht, sondern um korrelierte zufällige Effekte. Sie können das unkorrelierte auch passen; Das heißt, ein Modell mit gemischten Effekten mit unabhängigen zufälligen Effekten:

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

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

Patrick McCann
quelle
8

Hier ist (zumindest die meisten) eine Lösung mit MCMCglmm.

Passen Sie zuerst das äquivalente Intercept-Varianz-Only-Modell an MCMCglmm:

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

Vergleichen von Übereinstimmungen zwischen MCMCglmmund lmer, erstmaliges Abrufen meiner gehackten Version von 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)

Versuchen Sie es jetzt mit zufälligen Steigungen:

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

Dies gibt eine Art "MCMC p-Werte" ... Sie müssen selbst herausfinden, ob das Ganze Sinn macht ...

Ben Bolker
quelle
Vielen Dank, Ben. Das sieht so aus, als würde es mich in die richtige Richtung lenken. Ich muss nur ein bisschen Zeit damit verbringen, die Hilfe und die zugehörigen Artikel für MCMCglmm zu lesen, um zu sehen, ob ich mich mit dem Geschehen beschäftigen kann.
Russellpierce
1
Hat das Zufallssteigungsmodell in diesem Fall eine Korrelation zwischen der Zufallssteigung und dem zufälligen Achsenabschnitt oder ist eine solche Idee in diesem Rahmen unsinnig?
Russellpierce
Ich habe deinen Code hier optimiert, um das Lesen zu vereinfachen, @Ben; Wenn es dir nicht gefällt, kannst du es gerne mit meiner Entschuldigung zurücksetzen.
gung - Wiedereinsetzung von Monica