[R] random effects with lmer() and lme(), three random factors
Spencer Graves
spencer.graves at pdf.com
Thu Aug 3 13:39:41 CEST 2006
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 at 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.
More information about the R-help
mailing list