From: Douglas Bates <dmbates_at_gmail.com>

Date: Mon 27 Jun 2005 - 04:40:36 EST

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 04:57:15 2005

Date: Mon 27 Jun 2005 - 04:40:36 EST

On 6/26/05, Spencer Graves <spencer.graves@pdf.com> wrote:

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

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 04:57:15 2005

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