Divaker,
Thanks for the data.
For me,
> summary(aov(Purity~Supplier/Batch, process))
gives exactly the same results for mean squares as
> aov(Purity~Supplier+Error(Supplier/Batch), process)
except that the latter gives no p-values (because Supplier appears as
both error term and fixed effect, there isn't anything left to test
supplier against)
For compactness, I'll use the mean squares in
summary(aov(Purity~Supplier/Batch,data=process))
Df Sum Sq Mean Sq F value Pr(>F)
Supplier 2 15.056 7.528 2.8526 0.07736 .
Supplier:Batch 9 69.917 7.769 2.9439 0.01667 *
Residuals 24 63.333 2.639
The component of variance for batch is
(7.769-2.639)/3 = 1.71(1.307 as SD)
and for supplier, the between-supplier component of variance comes out
negative because it's 7.53-7.77 etc... and if you were testing the
supplier as a fixed effect, testing against the batch MS would imply
that supplier is not a significant effect.
For a mixed-effects model, I normally expected lme's syntax to be
lme(Purity~Supplier, random=~1|Batch, data=P),
but in this case this does not treat the batch as nested, because the
same batch levels appear in all Suppliers. I get the classical result
from
> process$BS<-factor(process$Supplier:process$Batch)
> lme(Purity~Supplier, random=~1|BS, data=process)
Linear mixed-effects model fit by REML
Data: P
Log-restricted-likelihood: -71.42198
Fixed: Purity ~ Supplier
(Intercept) SupplierT2 SupplierT3
-0.4166667 0.7500000 1.5833333
Random effects:
Formula: ~1 | BS
(Intercept) Residual
StdDev: 1.307622 1.624466
Which returns the expected 1.307 between-batch SD.
Lesson, subject to correction from the gods of lme: lme's random
effects are not treated as nested within fixed effects groups unless the
level identifiers for the random effects are unique to the fixed effects
group they are in.
And I obviously need to read Pinheiro and Bates again, or I'd have
known that without thinking about it.
Steve E
______________________________________________
R-help@r-project.org 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.
