[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
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> 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.
Type 'contributors()' for more information and
'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)
Loading required package: Matrix
Loading required package: lattice
> 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                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=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>


More information about the R-help mailing list