# Re: [R] random effects with lmer() and lme(), three random factors

From: Spencer Graves <spencer.graves_at_pdf.com>
Date: Thu 03 Aug 2006 - 21:39:41 EST

In theory, 'lme' may be able to handle crossed random effects, but I don't know how to do it. I've used 'lmer' for such models, and your 'lmer' code looks plausible to me. Therefore, I would be inclined to believe the 'lmer' results.

To remove any lingering doubt, I suggest you construct Monte Carlo simulations where you know the structure. Then feed the simulations to 'lme' and 'lmer' and see what you get.

Hope this helps.
Spencer Graves

Xianqun (Wilson) Wang wrote:
> Hi, all,
>
>
>
> I have a question about random effects model. I am dealing with a
> three-factor experiment dataset. The response variable y is modeled
> against three factors: Samples, Operators, and Runs. The experimental
> design is as follow:
>
>
>
> 4 samples were randomly chosen from a large pool of test samples. Each
> of the 4 samples was analyzed by 4 operators, randomly selected from a
> group of operators. Each operator independently analyzed same samples
> over 5 runs (runs nested in operator). I would like to know the
> following things:
>
>
>
> (1) the standard deviation within each run;
>
> (2) the standard deviation between runs;
>
> (3) the standard deviation within operator
>
> (4) the standard deviation between operator.
>
>
>
> With this data, I assumed the three factors are all random effects. So
> the model I am looking for is
>
>
>
> Model: y = Samples(random) + Operator(random) + Operator:Run(random) +
> Error(Operator) + Error(Operator:Run) + Residuals
>
>
>
> I am using lme function in nlme package. Here is the R code I have
>
>
>
> 1. using lme:
>
> First I created a grouped data using
>
> gx <- groupedData(y ~ 1 | Sample, data=x)
>
> gx\$dummy <- factor(rep(1,nrow(gx)))
>
>
>
> then I run the lme
>
>
>
> fm<- lme(y ~ 1, data=gx,
> random=list(dummy=pdBlocked(list(pdIdent(~Sample-1),
>
> pdIdent(~Operator-1),
>
> pdIdent(~Operator:Run-1)))))
>
>
>
> finally, I use VarCorr to extract the variance components
>
>
>
> vc <- VarCorr(fm)
>
>
>
> Variance StdDev
>
> Operator:Run 1.595713e-10:20 1.263215e-05:20
>
> Sample 5.035235e+00: 4 2.243933e+00: 4
>
> Operator 5.483145e-04: 4 2.341612e-02: 4
>
> Residuals 8.543601e-02: 1 2.922944e-01: 1
>
>
>
>
>
> 2. Using lmer in Matrix package:
>
>
>
> fm <- lmer(y ~ (1 | Sample) + (1 | Operator) +
>
> (1|Operator:Run), data=x)
>
> summary(fm)
>
>
>
> Linear mixed-effects model fit by REML
>
> Formula: H.I.Index ~ (1 | Sample.Name) + (1 | Operator) + (1 |
> Operator:Run)
>
> Data: x
>
> AIC BIC logLik MLdeviance REMLdeviance
>
> 96.73522 109.0108 -44.36761 90.80064 88.73522
>
> Random effects:
>
> Groups Name Variance Std.Dev.
>
> Operator:Run (Intercept) 4.2718e-11 6.5359e-06
>
> Operator (Intercept) 5.4821e-04 2.3414e-02
>
> Sample (Intercept) 5.0352e+00 2.2439e+00
>
> Residual 8.5436e-02 2.9229e-01
>
> number of obs: 159, groups: Operator:Run, 20; Operator, 4; Sample.Name,
> 4
>
>
>
> Fixed effects:
>
> Estimate Std. Error t value
>
> (Intercept) 0.0020818 1.1222683 0.001855
>
>
>
>
>
> There is a difference between lmer and lme is for the factor
> Operator:Run. I cannot find where the problem is. Could anyone point me
> out if my model specification is correct for the problem I am dealing
> with. I am pretty new user to lme and lmer. Thanks for your help!
>
>
>
>
>
> Wilson Wang
>
>
>
>
>
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> and provide commented, minimal, self-contained, reproducible code.

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 and provide commented, minimal, self-contained, reproducible code. Received on Thu Aug 03 21:46:10 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Thu 03 Aug 2006 - 22:17:31 EST.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.