Re: [R] Components of variance

From: Spencer Graves <spencer.graves_at_pdf.com>
Date: Sun 26 Jun 2005 - 12:13:25 EST

          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:"

          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@grecc.umaryland.edu
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help@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@pdf.com
www.pdf.com <http://www.pdf.com>
Tel:  408-938-4420
Fax: 408-280-7915

______________________________________________
R-help@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
Received on Sun Jun 26 12:18:13 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:33:02 EST