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
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
