Re: [R] Components of variance

From: Spencer Graves <spencer.graves_at_pdf.com>
Date: Mon 27 Jun 2005 - 02:10:52 EST

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@reSumry$lot@stdDev))
[1] 0.4356044 3.4763815
> 0.98007*sqrt(exp(qnorm(c(0.025, 0.975))*vcFit@reSumry$wafer@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@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 Mon Jun 27 02:14:57 2005

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