[R] Problem with lmer and wiki example

Fri Feb 13 14:08:21 CET 2009

Philippe Grosjean wrote:
> 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
Phillipe,
My output from R ver 2.8.1 on a different PC running Windows XP agrees 
with your output verbatim, including the df difference.
I used the kubuntu repository when I installed R on the other machine, I 
guess I'll have to remove it and reinstall a version from CRAN.
Best Regards,
Jim