[R] fit.contrast and interaction terms

[Gregory Warnes]

Hello Berta,

gmodels::fit.contrasts() simply performs a single-variable contrast  
using the model you have specified.   To perform the more involved  
contrasts that you are describing, there are two approaches:

1) use the estimable() function in the gmodels package.   
gmodels::estimable() allows you to provides an arbitrary function to  
be applied to the estimated model parameters, allowing you to perform  
any appropriate (or inappropriate) contrast calculation.

2) Use the contrast functions provided by Frank Harrell's Hmisc  
package.   Frank's functions allow you to specify the desired  value  
or level of each parameter for the contrast, and handle unspecified  
parameters.

I hope this helps,

-Greg



On Oct 9, 2007, at 6:10AM , Berta wrote:

> Dear R-users,
> I want to fit a linear model with Y as response variable and X a  
> categorical variable (with 4 categories), with the aim of comparing  
> the basal category of X (category=1) with category 4.  
> Unfortunately, there is another categorical variable with 2  
> categories which interact with x and I have to include it, so my  
> model is s  "reg3: Y=x*x3". Using fit.contrast to make the contrast  
> (category 1 vs category 4) with options(contrasts=c 
> ("contr.treatment", "contr.poly")), it makes the contrast but just  
> for the basal category of x3, (coincident with the corresponding  
> result of summary(reg3)), so that it is not what I am looking for,  
> and it seems that when I write (contrasts=c("contr.sum",  
> "contr.poly")) before fit.contrast, it adjust for x3. Here I send a  
> SMALL EXAMPLE that tries to imiate my problem.
>
> library(gmodels)
> set.seed(100)
> options(contrasts=c("contr.treatment", "contr.poly"))
> y <- rnorm(100)
> x <-  cut(rnorm(100, mean=y, sd=0.25),c(-4,-1.5,0,1.5,4)) # 4  
> CATEGORIES
> x3 <-  cut(rnorm(100, mean=y, sd=8),c(-50,0,50)) # 2 CATEGORIES
> reg3<-lm(y~ x * x3 )
> summary(reg3) # category 1 vs category 4 estimate: 3.84776 , for  
> basal category of x3
> fit.contrast(reg3,x,c(-1,0,0,1)) # g4 vs g1  # category 1 vs  
> category 4 estimate: 3.84776 ,
> options(contrasts=c("contr.sum", "contr.poly"))
> fit.contrast(reg3,x,c(-1,0,0,1)) # g4 vs g1 #  category 1 vs  
> category 4 estimate:  4.010439
>
> I deduce that the global comparison of category 1 vs category 4 is  
> 4.01, and not 3.84.
>
> Unfortunately, the real variable that interact with x is not  
> categorical but continuous in my real case, and the new model is  
> Y=x*x2. Again, with the contr.treatment option, the comparison of  
> category 1 vs category 4 is done for the intercept, that is, for  
> the numerical variable assumed to be 0. As i am interested in  
> comparing category 1 vs category 4 adjusting for x2 in general  
> terms, I use contr.sum as before, but it does not seem to produce  
> any modification:
>
> x2 <- rnorm(100,mean=y,sd=0.5) # NUMERIC
> reg2 <- lm(y ~ x * x2 )
> summary(reg2) # category 1 vs category 4 estimate: 3.067346 for x2=0
> options(contrasts=c("contr.treatment", "contr.poly"))
> fit.contrast(reg2,x,c(-1,0,0,1)) # g4 vs g1 estimate: 3.067346 for  
> x2=0
> options(contrasts=c("contr.sum", "contr.poly"))
> fit.contrast(reg2,x,c(-1,0,0,1)) # g4 vs g1 estimate: 3.067346 for  
> x2=0
>
> The question is: How could I contrast simulaneously in global terms  
> (not intercept and slope separately) if there are differences among  
> category 1 vs category 4 but adjusting for a continuous variable?  
> Or it is not possible to do so, and I have to contrast both  
> (difference of intercepts and difference of slopes separately) and  
> obtain a global conclussion?
>
> Thanks a lot in advance,
>
> Berta
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> 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.