[R] question about custom contrasts in ANOVA

Spencer Graves spencer.graves at pdf.com
Tue Aug 30 19:09:15 CEST 2005


Dear Prof. Ripley:

	  Agreed.  I thought that coding the factor levels manually and / or 
looking at model.matrix(y~xma) after further simplifying the problem 
might help Scot clarify his question.

	  Thanks,
	  Spencer Graves

Prof Brian Ripley wrote:

> I was puzzled as to what the question actually was. If you set
> 
> 
>>options(contrasts=c("contr.sum", "contr.poly"))
> 
> 
> the interaction contrasts are precisely those created manually.  Nothing 
> fancier is required.  But I am not sure what you want to do with them once 
> you have them.
> 
> On Tue, 30 Aug 2005, Spencer Graves wrote:
> 
> 
>>	  Since I have not seen a reply to this post, I will attempt a brief
>>comment:  I won't comment on the specifics, but your general approach
>>seems appropriate.  You already seem to know that a factor with k levels
>>is converted into (k-1) separate numerical variables, and a separate
>>regression coefficient is estimated for each one.
>>
>>	  To see in more detail what you were estimating, I looked at the
>>following:
>>
>>	  model.matrix(y~xma)
>>	  fit2 <- aov(y~xma)
>>	  attributes(fit2)
>>
>>	  From the latter, I identified "contrasts" as something that might be
>>interesting to examine, as follows:
>>
>>	  fit2$contrasts
>>
>>	  For clarity, I think I might reduce the number of observations
>>substantially, limiting myself to only 2 or 3 schools and paramterize
>>the problem manually, but preserving imbalance.  Then I'd use "lm",
>>specifying the terms in different orders.  With imbalance, the answer
>>depends on the order unfortunately.  When in doubt, I often experiment
>>with changing the order:  If the changes do not affect the conclusions,
>>I pick the simplest case to present to my audience.  If the changes do
>>affect the conclusions, I know I need to worry about which answer seems
>>most correct, and I also know something about the limits of the
>>conclusions.
>>
>>	  I know this doesn't answer your question, but I hope it helped with a
>>solution methodology.
>>
>>	  Best Wishes,
>>	  Spencer Graves
>>
>>Scot W McNary wrote:
>>
>>
>>>Hi,
>>>
>>>I have a problem in which I have test score data on students from a number
>>>of schools.  In each school I have a measure of whether or not they
>>>received special programming.  I am interested in the interaction between
>>>school and attendance to the programming, but in a very select set of
>>>comparisons.  I'd like to cast the test as one in which students in each
>>>school who attend are compared with students who don't across all schools.
>>>So, I would be comparing school 1 attenders with school 1 non-attenders,
>>>school 2 attenders with school 2 non-attenders, etc.  The reason for the
>>>custom contrast is that the between school comparisons (e.g., school 1
>>>attenders vs. school 2 non-attenders) are of less interest.
>>>
>>>This seems to require a custom contrast statement for the interaction
>>>term.  I have a toy example that seems to work as it should, but wonder if
>>>I've correctly created the contrast needed.
>>>
>>>Here is a toy example (code put together from bits taken from MASS ch 6,
>>>and various R-help postings, (e.g.,
>>>http://finzi.psych.upenn.edu/R/Rhelp02a/archive/49077.html)):
>>>
>>># toy interaction contrast example, 10 schools, 100 kids, 5 attenders (1)
>>># and 5 non-attenders (2) in each school
>>>
>>># make the data
>>>school <- gl(10, 10)
>>>attend <- gl(2, 5, 100)
>>># creates an interaction with schools 6 and 7
>>>y <- c(sample(seq(450, 650, 1), 50), rep(c(rep(650, 5), rep(450, 5)), 2),
>>> 	sample(seq(450, 650, 1), 30))
>>>
>>># anova
>>>summary(aov(y ~ school * attend))
>>>
>>># graphically
>>>Means <- tapply(y, list(school, attend), mean)
>>>
>>>plot(Means[,1], col="red", type = "l", ylim = c(400,700))
>>>
>>>points(Means[,2], col="blue", type = "l")
>>>
>>># create contrasts for hypothesis of interest
>>># school i attend j - school i attend j'
>>># for all schools
>>>sxa <- interaction(school, attend)
>>>sxam <- as.matrix(rbind(diag(1,10), diag(1,10) * -1))
>>>contrasts(sxa) <- sxam
>>>
>>>summary(aov(y ~ sxa), split=list(sxa=1:10), expand.split = T)
>>>
>>>The actual problem has a few more schools, other covariates, considerably
>>>more students, and is somewhat unbalanced.
>>>
>>>Thanks,
>>>
>>>Scot
>>>
>>>
>>>--
>>>   Scot W. McNary  email:smcnary at charm.net
>>>
>>>______________________________________________
>>>R-help at stat.math.ethz.ch mailing list
>>>https://stat.ethz.ch/mailman/listinfo/r-help
>>>PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>>
>>-- 
>>Spencer Graves, PhD
>>Senior Development Engineer
>>PDF Solutions, Inc.
>>333 West San Carlos Street Suite 700
>>San Jose, CA 95110, USA
>>
>>spencer.graves at pdf.com
>>www.pdf.com <http://www.pdf.com>
>>Tel:  408-938-4420
>>Fax: 408-280-7915
>>
>>______________________________________________
>>R-help at stat.math.ethz.ch mailing list
>>https://stat.ethz.ch/mailman/listinfo/r-help
>>PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>>
> 
> 

-- 
Spencer Graves, PhD
Senior Development Engineer
PDF Solutions, Inc.
333 West San Carlos Street Suite 700
San Jose, CA 95110, USA

spencer.graves at pdf.com
www.pdf.com <http://www.pdf.com>
Tel:  408-938-4420
Fax: 408-280-7915




More information about the R-help mailing list