Susanne Schwenke <ml-r-help at epigenomics.com> writes:
Dear R-help group,
I have a two-way ANOVA with two crossed random factors, no
nesting. Each factor has three levels, resulting in 9 cells for the
experiment. Each cell contains 10 repetitions. According to the
ANOVA model I assume equal variances for all levels per factor.
I would like to get REML-estimates for the variances of the two
factors and moreover get confidence intervals for these estimates,
so the use of the nlme-package seems to be a good idea.
My problem in the first place is to formulate the model itself for
the lme-function. The fixed part would at most consist of the
intercept, resulting in
fixed= response ~ 1
and the random part would be
random = ~ a + b
but I have no idea what my gouping factor there should be. Could
somebody please point me in the right direction ?
Sorry if this turns out to be an extremely simple question, I'm a
newbie to R ...
Many greetings,
Susanne
----
Susanne Schwenke
Epigenomics AG www.epigenomics.com Kastanienallee 24
+4930243450 10435 Berlin
The answer to your question is not "extremely simple". It happens
that the lme function is much better suited to nested random effects
than to crossed random effects. To estimate crossed random effects
you must create an awkward formulation with a grouping factor that has
one level and the random effects model matrix based on the indicators
for factor a and the indicators for factor b. These two sets of
random effects each have variance-covariance matrices that are
multiples of an identity and the are grouped together as a
block-diagonal matrix. The whole formulation is
lme(response ~ 1, data = myData,
random = pdBlocked(list(pdIdent(~ a - 1), pdIdent(~ b - 1))))
and myData must be a groupedData object with a grouping factor that
has only one level.
There is an example of fitting crossed random effects in section 4.2
of Pinheiro and Bates (2000) "Mixed-effects Models in S and S-PLUS",
Springer. That example happens to have two blocks and a random effect
for the block is added. If we fit a single block it would look like