[R] Problem with lmer and wiki example
James Widman
jwidman at mi.nmfs.gov
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
More information about the R-help
mailing list