[R] Components of variance

[Renaud Lancelot]

Spencer Graves a écrit :
> 	  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.

The issue of variance components in lme4 was recently by D. Bates:

http://tolstoy.newcastle.edu.au/R/help/05/06/7291.html

Also, a colleague wrote (in French) a nice report - with data and R code 
on how to compute variance components, their bias and distribution with 
lme (normal approximation , Monte Carlo simulation and delta method). 
Let me know if you want it.

Best,

Renaud


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


-- 
Dr Renaud Lancelot, vétérinaire
Projet FSP régional épidémiologie vétérinaire
C/0 Ambassade de France - SCAC
BP 834 Antananarivo 101 - Madagascar

e-mail: renaud.lancelot at cirad.fr
tel.:   +261 32 40 165 53 (cell)
         +261 20 22 665 36 ext. 225 (work)
         +261 20 22 494 37 (home)