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