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
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
