Re: [R] Checking for orthogonal contrasts
On Dec 3, 2010, at 17:17 , S Ellison wrote:
> 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.
You'll find that the block and treatment terms are NOT orthogonal. That's where
all the stuff about "efficiency factors" and "recovery of interblock
information" comes from.
>
> 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... ?
a) basically, yes, I think you do
b) yes, many, but there is an amazing amount of sloppily thought out "folklore"
going around, including the common misconception that somehow sum-to-zero
contrasts are inherently better than the other types. What does seem to be the
case is just that they have computational advantages in completely balanced
designs, because they then imply orthogonality of COLUMNS of the design matrix.
That in turn means that you can construct the sum of squares for each model
term based on its own columns only. In unbalanced designs, they just tend to
give incorrect results...
--
Peter Dalgaard
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd@cbs.dk Priv: pda...@gmail.com
__
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.
Re: [R] Checking for orthogonal contrasts
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.
Re: [R] Checking for orthogonal contrasts
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 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}}
__
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.
Re: [R] Checking for orthogonal contrasts
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 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}}
__
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.
[R] 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}}
__
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.
