package lme4
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-0002.pdf>