[R] problem with se.contrast()

[Prof Brian Ripley]

On Mon, 21 Feb 2005, Peter Dalgaard wrote:

> Prof Brian Ripley <ripley at stats.ox.ac.uk> writes:
>
>>>> test.aov <- with(testdata,aov(Measurement ~ Material + Error(Lab/Material)))
>>>> se.contrast(test.aov,
>>>>             list(Material=="A",Material=="B",Material=="C",Material=="D"),
>>>>             coef=c(0.5,0.5,-0.5,-0.5),data=testdata)
>>>> [1] 0.1432572
>>>
>>> I got a different result and I have admit that I didn't understand why
>>> there is a differnce between the lme model and this one. There are some
>>> comments in the help pages but I'm not sure if this is the answer.
>>
>> It is.  You used the default `contrasts', which are not actually
>> contrasts.  With
>>
>> options(contrasts=c("contr.helmert", "contr.poly"))
>>
>> it gives the same answer as the other two.  Because you used
>> non-contrasts there was an efficiency loss (to the Intercept stratum).
>
> Brian,
>
> I'm not sure how useful that contrasts-that-are-not-contrasts line is.

I agree, it was not precise enough (too late at night for me).  Try 
`non-orthogonal contrasts'.  The issue was correct though, it is the 
choice of contrasts, and as I would automatically use orthogonal contrasts 
with aov() I had not encountered it and it took me a while to pick up on 
what Christoph had done differently from me (I had run the example and got 
the same result as the randomized-block analysis before my original post).

There's a comment in the code for aov:

         ##  helmert contrasts can be helpful: do we want to force them?
         ##  this version does for the Error model.

and perhaps we should make them the default for the per-stratum fits.

> It certainly depends on your definition of "contrasts". Contrast
> matrices having zero column sums was not part of the definition I was
> taught. I have contrasts as "representations of the group mean
> structure that are invariant to changes of the overall level", so
> treatment contrasts are perfectly good contrasts in my book.

I don't think Yates thought in those terms, and the whole idea of dividing 
into strata (and the generalized Yates algorithm) is based on contrasts 
being orthogonal to the overall mean (and to things in other strata).

It was that, not zero-sum, that I was taught, but in balanced cases such 
as here it is the same thing.

> The zero-sum condition strikes me as a bit arbitrary: after all there
> are perfectly nice orthogonal designs where some levels of a factor
> occur more frequently than others.

Balance and orthogonality are not quite the same, though.

> This in turn makes me a bit wary of what is going on inside 
> se.contrasts, but it's gotten too late for me to actually study the code 
> tonight.
>
> Could you elaborate on where precisely the condition on the contrast
> matrices comes into play?

In finding the projection lengths in eff.aovlist, here

     proj.len <-
 	lapply(aovlist, function(x)
 	   {
 	       asgn <- x$assign[x$qr$pivot[1:x$rank]]
 	       sp <- split(seq(along=asgn), attr(terms(x), "term.labels")[asgn])
 	       sapply(sp, function(x, y) sum(y[x]), y=diag(x$qr$qr)^2)
 	   })

using only the diagonal requires orthogonality of the coding.

It is many years since I looked at this, and it is not immediately clear 
to me how best to calculate efficiencies without this assumption (which 
could be tested inside eff.aovlist, of course).

[People wondering why this was ever useful given lme should be aware that
1) It predates the availability of lme.
2) When I first used this in 1998, lme appeared very slow.
3) People with a classical aov training (not me, but e.g. Bill Venables) 
think in these terms.]

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595