[R] Checking for orthogonal contrasts

[S Ellison]

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}}