3

I have two dependent variable Prime Type (five levels), and Prime Relatedness (two levels), and one dependent variable; Reaction Time (RT). I have ran a linear mixed-effect model in R using the following formula for lme4 package:

data.mod1=lmer(RT~Related*PrimeType*(1|Subject)+(1|Item),mydata)

Then I ran a similar formula for nlme package

data.mod1.lme=lme(RT~Related*PrimeType, random=~1|Subject/Item,mydata)

Considering that these models are analyzing the same data set using linear mixed effect models, the t-test show different values! In fact, for one of the variables it shows a significant p-value in nlme model but not in lme4 package:

nlme.variable t=-1.98707 p-value=0.0470
lme4.variable t=-1.14    p-value=0.2917866*

The p-value in lme4 is not calculated and sampMCMC is no longer available, I had to calculated using the following formula 2(1 - pt(abs(x), df)) --> 2*(1 - pt(abs(-1.14), 7))

My question is that why are the t-values differ between the two packages. Do I need to change or add something in R to have both models show the same t-values? Also, is my calculation of the p-value wrong? If yes, what is the correct way for it to be calculated it?

EDIT: nlme experts - what is the correct code for the above mentioned formula to add cross random effect?

ama
  • 63
  • 1
  • 1
  • 5
  • I think you should probably submit the "how do I specify crossed random effects" as a separate question on StackOverflow, if http://www.mail-archive.com/r-help@stat.math.ethz.ch/msg10849.html doesn't answer your question. – Ben Bolker Dec 23 '13 at 19:03

2 Answers2

3

You haven't given a reproducible example, but:

library(lme4)
library(nlme)
library(pbkrtest)

fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
## see definition of KRSumFun below, taken from ?drop1.merMod:
## requires recent [development?] version of `lme4`:

packageVersion("lme4")  ## 1.1.2
drop1(fm1,test="user",sumFun=KRSumFun)

## 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

fm2 <- lme(Reaction ~ Days , random = ~Days | Subject, sleepstudy)
anova(fm2)
##             numDF denDF   F-value p-value
##(Intercept)      1   161 1454.0766  <.0001
## Days            1   161   45.8534  <.0001

Note that this is a case (random-slopes model) where lme actually gets the wrong answer for the denominator degrees of freedom.

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
    }
Ben Bolker
  • 34,308
  • 2
  • 93
  • 126
  • library(lme4) # I updated it and it is on version'1.0.4' library(pbkrtest) mydata.mod1=lmer(RT~Related*PrimeType*(1|Item)+(1|Subject),mydata) drop1(mydata.mod1,test="user",sumFun=KRSumFun) # I recieved the following message "Error in match.arg(test) : 'arg' should be one of "none", "Chisq"" # @BenBolker what do you think the issue is? – ama Dec 23 '13 at 19:53
  • http://artax.karlin.mff.cuni.cz/r-help/library/afex/html/mixed.html – ama Dec 23 '13 at 19:59
  • # Also, I have tried the following the above post: mixed(mydata.mod1, mydata, type = 3, method = c("KR")) # and m1 – ama Dec 23 '13 at 19:59
  • `mixed` is from the `afex` package -- you need to install and load it. You need a more recent version of `lme4` (are you running R < 3.0?); try `install.packages("lme4",repos="http://lme4.r-forge.r-project.org/repos")` (if you are running R >3.0 -- if not this will be more difficult) – Ben Bolker Dec 23 '13 at 20:58
  • You were right, I had to update R and re-install lme4. Both the mixed function from afex package and drop1 function worked well. Thank you! – ama Dec 29 '13 at 13:26
1

I primarily use lme4, so you'll have to excuse me if this answer is incorrect. I believe you specified crossed-random effects in the lmer function call, but nested random effects in the lme function call. I'm not sure how to do crossed-random effects in lme, but perhaps this or this might yield some insight.

One other thing (unrelated), is that RTs usually possess a distribution with positive skew and could benefit from an inverse transformation. See my answer here for a little more info.

Community
  • 1
dmartin
  • 3,010
  • 3
  • 22
  • 27
  • I see, then I will look more into how the crossed-random effect is written in lme. @dmartin Since you use lme4, could you tell me how do you calculate the p-values? Also, what are the terms used in lme4 package to included nested random effects? I really appreciate it. – ama Dec 22 '13 at 21:19
  • 1
    1) Lme4 doesn't give you p-values, because for unbalanced, multi-level data, denominator degrees of freedom are unknown. Most packages give an approximation. Douglas Bates (lme4 author) instead chooses not to. More information found [here](https://stat.ethz.ch/pipermail/r-help/2006-May/094765.html) 2) Nested effects really depend on the structure of the data, so it's hard for me to see what you want to specify from the example. [Here](http://stats.stackexchange.com/questions/13166/rs-lmer-cheat-sheet) would be a good place to start – dmartin Dec 23 '13 at 16:54