# 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