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

From: Xianqun (Wilson) Wang <>
Date: Tue 01 Aug 2006 - 05:27:39 EST

Dr. Bates,

Thanks for the notes! It helps. So now I see consistent resluts from both lme and lmer. Since I have several response variables to look at, I will reduce the model separately.

Speaking of the model reduction, it is clear in this example that the trivial variance of Operator:Run could be ignored. For a non-trivial-variance multivariate model, I wonder if there is any function (like step function for lm and glm) would work with lme or lmer and allow me to do the AIC-based model selection.


-----Original Message-----
From: [] On Behalf Of Douglas Bates
Sent: Saturday, July 29, 2006 7:51 AM
To: Xianqun (Wilson) Wang
Subject: Re: [R] random effects with lmer() and lme(), three random factors

On 7/28/06, 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)

That model specification can now be written as

fm <- lmer(y ~ (1|Sample) + (1|Operator/Run), 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.

It's just a matter of round-off. It is possible for the ML or REML estimates of a variance component to be zero, as is the case here, but the current computational methods do not allow the value zero because this will cause some of the matrix decompositions to fail. In lmer we use a constrained optimization with the relative variance (variance of a random effect divided by the residual variance) constrained to be greater than or equal to 5e-10, which is exactly the value you have here.

I'll add code to the model fitting routine to warn the user when convergence to the boundary value occurs. I haven't done that in the past because it is not always easy to explain what is occurring. For a model with variance components only, like yours, convergence on the boundary means that an estimated variance component is zero. In the case of bivariate or multivariate random effects the variance-covariance matrix can be singular without either of the variances being zero.

The bottom line for you is that the estimated variance for Operator:Run is zero and you should reduce the model to y ~ (1|Sample) + (1|Operator)

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]]
> ______________________________________________
> mailing list
> PLEASE do read the posting guide
> and provide commented, minimal, self-contained, reproducible code.
> mailing list PLEASE do read the posting guide and provide commented, minimal, self-contained, reproducible code. Received on Tue Aug 01 08:43:54 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 Tue 01 Aug 2006 - 10:18:25 EST.

Mailing list information is available at Please read the posting guide before posting to the list.