On Sun, Nov 1, 2009 at 9:01 AM, wenjun zheng <wjzhen...@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@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. >
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
Rplots.pdf
Description: Adobe PDF document
______________________________________________ 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.
