R help - Dec 2010 - Checking for orthogonal contrasts

A common point made in discussion of contrasts, type I, II, III SS etc
is that for sensible comparisons one should use contrasts that are 
'orthogonal in the row-basis of the model matrix' (to quote from
http://finzi.psych.upenn.edu/R/Rhelp02/archive/111550.html)

Question: How would one check, in R, that this is so for a particular
fitted linear model object?

Steve Ellison



*******************************************************************
This email and any attachments are confidential. Any use...{{dropped:8}}

Hi Steve,
The short answer is that there is typically no reason to check (beyond
looking to see what contrasts have been defined for the factors in the
model; see ?contrasts) because the rules are pretty simple.
1) the design matrix is orthogonal on the row-basis (i.e., the columns
sum to zero) when using contr.sum, contr.helmert, or manually set
contrasts in which the contrasts sum to zero.
2) if the design matrix is orthogonal and cell sizes are equal, the
model matrix will also be orthogonal.

However, you can of course check, as in the following example:

dat <- data.frame(y=rnorm(100), x=factor(rep(c("a", "a",
"b", "c"), 25)))
m <- lm(y ~ x, data=dat)
## model matrix orthoganal on row-basis?
apply(model.matrix(m)[,grep("x", colnames(model.matrix(m)))], 2, sum)
## No, sum is non-zero
## design matrix orthoganal on row-basis?
apply(contrasts(dat$x), 2, sum) ## No, sum is non-zero
## change contrasts
contrasts(dat$x) <- contr.sum(n=3)
m2 <- lm(y ~ x, data=dat)
## model matrix orthoganal on row-basis?
apply(model.matrix(m2)[,grep("x", colnames(model.matrix(m2)))], 2,
sum) ## No, sum is non-zero (because of unbalenced cell sizes)
## model matrix orthoganal on row-basis?
apply(contrasts(dat$x), 2, sum) ## Yes, sum is zero.

HTH,
Ista

On Fri, Dec 3, 2010 at 8:47 AM, S Ellison <S.Ellison at lgc.co.uk>
wrote:
>
> A common point made in discussion of contrasts, type I, II, III SS etc
> is that for sensible comparisons one should use contrasts that are
> 'orthogonal in the row-basis of the model matrix' (to quote from
> http://finzi.psych.upenn.edu/R/Rhelp02/archive/111550.html)
>
> Question: How would one check, in R, that this is so for a particular
> fitted linear model object?
>
> Steve Ellison
>
>
>
> *******************************************************************
> This email and any attachments are confidential. Any u...{{dropped:18}}

There's a pretty good section in the R book by Crawley on contrast
statements in R, including some discussion of the contrasts being
orthogonal.

I would say you should just make your own table and sort it out there- if
you have equal sample sizes, then the contrast coefficients along the row
should sum to 0. Further, if the column sums also sum to zero, then you know
you are within your experimentwise error rate as well. Sokal and Rohlf
(1995) has a good discussion about this as well, as does Zar 1999 (I
think... can't be  sure).

Alternatively, there's some good stuff using hypothesis() in John Fox's
car() package for setting up appropriate contrast statements in linear
models.

Mike

On Fri, Dec 3, 2010 at 7:47 AM, S Ellison <S.Ellison@lgc.co.uk> wrote:
>
> A common point made in discussion of contrasts, type I, II, III SS etc
> is that for sensible comparisons one should use contrasts that are
> 'orthogonal in the row-basis of the model matrix' (to quote from
> http://finzi.psych.upenn.edu/R/Rhelp02/archive/111550.html)
>
> Question: How would one check, in R, that this is so for a particular
> fitted linear model object?
>
> Steve Ellison
>
>
>
> *******************************************************************
> This email and any attachments are confidential. Any u...{{dropped:20}}

David,

Thanks for the comments.

I think, though, that I have found the answer to my own post.
>Question: How would one check, in R, that [contrasts .. are 
>'orthogonal in the row-basis of the model matrix'] for a particular
>fitted linear model object?
?lm illustrates the use of crossprod() for probing the orthogonality of
a model matrix. If I understand correctly, the necessary condition is
essentially that all between-term off-diagonal elements of crossprod(m)
are zero if the contrasts are orthogonal, where 'term' refers to the
collection of columns related to a single term in the model formula.

Example:

y<-rnorm(27)
g <- gl(3, 9)
h <- gl(3,3,27)

m1 <- model.matrix(y~g*h,  contrasts = list(g="contr.sum",
h="contr.sum"))
crossprod(m1)

#Compare with
m2 <- model.matrix(y~g*h,  contrasts = list(g="contr.treatment",
h="contr.treatment"))
crossprod(m2)
	#Note the nonzero off-diagonal elements between, say, g and h or
g, h and the various gi:hj elements


That presumably implies that one could test a linear model explicitly
for contrast orthogonality (and, implicitly, balanced design?) using
something like

model.orthogonal.lm <- function(l) {
	#l is a linear model 
	m <- model.matrix(l)
	a <- attr(m, "assign")
	a.outer <- outer(a, a, FUN="!=")
	m.xprod <- crossprod(m) 
	all( m.xprod[a.outer] == 0 )
}

l1 <- lm(y~g*h,  contrasts = list(g="contr.sum",
h="contr.sum"))

l2 <- lm(y~g*h,  contrasts = list(g="contr.treatment",
h="contr.treatment"))

model.orthogonal.lm(l1) 
	#TRUE

model.orthogonal.lm(l2)
	#FALSE

Not sure how it would work on balanced incomplete block designs,
though. I'll have to try it.

Before I do, though, a) do I have the stats right? and b) this now
seems so obvious that someone must already have done it somewhere... ?


*******************************************************************
This email and any attachments are confidential. Any\ us...{{dropped:19}}