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

______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Reply via email to