[R] package lme4

Douglas Bates bates at stat.wisc.edu
Mon Nov 2 16:39:48 CET 2009

```On Sun, Nov 1, 2009 at 9:01 AM, wenjun zheng <wjzheng09 at gmail.com> wrote:
> Hi R Users,
>     When I use package lme4 for mixed model analysis, I can't distinguish
> the significant and insignificant variables from all random independent
> variables.
>     Here is my data and result:
> Data:
>
>  Rice<-data.frame(Yield=c(8,7,4,9,7,6,9,8,8,8,7,5,9,9,5,7,7,8,8,8,4,8,6,4,8,8,9),
>                 Variety=rep(rep(c("A1","A2","A3"),each=3),3),
>                 Stand=rep(c("B1","B2","B3"),9),
>                 Block=rep(1:3,each=9))
>    Rice.lmer<-lmer(Yield ~ (1|Variety) + (1|Stand) + (1|Block) +
> (1|Variety:Stand), data = Rice)
>
> Result:
>
> Linear mixed model fit by REML
> Formula: Yield ~ (1 | Variety) + (1 | Stand) + (1 | Block) + (1 |
> Variety:Stand)
>   Data: Rice
>   AIC   BIC logLik deviance REMLdev
>  96.25 104.0 -42.12    85.33   84.25
> Random effects:
>  Groups        Name        Variance Std.Dev.
>  Variety:Stand (Intercept) 1.345679 1.16003
>  Block         (Intercept) 0.000000 0.00000
>  Stand         (Intercept) 0.888889 0.94281
>  Variety       (Intercept) 0.024691 0.15714
>  Residual                  0.666667 0.81650
> Number of obs: 27, groups: Variety:Stand, 9; Block, 3; Stand, 3; Variety, 3

> Fixed effects:
>            Estimate Std. Error t value
> (Intercept)   7.1852     0.6919   10.38

> Can you give me some advice for recognizing the significant variables among
> random effects above without other  calculating.

Well, since the estimate of the variance due to Block is zero, that's
probably not one of the significant random effects.

Why do you want to do this without other calculations?  In olden days
when each model fit involved substantial calculations by hand one did
try to avoid fitting multiple models but now that is not a problem.
You can get a hint of which random effects will be significant by
looking at their precision in a "caterpillar plot" and then fit the
reduced model and use anova to compare models.  See the enclosed

>    Any suggestions will be appreciated.
> Wenjun
>
>        [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> and provide commented, minimal, self-contained, reproducible code.
>
-------------- next part --------------

R version 2.10.0 Patched (2009-11-01 r50288)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

Natural language support but running in an English locale

R is a collaborative project with many contributors.
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(lme4)
> Rice <- data.frame(Yield = c(8,7,4,9,7,6,9,8,8,8,7,5,9,9,5,7,
+                    7,8,8,8,4,8,6,4,8,8,9),
+                    Variety = gl(3, 3, 27, labels = c("A1","A2","A3")),
+                    Stand = gl(3, 1, 27, labels = c("B1","B2","B3")),
+                    Block = gl(3, 9))
> print(fm1 <- lmer(Yield ~ 1 + (1|Block) + (1|Stand) +
+                   (1|Variety) + (1|Variety:Stand), Rice))
Linear mixed model fit by REML
Formula: Yield ~ 1 + (1 | Block) + (1 | Stand) + (1 | Variety) + (1 |      Variety:Stand)
Data: Rice
AIC   BIC logLik deviance REMLdev
96.25 104.0 -42.12    85.33   84.25
Random effects:
Groups        Name        Variance Std.Dev.
Variety:Stand (Intercept) 1.34568  1.16003
Variety       (Intercept) 0.02469  0.15713
Stand         (Intercept) 0.88889  0.94281
Block         (Intercept) 0.00000  0.00000
Residual                  0.66667  0.81650
Number of obs: 27, groups: Variety:Stand, 9; Variety, 3; Stand, 3; Block, 3

Fixed effects:
Estimate Std. Error t value
(Intercept)   7.1852     0.6919   10.38
> ## Estimate of Block variance is zero, remove that term
> print(fm2 <- lmer(Yield ~ 1 + (1|Stand) + (1|Variety) + (1|Variety:Stand),
+                   Rice))
Linear mixed model fit by REML
Formula: Yield ~ 1 + (1 | Stand) + (1 | Variety) + (1 | Variety:Stand)
Data: Rice
AIC   BIC logLik deviance REMLdev
94.25 100.7 -42.12    85.33   84.25
Random effects:
Groups        Name        Variance Std.Dev.
Variety:Stand (Intercept) 1.345679 1.16003
Variety       (Intercept) 0.024692 0.15714
Stand         (Intercept) 0.888888 0.94281
Residual                  0.666667 0.81650
Number of obs: 27, groups: Variety:Stand, 9; Variety, 3; Stand, 3

Fixed effects:
Estimate Std. Error t value
(Intercept)   7.1852     0.6919   10.38
> ## Very small estimate of variance for Variety
> ## Check the precision of the random effects
> pl2 <- dotplot(ranef(fm2, postVar = TRUE))
> pl2[[1]]                                # perhaps significant
> pl2\$Variety                             # not significant
> pl2\$Stand                               # probably not significant
>
> print(fm3 <- lmer(Yield ~ 1 + (1|Stand) + (1|Variety:Stand),
+                   Rice))
Linear mixed model fit by REML
Formula: Yield ~ 1 + (1 | Stand) + (1 | Variety:Stand)
Data: Rice
AIC   BIC logLik deviance REMLdev
92.25 97.43 -42.12    85.31   84.25
Random effects:
Groups        Name        Variance Std.Dev.
Variety:Stand (Intercept) 1.37037  1.17063
Stand         (Intercept) 0.88066  0.93844
Residual                  0.66667  0.81650
Number of obs: 27, groups: Variety:Stand, 9; Stand, 3

Fixed effects:
Estimate Std. Error t value
(Intercept)   7.1852     0.6859   10.48
> anova(fm3, fm2)
Data: Rice
Models:
fm3: Yield ~ 1 + (1 | Stand) + (1 | Variety:Stand)
fm2: Yield ~ 1 + (1 | Stand) + (1 | Variety) + (1 | Variety:Stand)
Df    AIC     BIC  logLik Chisq Chi Df Pr(>Chisq)
fm3  4 93.315  98.498 -42.657
fm2  5 95.331 101.810 -42.665     0      1          1
> pl3 <- dotplot(ranef(fm3, postVar = TRUE))
> pl3[[1]]                                # perhaps significant
> pl3\$Stand                               # probably not significant
>
> print(fm4 <- lmer(Yield ~ 1 + (1|Variety:Stand), Rice))
Linear mixed model fit by REML
Formula: Yield ~ 1 + (1 | Variety:Stand)
Data: Rice
AIC   BIC logLik deviance REMLdev
91.07 94.96 -42.53     85.5   85.07
Random effects:
Groups        Name        Variance Std.Dev.
Variety:Stand (Intercept) 2.03086  1.4251
Residual                  0.66667  0.8165
Number of obs: 27, groups: Variety:Stand, 9

Fixed effects:
Estimate Std. Error t value
(Intercept)   7.1852     0.5003   14.36
> anova(fm4, fm3)
Data: Rice
Models:
fm4: Yield ~ 1 + (1 | Variety:Stand)
fm3: Yield ~ 1 + (1 | Stand) + (1 | Variety:Stand)
Df    AIC    BIC  logLik  Chisq Chi Df Pr(>Chisq)
fm4  3 91.504 95.391 -42.752
fm3  4 93.315 98.498 -42.657 0.1888      1     0.6639
> sessionInfo()
R version 2.10.0 Patched (2009-11-01 r50288)
x86_64-unknown-linux-gnu

locale:
[1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=C              LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8       LC_NAME=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

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

other attached packages:
[1] lme4_0.999375-32   Matrix_0.999375-32 lattice_0.17-26

loaded via a namespace (and not attached):
[1] grid_2.10.0
>
> proc.time()
user  system elapsed
13.700   0.170  13.846
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Rplots.pdf
Type: application/pdf
Size: 14503 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20091102/8dfa9a72/attachment.pdf>
```