[R] Strata and Degrees of freedom in anova and multi-level modeling

From: John Maindonald <john.maindonald_at_anu.edu.au>
Date: Thu 16 Feb 2006 - 12:04:43 EST

I am changing the title because this is really about the history of anova,
and about strata in analysis of variance. As this kind of question has been
arising very frequently, an extended comment may be in order.

The ideas, and the sums of squares breakdowns, go back to Fisher; see in particular his "Design of Experiments", first published in 1935. This
book is still a good read. I'm not sure that Fisher talked about error strata,
but he certainly did emphasize the need to choose the error sum of squares
that was appropriate to the inference that was of interest.

It is unfortunate that, in later developments, analysis of variance was often
presented in a manner that forgot important parts of Fisher's insights. One
of my complaints about the data mining community is that it seems remarkably resistant to the notion that the distribution of the error term in a
model can matter a lot.

Thus in an experiment that has blocks, plots within blocks and subplots (I'll use as an example the kiwishade data set in the DAAG package), treatments that are assigned at random to subplots are compared using the
within plot residual (this ignores the possibility of a plot:subplotTreatment
interaction, so that such designs must be used with care and/or finesse),
treatments that are applied to whole plots are compared using the within block residual, and so on.

In the kiwishade experiment, there were three blocks, four plots per block
(one for each of 4 treatments) and four subplots (vines) per treatment. See Maindonald and Braun (the book we've been calling DAAGUR), pp.230-239, for a detailed exposition of the analysis of these data, comparing the aov() output with the lme() output. [Details of this and other
R-related books can be obtained from a CRAN site.]

For this experiment, the error structure is hierarchical, and it thus makes sense
to equate strata with levels, in the sense in which this word is used in the
documentation for lme(). The strata are: within plots [level 0 in lme ()], within
blocks between plots (level 1), and between blocks (level 2). To see the
aov() output, type:

   aov(yield ~ shade + Error(block/plot), data=kiwishade) [Blocks might alternatively be regarded as a fixed effect.]

To understand the use of the between plot sum of squares, observe that there
are 12 independent pieces of information at the plot level, i.e. 12 plots.
However, 2 (3-1) of these can be explained as due to differences between blocks,
and 3 of them can be explained as due to differences between treatments. Thus, we are left with 7 independent (in an algebraic sense) items of information.
The corresponding error variance [the usual sum of squares about the mean
divided by number of values less 1] has 6 degrees of freedom.

The within plot residual mean square has 36 degrees of freedom. This experiment
did not have treatment that was applied at the subplot level, and this is not the
correct mean squares for any treatment comparison. The analysis is readily
simplified by calculating plot means, and basing all analyses on those. In just the
same way:
1) treatments may be applied to whole apples, there may be multiple measurements on an apple, and a reasonable way to do the analysis may be to calculate means for each apple.
2) Commonly, in laboratory experiments, treatments are applied to whole samples, and multiple measurements or assays are made on the one sample. 3) In a comparison of prescribing habits between registrars and specialists,
there might be data from 50 registrars and 50 specialists in each hospital,
but perhaps only 3 hospitals. So, do we have

   (a) 3 pieces of information (3 hospitals with two items to compare per hospital;

    there'll be 2 degrees of freedom for the mean square), or     (b) closer to 300 [a split plot type comparison will give a mean square with

    (50+50-2)*3 degrees of freedom].

John Maindonald email: john.maindonald@anu.edu.au phone : +61 2 (6125)3473 fax : +61 2(6125)5549 Mathematical Sciences Institute, Room 1194, John Dedman Mathematical Sciences Building (Building 27) Australian National University, Canberra ACT 0200.

Genstat's achievement was to push to the limit, based on ideas of JA Nelder
(a theory) and GN Wilkinson (an algorithm), ideas of balance that made it
possible to retain the conceptual simplicity of analysis of variance type strata.
Estimates of effects of interest may for crossed designs use an error term that
combines variances across strata, so that "t"-statistics (with a fudge used to get
degrees of freedom) may at best be reasonably approximated by a t- distribution.

The methodology that is used in lme(), and other software that uses a similar
approach, gives up the attempt to form an analysis of variance table. The
model specifies the variance-covariance structure, and an important part of the
analysis is the estimation of that variance-covariance structure. Estimates of
effects that are of interest no longer have a t-distribution under the standard
normality etc assumptions. Standard fudges (calculate a df for a t- distribution
approximation) often work quite well. Doug Bates is right to keep reminding
everyone that they are fudges.

The fudges are important for another more pervasive reason. It is important
to know whether a variance is estimated with 2 or with ~300 degrees of freedom. If there are only 2 pieces of information the normality assumptions
matter a lot. Was one of those hospitals an outlier? Problems with the fudge
used for calculate degrees of freedom for a t-approximation may be the least
of the matters that should be a concern.

There's a nice account of Fisher's contribution to experimental design (as well as a few complaints about some of Fisher's comments on inference) in:
R. A. Fisher and Experimental Design: A Review D. A. Preece
Biometrics, Vol. 46, No. 4 (Dec., 1990) , pp. 925-935

John Maindonald

On 15 Feb 2006, at 10:00 PM, r-help-request@stat.math.ethz.ch wrote:

> From: Peter Dalgaard <p.dalgaard@biostat.ku.dk>
> Date: 15 February 2006 3:26:27 AM
> To: WPhantom <wp1@tiscali.fr>
> Cc: r-help@stat.math.ethz.ch
> Subject: Re: [R] A concrete type I/III Sum of square problem
> WPhantom <wp1@tiscali.fr> writes:
>> Thanks Brian for the reference.
>> I just discover that it is available in our
>> library so I going to take it & read it soon.
>> Actually, I don't even know the difference
>> between a multistratum vs a single-stratum AOV. A
>> quick search on google returned me the R materials so that I imagine
>> that these concepts are quite specific to R.
> You have to be careful not to confuse Google's view of the world with
> Reality...
> The concept of error strata is much older than R, and existed for
> instance in Genstat, anno 1977 or so. However, Genstat seems to have
> left little impression on the Internet.
>> I will read the book first before asking for more informations.
> The executive summary is that the concept of error strata relies
> substantially on having a balanced design (at least for the random
> effects), so that the analysis can be decomposed into analyses of
> means, contrasts, and contrasts of means. For unbalanced designs, you
> usually get meaningless analyses.
>> Thanks
>> Sylvain Clément
>> At 12:38 14/02/2006, you wrote:
>>> More to the point, you are confusing
>>> multistratum AOV with single-stratuam AOV. For
>>> a good tutorial, see MASS4 (bibliographic
>>> information in the R FAQ). For unbalanced data
>>> we suggest you use lme() instead.
>>> --
>>> Brian D. Ripley, ripley@stats.ox.ac.uk
>> ______________________________________________
>> 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
> --
> O__ ---- Peter Dalgaard Øster Farimagsgade 5, Entr.B
> c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
> (*) \(*) -- University of Copenhagen Denmark Ph: (+45)
> 35327918
> ~~~~~~~~~~ - (p.dalgaard@biostat.ku.dk) FAX: (+45)
> 35327907

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 Received on Thu Feb 16 12:09:47 2006

This archive was generated by hypermail 2.1.8 : Thu 16 Feb 2006 - 16:35:11 EST