[R] Checking for orthogonal contrasts
S Ellison
S.Ellison at lgc.co.uk
Fri Dec 3 17:17:58 CET 2010
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}}
More information about the R-help
mailing list