[R] Components of variance

[Spencer Graves]

	  As far as I know, the best available is lme in library(nlme).  For 
more information, see the the following:

	  Pinheiro and Bates (2000) Mixed-Effects Models in S and S-Plus 
(Springer)

	  Consider the following example:

 > set.seed(2)
 > lot <- rep(LETTERS[1:3], each=9)
 > lot.e <- rep(rnorm(3), each=9)
 > wf <- paste(lot, rep(1:9, each=3))
 > wf.e <- rep(rnorm(9), each=3)
 > DF <- data.frame(lot=lot, wafer=wf,
+ Vt=lot.e+wf.e+rnorm(27))
 > (fit <- lme(Vt~1, random=~1|lot/wafer, DF))
Linear mixed-effects model fit by REML
   Data: DF
   Log-restricted-likelihood: -48.44022
   Fixed: Vt ~ 1
(Intercept)
   0.6083933

Random effects:
  Formula: ~1 | lot
         (Intercept)
StdDev:    1.230572

  Formula: ~1 | wafer %in% lot
         (Intercept) Residual
StdDev:   0.9801403 1.161218

Number of Observations: 27
Number of Groups:
            lot wafer %in% lot
              3              9
 > (CI.fit <- intervals(fit))
Approximate 95% confidence intervals

  Fixed effects:
                 lower      est.    upper
(Intercept) -1.100281 0.6083933 2.317068
attr(,"label")
[1] "Fixed effects:"

  Random Effects:
   Level: lot
                     lower     est.    upper
sd((Intercept)) 0.3368174 1.230572 4.495931
   Level: wafer
                    lower      est.    upper
sd((Intercept)) 0.426171 0.9801403 2.254201

  Within-group standard error:
     lower      est.     upper
0.8378296 1.1612179 1.6094289
 > str(CI.fit)
List of 3
  $ fixed   : num [1, 1:3] -1.100  0.608  2.317
   ..- attr(*, "dimnames")=List of 2
   .. ..$ : chr "(Intercept)"
   .. ..$ : chr [1:3] "lower" "est." "upper"
   ..- attr(*, "label")= chr "Fixed effects:"
  $ reStruct:List of 2
   ..$ lot  :`data.frame':       1 obs. of  3 variables:
   .. ..$ lower: num 0.337
   .. ..$ est. : num 1.23
   .. ..$ upper: num 4.5
   ..$ wafer:`data.frame':       1 obs. of  3 variables:
   .. ..$ lower: num 0.426
   .. ..$ est. : num 0.98
   .. ..$ upper: num 2.25
   ..- attr(*, "label")= chr "Random Effects:"
  $ sigma   : atomic [1:3] 0.838 1.161 1.609
   ..- attr(*, "label")= chr "Within-group standard error:"
  - attr(*, "level")= num 0.95
  - attr(*, "class")= chr "intervals.lme"
 > diff(log(CI.fit$sigma))
    est.   upper
0.32641 0.32641

	  The last line combined with help for intervals.lme shows that the 
confidence interval for sigma (and doubtless also for lot and wafer 
variance components is based on a normal approximation for the 
distribution of log(sigma).

	  The state of the art is reflected in "lmer" in library(lme4), 
described in the following:

	  Doug Bates (2005) "Fitting linear mixed models in R" in R News 5/1 
available from "www.r-project.org" -> Newsletter

	  However, an "intervals" function is not yet available for "lmer" 
objects.
	
	  spencer graves

John Sorkin wrote:

> Could someone identify a function that I might use to perform a
> components of variance analysis? In addition to the variance
> attributable to each factor, I would also like to obtain the SE of the
> variances.
> Thank you,
> John
>  
> John Sorkin M.D., Ph.D.
> Chief, Biostatistics and Informatics
> Baltimore VA Medical Center GRECC and
> University of Maryland School of Medicine Claude Pepper OAIC
>  
> University of Maryland School of Medicine
> Division of Gerontology
> Baltimore VA Medical Center
> 10 North Greene Street
> GRECC (BT/18/GR)
> Baltimore, MD 21201-1524
>  
> 410-605-7119 
> -- NOTE NEW EMAIL ADDRESS:
> jsorkin at grecc.umaryland.edu
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

-- 
Spencer Graves, PhD
Senior Development Engineer
PDF Solutions, Inc.
333 West San Carlos Street Suite 700
San Jose, CA 95110, USA

spencer.graves at pdf.com
www.pdf.com <http://www.pdf.com>
Tel:  408-938-4420
Fax: 408-280-7915