wenqing li
2005-Apr-16 18:38 UTC
[R] help on extract variance components from the fitted model by lm
Hey, all: Do we have a convenient command(s) to extract the variance components from a fitted model by lm (actually it's a nexted model)? e.g.: using the following codes we could get MSA,MSB(A) and MSE. How to get the variance component estimates by command in R rather than calculations by hand? A<-as.vector(rep(c(rep(1,5), rep(2,5), rep(3,5), rep(4,5), rep(5,5)),2)) B<-as.vector(rep(c(rep(c(1,2,3,4,5),5)),2)) y<-as.vector(c(15.5,15.2,14.2,14.3,15.8,6.2,7.2,6.6,6.2,5.6,15.4,13.9,13.4,12.5,13.2,10.9,12.5,12.3,11.0,12.3,7.5,6.7,7.2,7.6,6.3,14.9,15.2,14.2,14.3,16.4,7,8.4,7.8,7.6,7.4,14.4,13.3,14.8,14.1,15,11.3,12.7,11.7,12,13.3,6.7,7.3,6,7.6,7.1)) lm1<-lm(y~factor(A)/factor(B)) anova(lm1) Thanks a lot! And have a good weekend! Regards, Wenqing
Marc Schwartz
2005-Apr-16 20:31 UTC
[R] help on extract variance components from the fitted model by lm
On Sun, 2005-04-17 at 02:38 +0800, wenqing li wrote:> Hey, all: Do we have a convenient command(s) to extract the variance > components from a fitted model by lm (actually it's a nexted model)? > > e.g.: using the following codes we could get MSA,MSB(A) and MSE. How > to get the variance component estimates by command in R rather than > calculations by hand? > > A<-as.vector(rep(c(rep(1,5), rep(2,5), rep(3,5), rep(4,5), rep > (5,5)),2)) > B<-as.vector(rep(c(rep(c(1,2,3,4,5),5)),2)) > y<-as.vector(c > (15.5,15.2,14.2,14.3,15.8,6.2,7.2,6.6,6.2,5.6,15.4,13.9,13.4,12.5,13.2, > 10.9,12.5,12.3,11.0,12.3,7.5,6.7,7.2,7.6,6.3,14.9,15.2,14.2,14.3,16.4,7,8.4, > 7.8,7.6,7.4,14.4,13.3,14.8,14.1,15,11.3,12.7,11.7,12,13.3,6.7,7.3,6,7.6,7.1)) > lm1<-lm(y~factor(A)/factor(B)) > anova(lm1) > > Thanks a lot! And have a good weekend! > Regards, > WenqingFirst, the use of as.vector() above is redundant:> A <- as.vector(rep(c(rep(1,5), rep(2,5), rep(3,5), rep(4,5),rep(5,5)),2))> str(A)num [1:50] 1 1 1 1 1 2 2 2 2 2 ...> A1 <- rep(c(rep(1,5), rep(2,5), rep(3,5), rep(4,5), rep(5,5)),2)> str(A1)num [1:50] 1 1 1 1 1 2 2 2 2 2 ...> all.equal(A, A1)[1] TRUE On your question, a couple of options: # First create a data frame df <- data.frame(factor(A), factor(B), y) library(nlme) # Set your factors as the nested random effects lme1 <- lme(y ~ 1, random = ~ 1 | A / B, data = df)> summary(lme1)Linear mixed-effects model fit by REML Data: df AIC BIC logLik 147.4539 155.0211 -69.72693 Random effects: Formula: ~1 | A (Intercept) StdDev: 3.797784 Formula: ~1 | B %in% A (Intercept) Residual StdDev: 0.3768259 0.6957023 Fixed effects: y ~ 1 Value Std.Error DF t-value p-value (Intercept) 11 1.702937 25 6.45943 0 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -1.7696293 -0.7042397 0.1521976 0.5708701 1.5679386 Number of Observations: 50 Number of Groups: A B %in% A 5 25 Also: library(ape)> varcomp(lme1)A B Within 14.4231663 0.1419978 0.4840016 attr(,"class") [1] "varcomp" Note in both cases, lme() is used, not lm(). These are referenced in Section 10.2 (pg 279) of MASS4 by Venables & Ripley and in Section 4.2.3 (pg 167) of MEMSS by Pinheiro and Bates. HTH, Marc Schwartz