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