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

[Spencer Graves]

	  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.