My problem can be summarized very simply: I'm using a linear mixed-effects model and I am trying to get p-values using pvals.fnc(). The problem is that this function seems to have some trouble estimating p-values directly from the t-values associated with model coefficients (Baayen et al., 2008), and I don't know what is going wrong with the way I do it (i.e. according to what I have read, it should work). So, I'm explaining my model below, and if you can point out what I am doing wrong and suggest changes I would really appreciate it!

DESCRIPTION: I have a 2 by 2 within subjects design, fully crossing two categorical factors, "Gram" and "Number", each with two levels. This is the command I used to run the model:

>m <- lmer(RT ~ Gram*Number + (1|Subject) + (0+Gram+Number|Subject) + (1|Item),data= data)

If I understand this code, I am getting coefficients for the two fixed effects (Gram and Number) and their interaction, and I am fitting a model that has by-subject intercepts and slopes for the two fixed effects, and a by-item intercept for them. Following Barr et al. (2013), I thought that this code gets rid of the correlation parameters. I don't want estimate the correlations because I want to get the p-values using pvals.fnc (), and I read that this function doesn't work if there are correlations in the model.

The command seems to work:

Linear mixed model fit by REML 
Formula: RT ~ Gram * Number + (1 | Subject) + (0 + Gram + Number | Subject) + (1 |Item) 
   Data: mverb[mverb$Region == "06v1", ] 
   AIC   BIC logLik deviance REMLdev
 20134 20204 -10054    20138   20108
Random effects:
 Groups      Name        Variance  Std.Dev. Corr          
 Item       (Intercept)   273.508  16.5381               
 Subject     Gramgram        0.000   0.0000               
             Gramungram   3717.213  60.9689    NaN        
             Number1        59.361   7.7046    NaN -1.000 
 Subject     (Intercept) 14075.240 118.6391               
 Residual                35758.311 189.0987               
Number of obs: 1502, groups: Item, 48; Subject, 32

Fixed effects:
             Estimate Std. Error  t value
(Intercept)    402.520     22.321  18.033
Gram1          -57.788     14.545  -3.973
Number1         -4.191      9.858  -0.425
Gram1:Number1   15.693     19.527   0.804

Correlation of Fixed Effects:
            (Intr) Gram1  Numbr1
Gram1       -0.181              
Number1     -0.034  0.104       
Gram1:Nmbr1  0.000 -0.002 -0.011

However, when I try to calculate the p-values I still get an error message:

>pvals.fnc(m, withMCMC=T)$fixed
Error in pvals.fnc(m, withMCMC = T) : 
MCMC sampling is not implemented in recent versions of lme4
  for models with random correlation parameters

Am I making a mistake when I specify my model? Shouldn't pvals.fnc() work if I removed the correlations?

  • $\begingroup$ This may be helpful. $\endgroup$
    – martin
    Jul 25, 2013 at 6:18

2 Answers 2


(Italics represent corrected text)

You are making a 'mistake' in your model specification given what you say you want.

Random effects:
 Groups      Name        Variance  Std.Dev. Corr          
 Item       (Intercept)   273.508  16.5381               
 Subject     Gramgram        0.000   0.0000               
             Gramungram   3717.213  60.9689    NaN        
             Number1        59.361   7.7046    NaN -1.000 

You see the numbers there under Corr for Subject? That shows that you are estimating correlations between the random slope of Gramgungram and Gramgram, Number1 and Gramgram, and Number1 and Gramungram by subject. If Gram was a numeric , you could eliminate the random correlation between Gram and Number1 with a model specified by:

m <- lmer(RT ~ Gram*Number + (1|Subject) + (0+Gram|Subject) (0+Number|Subject) + (1|Item),data= data)

You'll notice that any random effect specified in the same set of parentheses yields a random effect correlation. At least that is true for models without the / symbol, I'm not really familiar with the notion for lmer when that is in the mix.

However, given what we saw from a model where you estimated this parameter, I'd advise caution. Moreover, you'll probably note that my code above doesn't work for you.


For those of you just now joining our program... for these examples I'll refer to primingHeid as OP did in the comments, this dataset can be found in languageR.


Why doesn't my code work? It doesn't work because Gram is a factor. Think about it for a little bit... and look at your fixed effects. If a factor has two levels how many parameters must you estimate to explain its effects? Two. Of course, one of the parameters you estimate is the intercept. The interpretation of the intercept will depend on how your factors are coded. In treatment coding (the default in R), the intercept represents the value for a case when all variables are at level 1 (cf. a regression textbook for details of other contrasts). Regardless of your contrasts two parameters are estimated for two levels of a factor. I think what is happening is that when you fail to specify an intercept R is protecting you from yourself and going ahead and estimating two parameters anyway. Try summary(lm(RT ~ 0 + Condition,data=primingHeid)) and you'll see that it went right on ahead and estimated two parameters. So, back to the context of lmer... if you have a factor with two levels, R will gladly estimate two parameters and then correlate them all under the hood. Back again to your comments... estimate lmer(RT ~ Condition +(0+Condition|Subject),data=primingHeid) and look and the ranef of that model and you'll see yet again that this is exactly what R has done.

If you wanted to force R to stop that, you'd have to do the factor coding manually by turning Condition into a numeric. The assumptions you'd have to make about the mean value of RT when Condition was at the level you coded as 0 are likely untenable (i.e. that RT is really 0). I won't exclude the possibility that with some careful thought, transformation of the DV (centering on the mean of condition you are setting equal to 0?), and good model specification you might work your way somewhere that made some sense... but, that would be an entirely different question and I can't speak to it at the moment.


I think you probably should step back and think about your model structure a little more (which really is one of the great messages in Barr et al., 2013). Are items crossed with gram and number? How many items occur within a unique arrangement of gram and number per subject?

More general issues now...

I have a huge amount of respect for Barr (no surprises there). However, he is not entirely mainstream on issues related to fitting this type of model. That isn't a bad thing ... but time will tell whether his approach for these models will become the next big thing. I have little doubt that 'keeping it maximal' is great if your data will tolerate it. But sometimes it won't. The backwards selection procedure he published involving the use of non-converged models is a bit unexpected. However, I have to admit now that I've seen his appendices I'm a little less sour on the idea than I was when I first read it. All the same, I'd like to see it be vetted a bit more.

You'll note Barr specifically does not use pvals.fnc() for models that have random correlations. So, only having skimmed the published version of his paper, I'd guess you can only use it under his approach if you can backwards step to a point where you don't have any.

Going now to my training with other stats gurus I feel compelled to say that almost all of this worry is an exercise in p-value fetishism that may be entirely misplaced - especially if you consider that this level of nested decision making yields a test that has a definition that is difficult to define.

  • $\begingroup$ +1. I especially like the "p-value fetishism" scepticism. $\endgroup$
    – usεr11852
    Jul 26, 2013 at 1:10
  • $\begingroup$ @rpierce: Thanks for the suggestions! I tried refitting your model, and the same message about the random effects correlation occured with pvals.fnc(). I also did length(levels(data$Gram)) for both Gram and Number and the answer is always 2. I think it might be a more general problem than my own data coding. To double-check this and to give you something tha you can verify yourself, I tried using a built in dataframe in R called "primingHeid". $\endgroup$
    – Sol
    Jul 26, 2013 at 15:34
  • $\begingroup$ This dataframe has RTs to a word depending on a preceding prime, and I chose two factors that I could check had two levels only: Condition ("baseheid/heid") and ResponseToPrime ("incorrect/correct"). You can check that the contrasts are coded correctly doing contrasts(primingHeid$Condition). This is the code I tried (it is identical to the one you gave me, just applied to this dataset): m <- lmer(RT ~ Condition*ResponseToPrime +(1|Subject) + (0+Condition|Subject) + (0+ResponseToPrime|Subject) + (1|Word) + (0+Condition|Word) + (0+ResponseToPrime|Word), data = primingHeid) $\endgroup$
    – Sol
    Jul 26, 2013 at 15:34
  • $\begingroup$ You can check yourself: there are still numbers under Corr in the model output, and the same error occurs when I do pvals.fnc(m, withMCMC=T)$fixed). Do you think this error message might be due to a more general issue with lmer()? I'm not sure what else to try. Thanks! $\endgroup$
    – Sol
    Jul 26, 2013 at 15:36
  • 1
    $\begingroup$ @SolLago: Edited the answer to address your question and deleted my comments that lead you down the wrong path. $\endgroup$ Jul 27, 2013 at 0:53

Luke (2016) Evaluating significance in linear mixed-effects models in R reports that the most optimal (certainly most conservative) test is p values based on Kenward-Roger approximation for degrees of freedom (in lmer). With large samples, p values based on the likelihood ratio (through anova()) are just as good.


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.