[R] Components of variance

Spencer Graves spencer.graves at pdf.com
Sun Jun 26 18:10:52 CEST 2005 

Hi, Renaud:

	  Thanks for the link to Doug Bates' recent comment on this.
Unfortunately, the information contained therein is not enough for me to
easily obtain all the required standard errors and write code to make
confidence intervals.  To be more specific, I'm able to produce the same
answers as before using lmer in place of lme as follows:

> library(lme4)
> 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))
> (fitr <- lmer(Vt~1+(1|lot)+(1|wafer), DF))
<snip>
Random effects:
  Groups   Name        Variance Std.Dev.
  wafer    (Intercept) 0.96054  0.98007
  lot      (Intercept) 1.51434  1.23058
  Residual             1.34847  1.16124
# of obs: 27, groups: wafer, 9; lot, 3

Fixed effects:
             Estimate Std. Error DF t value Pr(>|t|)
(Intercept)  0.60839    0.81329 26  0.7481   0.4611

	  These numbers agree with what I got with lme.  I then followed your
suggetion to "http://tolstoy.newcastle.edu.au/R/help/05/06/7291.html".
This suggested I try str(VarCorr(fitr)):

> vcFit <- VarCorr(fitr)
> str(vcFit)
Formal class 'VarCorr' [package "Matrix"] with 3 slots
   ..@ scale   : num 1.16
   ..@ reSumry :List of 2
   .. ..$ wafer:Formal class 'corrmatrix' [package "Matrix"] with 2 slots
   .. .. .. ..@ .Data : num [1, 1] 1
   .. .. .. .. ..- attr(*, "dimnames")=List of 2
   .. .. .. .. .. ..$ : chr "(Intercept)"
   .. .. .. .. .. ..$ : chr "(Intercept)"
   .. .. .. ..@ stdDev: Named num 0.844
   .. .. .. .. ..- attr(*, "names")= chr "(Intercept)"
   .. ..$ lot  :Formal class 'corrmatrix' [package "Matrix"] with 2 slots
   .. .. .. ..@ .Data : num [1, 1] 1
   .. .. .. .. ..- attr(*, "dimnames")=List of 2
   .. .. .. .. .. ..$ : chr "(Intercept)"
   .. .. .. .. .. ..$ : chr "(Intercept)"
   .. .. .. ..@ stdDev: Named num 1.06
   .. .. .. .. ..- attr(*, "names")= chr "(Intercept)"
   ..@ useScale: logi TRUE
>

	  Unfortunately, I've been so far unable to translate this into 
confidence intervals for the variance components that seem to match 
intervals(lme(...)).  The closest I came was as follows:
	
> 1.23058*sqrt(exp(qnorm(c(0.025, 0.975))*vcFit at reSumry$lot at stdDev))
[1] 0.4356044 3.4763815
> 0.98007*sqrt(exp(qnorm(c(0.025, 0.975))*vcFit at reSumry$wafer at stdDev))
[1] 0.4286029 2.2410887
>
	  The discrepancy between these numbers and intervals(lme(...)) is 25% 
for "lot", which suggest that I'm doing something wrong.  Moreover, I 
don't know how to get a similar interval for "scale", and I don't know 
how to access the estimated variance components other than manually.

	  Comments?
	  Thanks,
	  spencer graves
p.s.  R 2.1.0 patched under Windows XP.

Renaud Lancelot wrote:

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

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

More information about the R-help mailing list