[R] Problem with lmer and wiki example

[Philippe Grosjean]

James Widman wrote:
> I am trying to duplicate the example by Spencer Graves in the wiki, 
> using lmer with the Nozzle data.
> http://wiki.r-project.org/rwiki/doku.php?id=guides:lmer-tests
> However the Chisq value and the fitAB values that are calculated are 
> different compared to those in the example. I also get a warning message 
> when I attempt the fitAB. Does anyone have any guidance as to why this 
> might happen and how to correct it?
> I am using R on Kubutu in case that may be helpful.
> Thanks
> 
> ---- my code --
> [Previously saved workspace restored]
>  > rm(list=ls(all=TRUE))
>  > list=ls(all=TRUE)
>  > print(list)
> character(0)
>  > y <- c(6,6,-15, 26,12,5, 11,4,4, 21,14,7, 25,18,25,
> + 13,6,13, 4,4,11, 17,10,17, -5,2,-5, 15,8,1,
> + 10,10,-11, -35,0,-14, 11,-10,-17, 12,-2,-16, -4,10,24)
>  > Nozzle <- data.frame(Nozzle=rep(LETTERS[1:3], 
> e=15),Operator=rep(letters[1:5], e=3), flowRate=y)
>  > summary(Nozzle)
> Nozzle Operator flowRate
> A:15 a:9 Min. :-35.000
> B:15 b:9 1st Qu.: 0.000
> C:15 c:9 Median : 7.000
> d:9 Mean : 5.511
> e:9 3rd Qu.: 13.000
> Max. : 26.000
>  > library(lme4)
> Loading required package: Matrix
> Loading required package: lattice
>  > fitAB <- lmer(flowRate~Nozzle+(Nozzle|Operator),data=Nozzle, 
> method="ML")
> Warning messages:
> 1: In .local(x, ..., value) :
> Estimated variance-covariance for factor ‘Operator’ is singular
> 
> 2: In .local(x, ..., value) :
> nlminb returned message false convergence (8)
> 
>  > fitB <- lmer(flowRate~1+(1|Operator), data=Nozzle, method="ML")
>  > anova(fitAB, fitB)
> Data: Nozzle
> Models:
> fitB: flowRate ~ 1 + (1 | Operator)
> fitAB: flowRate ~ Nozzle + (Nozzle | Operator)
> Df AIC BIC logLik Chisq Chi Df Pr(>Chisq)
> fitB 2 359.36 362.98 -177.68
> fitAB 9 362.13 378.39 -172.06 11.237 7 0.1286
> ----------------------------------------------------------------------
> Output from Spencer Graves example
> fitAB 9 359.88 376.14 -170.94 13.479 7 0.06126

Now, on a Mac OS X (using the unstable, development version of R 2.9.0, 
and recompiled version of lme4_0.999375-28... so caution of course!), I 
got this:

First, the method = "ML" argument is deprecated and replaced by REML = 
TRUE/FALSE, but the doc at ?lmer does not tell exactly what is the 
equivalence to method = "ML" (and I don't know enough in this field to 
determine it by myself). Anyway, I tried both:

 > fitAB <- lmer(flowRate~Nozzle+(Nozzle|Operator),data=Nozzle, REML = TRUE)
#Warning messages:
#1: In .local(x, ..., value) :
#Estimated variance-covariance for factor Operator is singular

#2: In .local(x, ..., value) :
#nlminb returned message false convergence (8)

 > fitB <- lmer(flowRate~1+(1|Operator), data=Nozzle, REML = TRUE)
 > anova(fitAB, fitB)
Data: Nozzle
Models:
fitB: flowRate ~ 1 + (1 | Operator)
fitAB: flowRate ~ Nozzle + (Nozzle | Operator)
       Df     AIC     BIC  logLik  Chisq Chi Df Pr(>Chisq)
fitB   3  361.36  366.78 -177.68
fitAB 10  362.10  380.17 -171.05 13.261      7    0.06601 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 > fitAB <- lmer(flowRate~Nozzle+(Nozzle|Operator),data=Nozzle, REML = 
FALSE)
#Warning messages:
#1: In .local(x, ..., value) :
#Estimated variance-covariance for factor Operator is singular

#2: In .local(x, ..., value) :
#nlminb returned message false convergence (8)

 > fitB <- lmer(flowRate~1+(1|Operator), data=Nozzle, REML = FALSE)
 > anova(fitAB, fitB)
Data: Nozzle
Models:
fitB: flowRate ~ 1 + (1 | Operator)
fitAB: flowRate ~ Nozzle + (Nozzle | Operator)
       Df     AIC     BIC  logLik  Chisq Chi Df Pr(>Chisq)
fitB   3  361.36  366.78 -177.68
fitAB 10  361.96  380.03 -170.98 13.402      7     0.0629 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


So, I got same error messages as you. I got results closer to the one in 
the wiki page, BUT, I am puzzled by the degrees of freedom that are 
different 3/10 in my case, against 2/9 in yours and in the wiki page!

Could the authors of lmer(), and/or of the wiki page, or the code cited 
in the wiki page (in CC) provide some explanation to this? 
Corrections/updates of the wiki page so that it reflects latest lmer() 
version would be also very much appreciated.

All the best,

Philippe Grosjean

Just in case:

 > R.Version()
$platform
[1] "i386-apple-darwin8.11.1"

$arch
[1] "i386"

$os
[1] "darwin8.11.1"

$system
[1] "i386, darwin8.11.1"

$status
[1] "Under development (unstable)"

$major
[1] "2"

$minor
[1] "9.0"

$year
[1] "2009"

$month
[1] "01"

$day
[1] "22"

$`svn rev`
[1] "47686"

$language
[1] "R"

$version.string
[1] "R version 2.9.0 Under development (unstable) (2009-01-22 r47686)"

 > sessionInfo()
R version 2.9.0 Under development (unstable) (2009-01-22 r47686)
i386-apple-darwin8.11.1

locale:
en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] tcltk     stats     graphics  grDevices utils     datasets  methods
[8] base

other attached packages:
[1] lme4_0.999375-28   Matrix_0.999375-18 lattice_0.17-20
[4] svGUI_0.9-44       svSocket_0.9-43    svMisc_0.9-47

loaded via a namespace (and not attached):
[1] grid_2.9.0  tools_2.9.0