thr3ads.net - R help - [R] random effects with lmer() and lme(), three random factors [Jul 2006]

If this information is useful, please help other people find it:
Share via:

Xianqun (Wilson) Wang

2006-Jul-28 23:43 UTC

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

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]]

Xianqun (Wilson) Wang

2006-Jul-28 23:58 UTC

head link

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

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]]

Douglas Bates

2006-Jul-29 15:51 UTC

head link

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

On 7/28/06, Xianqun (Wilson) Wang <xwang at aviaradx.com>
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]]
>
> ______________________________________________
> 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.
>

Spencer Graves

2006-Aug-03 11:39 UTC

head link

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

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.

Reasonably Related Threads

Search for more reasonably related threads

R help - Jul 2006 - random effects with lmer() and lme(), three random factors

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

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

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

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

Reasonably Related Threads