[R] problem with se.contrast()

[Christoph Buser]

Dear Prof Ripley, Dear Prof Dalgaard 

Thank you both for your help. I tried it with helmert contrasts
and got a result that is consistent with lme. I didn't realize
that the parameterization of the model has an influence on the
contrasts that I tried to test. 
It seems that I should read a little bit more about this issue
for my better understanding.

I have one last point to propose:

You could include (as interim solution) a warning (that there might
be an in efficiency loss) in se.contrast() if one uses
non-orthogonal contrasts and a multi-stratum model.

Best regards,

Christoph Buser

--------------------------------------------------------------
Christoph Buser <buser at stat.math.ethz.ch>
Seminar fuer Statistik, LEO C11
ETH (Federal Inst. Technology)	8092 Zurich	 SWITZERLAND
phone: x-41-1-632-5414		fax: 632-1228
http://stat.ethz.ch/~buser/
--------------------------------------------------------------


Prof Brian Ripley writes:
 > 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