The library languageR provides a method (pvals.fnc) to do MCMC significance testing of the fixed effects in a mixed effect regression model fit using lmer. However, pvals.fnc gives an error when the lmer model includes random slopes.

Is there a way to do an MCMC hypothesis test of such models?
If so, how? (To be accepted an answer should have a worked example in R) If not, is there a conceptual/computation reason why there is no way?

This question might be related to this one but I didn't understand the content there well enough to be certain.

Edit 1: A proof of concept showing that pvals.fnc() still does 'something' with lme4 models, but that it doesn't do anything with random slope models.

#the example from pvals.fnc
# 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...
#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)

It says: Error in pvals.fnc(primingHeid.lmer.rs) : MCMC sampling is not yet implemented in lme4_0.999375 for models with random correlation parameters

Additional question: Is pvals.fnc performing as expected for random intercept model? Should the outputs be trusted?

  • 3
    $\begingroup$ (1) I'm surprised that pvals.fnc still works at all; I thought mcmcsamp (which pvals.fnc relies on) had been non-functional in lme4 for quite a while. What version of lme4 are you using? (2) There is no conceptual reason why having random slopes should break whatever one is doing to get a significance test (3) Combining significance testing with MCMC is a little bit incoherent, statistically, although I understand the urge to do so (getting confidence intervals is more supportable) (4) the only relationship between this Q & the other is 'MCMC' (i.e. none, practically) $\endgroup$
    – Ben Bolker
    Dec 15, 2010 at 16:42
  • $\begingroup$ The version of lme4 I use depends on the computer I am sitting at. This console has lme4_0.999375-32, but I seldom use this one for analysis. I noticed several months ago that pvals.fnc() was ripping apart the lme4 results after analysis - I built a work around for it at the time, but it doesn't seem to be an issue anymore. I'll have to ask another question on your 3rd point in the near future. $\endgroup$ Dec 16, 2010 at 4:12

2 Answers 2


It looks like your error message isn't about varying slopes, it is about correlated random effects. You can fit the uncorrelated as well; that is, a mixed-effects model with independent random effects:

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

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


Here's (at least most of) a solution with MCMCglmm.

First fit the equivalent intercept-variance-only model with 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:

  ## 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,
                                   data = primingHeid)        

This does give some sort of "MCMC p-values" ... you'll have to explore for yourself and see whether the whole thing makes sense ...

  • $\begingroup$ Thanks a lot Ben. This looks like it will point me in the right direction. I just need to spend some time reading through the help and related articles for MCMCglmm to see if I can wrap my head around what is happening. $\endgroup$ Dec 20, 2010 at 8:26
  • 1
    $\begingroup$ Does the random slopes model in this case have a correlation between the random slope and random intercept, or in this framework is such an idea nonsensical? $\endgroup$ Feb 12, 2013 at 15:34
  • $\begingroup$ I tweaked your code here to make it easier to read, @Ben; if you don't like it, feel free to roll it back w/ my apologies. $\endgroup$ Aug 24, 2013 at 15:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.