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

[John Maindonald]

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 at 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 at stat.math.ethz.ch wrote:

> From: Peter Dalgaard <p.dalgaard at biostat.ku.dk>
> Date: 15 February 2006 3:26:27 AM
> To: WPhantom <wp1 at tiscali.fr>
> Cc: r-help at stat.math.ethz.ch
> Subject: Re: [R] A concrete type I/III Sum of square problem
>
>
> WPhantom <wp1 at 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 at stats.ox.ac.uk
>>
>> ______________________________________________
>> R-help at 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 at biostat.ku.dk)                  FAX: (+45)  
> 35327907