[R] help on extract variance components from the fitted model by lm

[Marc Schwartz]

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,
> Wenqing

First, 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