# I want to use pvals.fnc() to get p-values for a lmer() model but cannot get rid of correlations between random factors

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:

>m
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) :
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?

• This may be helpful. Jul 25, 2013 at 6:18

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

EDIT

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.

library(languageR)
library(lme4)
data(primingHeid)


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.

\EDIT

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.

• +1. I especially like the "p-value fetishism" scepticism. Jul 26, 2013 at 1:10
• @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". – Sol Jul 26, 2013 at 15:34 • 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)
– Sol
Jul 26, 2013 at 15:34
• 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!
– Sol
Jul 26, 2013 at 15:36
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.