An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20091101/a2ec105c/attachment-0001.pl>
package lme4
6 messages · wenjun zheng, Douglas Bates
1 day later
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>
An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20091103/10966bed/attachment-0001.pl>
An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20091103/d17d5bc1/attachment-0001.pl>
On Tue, Nov 3, 2009 at 8:08 AM, wenjun zheng <wjzheng09 at gmail.com> wrote:
Thanks,Douglas, It really helps me a lot, but?is there any other way?if I want to show whether a random effect is significant in text file, like P value or other ?index. Thanks very much again. Wenjun.
Well there are p-values from the likelihood ratio tests in that transcript I sent. The point of those tests is that a p-value can only be calculated when you know both the null hypothesis and the alternative, which is why those p-values are the result of comparing two nested model fits.
2009/11/2 Douglas Bates <bates at stat.wisc.edu>
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.
On Tue, Nov 3, 2009 at 9:11 AM, wenjun zheng <wjzheng09 at gmail.com> wrote:
May be I can calculate p value by t testing approximately: ?1-qnorm(Variance/Std.Dev.)
That would be a z test, not a t test, wouldn't it? And it would only be meaningful if the distribution of the estimator is approximately normal, which we know it definitely is not. Estimates of variances have distributions like a chi-square, not like a normal. In particular, the estimators are not symmetrically distributed about the parameter value.
But which function can help me to extract Variance and Std.Dev values from the results below:
The VarCorr extractor produces the estimates of the variance. You would have to write your own functions to determine an approximate standard error of those estimates and I would not advise doing so. Firstly, the code would be rather complicated and secondly the answer would be of very limited usefulness, for the reasons discussed above.
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 2009/11/2 Douglas Bates <bates at stat.wisc.edu>
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.