Here's (at least most of) a solution with MCMCglmm
.
First fit the equivalent intercept-variance-only model with MCMCglmm
:
library(MCMCglmm)
primingHeid.MCMCglmm = MCMCglmm(fixed=RT ~ RTtoPrime * ResponseToPrime + Condition,
+ Condition,
random=~Subject+Word, data = primingHeid)
Comparing fits between MCMCglmm
and lmer
, first retrieving my hacked version of 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)
Now try it with random slopes:
primingHeid.rs.MCMCglmm = MCMCglmm(fixed=RT ~ RTtoPrime * ResponseToPrime
+ Condition,
random=~Subject+Subject:Condition+Word,
data = primingHeid)
summary(primingHeid.rs.MCMCglmm)
This does give some sort of "MCMC p-values" ... you'll have to explore for yourself and see whether the whole thing makes sense ...