[R] How to compute contrast where there are interaction terms in the linear model?

[Charles C. Berry]

On Fri, 4 Dec 2009, Peng Yu wrote:

> On Fri, Dec 4, 2009 at 12:18 PM, Charles C. Berry <cberry at tajo.ucsd.edu> wrote:
>> On Fri, 4 Dec 2009, Peng Yu wrote:
>>
>>> On Tue, Dec 1, 2009 at 4:04 PM, Charles C. Berry <cberry at tajo.ucsd.edu>
>>> wrote:
>>>>
>>>> On Tue, 1 Dec 2009, Peng Yu wrote:
>>>>
>>>>> Could somebody recommend some textbook how to compute contrast when
>>>>> there are interactions terms? "Applied Linear Regression Models"
>>>>> (book) mentioned contrast, but I cannot extend it to the case where
>>>>> there are interaction terms.
>>>>
>>>>
>>>> Textbook? Schmextbook!
>>>>
>>>> You have the power of R at your fingertips!
>>>>
>>>> Use it to explore concepts you are trying to wrap your brain around!
>>>>
>>>>
>>>>> df <- expand.grid(x1=1:20,x2=factor(letters[1:2]))
>>>>> vanilla <- model.matrix(~0+poly(x1,degree=3), df )
>>>>> matplot(row( vanilla) , vanilla, type='b')
>>>>> inter <-  model.matrix(~0+poly(x1,degree=3):x2, df )
>>>>> matplot(row( inter ) , inter, type='b')
>>>
>>> I don't understand. Where the contrast is?
>
> I don't understand what polynomial contrasts are. Would you please
> help me understand?
>

Look at the graphs, please. They are linear basis functions for 
polynomials. Construct a third order polynomial in x <- 1:20 (say) then 
fit it with poly( x , degree=1), poly ( x, degree =2 ) , ..., poly( x, 
degree=19 ). Compare the predict() ed values of each to the actual 
polynomial. Study the coefficients.


> May I ask you the following question to understand contrasts in a
> simpler examples?
>
> In the following example, I fit formula Y~x1+x2 and Y~x1*x2.
>
> For Y~x1+x2, I understand that coefficients(afit)['x1b'] is the
> contrast between the factor mean of a and the factor mean of b.
>
> But for Y~x1*x2, since there are the interaction terms, I'm not sure
> if it makes sense to compare the levels of a and b. Would you please
> let me know what coefficients(afit)['x1b'] means in this case?
>

Well, it is literally the coefficient of

 	model.matrix(afit)[, 'x1b' ]

and you can see from inspection of that matrix where it figures in.

And/or you can do some plots like:


 	repl1 <- aframe$x3==1 # just one replicate
 	matplot( 1:sum(repl1),
 		 predict(afit)[repl1] +
 		  outer(model.matrix(afit)[repl1,'x1b'],-1:1),
 		pch="-o+")

to get a sense of how that coefficient affects the means of the different 
groups.

As to what it 'means', in many discussions one does not mention that 
coefficient nor test it if interactions that involve it are present. Its 
status is much like the intercept in the main effects only model in the 
sense that it is not of much interest by itself, but its absence ( \equiv 
setting its coefficient to zero) usually implies an unwanted and unnatural 
constraint.


Chuck

>
>> n=3
>> aframe=expand.grid(x1=c('a','b'),x2=c('u','v','w'),x3=1:n)
>>
>> set.seed(0)
>> Y=rnorm(dim(aframe)[[1]])
>> aframe=cbind(aframe,Y)
>>
>> afit=aov(Y~x1+x2, aframe)
>> coefficients(afit)
> (Intercept)         x1b         x2v         x2w
> -0.14740700 -0.27974772  1.00234567 -0.01278949
>>
>> afit=aov(Y~x1*x2, aframe)
>> coefficients(afit)
> (Intercept)         x1b         x2v         x2w     x1b:x2v     x1b:x2w
> -0.27108992 -0.03238187  0.61269558  0.74790937  0.77930018 -1.52139771
>
>
>> The numbered curves are the terms of polynomial contrasts.
>>
>> See
>>
>>        ?poly
>>        ?contr.poly
>> and note
>>
>>        all( contr.poly(5) == poly(1:5,degree=4) )
>>
>>
>>
>>
>>>
>>> BTW, what does 'schmextbook' mean?
>>>
>>
>> Well, repeating a word and prepending 'schm' (in place of the initial
>> consonant sound, if any) is a slang idiom for asserting low regard for the
>> idea expressed by the word.
>>
>> As in 'Birthday-Schmirthday' when moaning about getting older:
>>
>>
>>  http://www.berkeleydailyplanet.com/issue/2009-01-29/article/32134?headline=First-Person-Birthday-Schmirthday-Whinings-on-Mortality
>>
>>
>> or  'Middle-schmiddle' which expresses low regard for the 'middle' policy,
>> see
>>
>>        http://aussiethule.blogspot.com/2006/03/olmerts-middle-schmiddle.html
>>
>> Or 'Nilsson-Schmilsson', a self-deprecating title for an album by Harry
>> Nilsson:
>>
>>        http://en.wikipedia.org/wiki/Nilsson_Schmilsson
>>
>> Or 'textbook-schmextbook', in which I express the idea that a textbook
>> really isn't needed.
>>
>> Chuck
>>
>>> ______________________________________________
>>> R-help at 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.
>>>
>>
>> Charles C. Berry                            (858) 534-2098
>>                                            Dept of Family/Preventive
>> Medicine
>> E mailto:cberry at tajo.ucsd.edu               UC San Diego
>> http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego 92093-0901
>>
>>
>
> ______________________________________________
> R-help at 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.
>

Charles C. Berry                            (858) 534-2098
                                             Dept of Family/Preventive Medicine
E mailto:cberry at tajo.ucsd.edu	            UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego 92093-0901