lme4 and languageR compatibility error: "input model is not a mer object”

Question

I have a dataset dat2 on which I want to fit a linear mixed-effects model. I used lmer() (package lme4) in the past in complement of pvals.fnc to compute the associated p-values.

However, since I reinstalled R 3.0.2 version with the new lme4 (1.0.4) and languageR (1.4) packages I obtain an error about the output of lmer function. It says that the output is not a mer object. Indeed its class is lmeRmod.

Here is the code I use:

names(dat2)<-c("auc","subj","decod","soa","vis")
attach(dat2)
mod1 <- lmer(auc ~ decod + (1 | subj))
mod2 <- lmer(auc ~ vis+ (1 | subj))
mod3 <- lmer(auc ~ decod + vis + (1 | subj))
mod4 <- lmer(auc ~ decod + vis + decod*vis + (1 | subj))
pvals.fnc(mod1)

And I get this error:

> pvals.fnc(mod1)
the input model is not a mer object
NULL

Indeed, when I look at mod1, I find it is a lmeRmod object and not a mer object.

> mod1
Linear mixed model fit by REML ['lmerMod']
Formula: auc ~ decod + (1 | subj) 
REML criterion at convergence: -213.3884 
Random effects:
 Groups   Name        Std.Dev.
 subj     (Intercept) 0.04187 
 Residual             0.11087 
Number of obs: 155, groups: subj, 6
Fixed Effects:
(Intercept)       decod2       decod3       decod4  
     0.9798      -0.1141      -0.3599      -0.3090 

This problem is very similar to the one discribed here. Any ideas of 1/ what the problem might be (why do I do not output a mer object) and 2/ how to work around it (I tried reinstalling older versions but I have compatibility problems between packages) ?

Any help would be great ! thanks !

Answer 1

I confirm that the pvals.fnc function doesn't work in the new languageR -- this is essentially because mcmcsamp wasn't implemented in the new version of lme4, which in turn because it was found to be unreliable in many cases.

We (the lme4 authors) are sorry to have left languageR users in the lurch this way, but it was somewhat unavoidable.

https://github.com/lme4/lme4/blob/master/man/pvalues.Rd offers some alternative suggestions for what to do about computing p-values.

https://github.com/lme4/lme4/blob/master/man/drop1.merMod.Rd gives a particular recipe (for the development version of lme4) about how to use pbkrtest::KRmodcomp to get p-values for all of the predictors in a model:

 fm1 <- lmer(Reaction~Days+(Days|Subject),sleepstudy)
 ## Likelihood ratio test
 drop1(fm1,test="Chisq")
 if (require(pbkrtest)) {
    KRSumFun <- function(object, objectDrop, ...) {
       krnames <- c("ndf","ddf","Fstat","p.value","F.scaling")
       r <- if (missing(objectDrop)) {
           setNames(rep(NA,5),krnames)
       } else {
          krtest <- KRmodcomp(object,objectDrop)
          unlist(krtest$stats[krnames])
       }
       attr(r,"method") <- c("Kenward-Roger via pbkrtest package")
       r
    }
    drop1(fm1,test="user",sumFun=KRSumFun)
}

This example produces:

Single term deletions

Model:
Reaction ~ Days + (Days | Subject)
Method: 
Kenward-Roger via pbkrtest package


       ndf ddf  Fstat    p.value F.scaling
<none>                                    
Days     1  17 45.853 3.2638e-06         1

Answer 2

You can use Package ‘lmerTest’ to get the p_values. See the example below:

#import lme4 package and lmerTest package
library(lmerTest)
# an object of class merModLmerTest
m <- lmer(Informed.liking ~ Gender+Information+Product +(1|Consumer), data=ham)
# gives summary of lmer object. The same as of class merMod but with
# additional p-values calculated based on Satterthwate's approximations
summary(m)

More about Package ‘lmerTest’, please see the link below: http://cran.r-project.org/web/packages/lmerTest/lmerTest.pdf

Answer 3

You might look into the mixed function in the afex package found here. It uses Kenward-Rogers for df.