From: Spencer Graves <spencer.graves_at_pdf.com>

Date: Mon 27 Jun 2005 - 02:10:52 EST

# of obs: 27, groups: wafer, 9; lot, 3

>

p.s. R 2.1.0 patched under Windows XP.

*>
*

*>
*

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

*>
*

*>
*

*>
*

*>
*

*>
*

Date: Mon 27 Jun 2005 - 02:10:52 EST

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

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

>> >> 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.htmlReceived 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
*