# This file is posted as pigment.ssc
pigment <- read.table("pigment.dat",
col.names=c("Batch","Sample","Test","Y"))
pigment$Batch <- as.factor(pigment$Batch)
pigment$Sample <- as.factor(pigment$Sample)
# The function raov() may be used for balanced
# designs when all factors are random. It gives
# a conventional analysis including the estimation
# of variance components.
# The function varcomp() is more general, and
# may be used to estimate variance components
# for balanced or unbalanced mixed models.
# raov(): Random Effects Analysis of Variance
pigment.raov <- raov(Y ~ Batch/Sample,
data=pigment)
summary(pigment.raov)
names(pigment.raov)
# List the number of replications for each
# level of the design.
pigment.raov$rep
# Print information on the expected mean
# squares
pigment.raov$ems.coef
# The same variance component estimates can be
# found using varcomp(), but as this allows mixed
# models we first must declare which facotrs are
# random using the is.random() function.
# All factors in the data frame are declared to
# be random effect by
is.random(pigment) <- T
is.random(pigment)
# The possible estimation methods are
# "minque0": minimum norm quadratic estimators
# (the default)
# "reml" : residual (or reduced or restricted)
maximum likelihood.
# "ml" : maximum likelihood.
varcomp(Y ~ Batch/Sample, data=pigment,
method="reml")$var
varcomp(Y ~ Batch/Sample, data=pigment,
method="ml")$var