[R] lme and ordering of terms

[Spencer Graves]

	  Since I have not seen a reply to this question, I will offer 
something:  The problem with testing interactions before main effects is 
that it's not clear what interactions even mean without the main 
effects.  This is true in virtually any context, not just lme.

	  Consider the following simulated data with two factors at two levels:

 > set.seed(1)
 > DF <- data.frame(x1=rep(LETTERS[1:2], 3),
+      x2=rep(letters[3:4], each=3), y=rnorm(6))

	  If you read the code for "lm", you will find that R translates the 
explanatory variables this into numbers using "model.matrix" something 
like the following:

 > options(contrasts=c("contr.treatment", "contr.poly"))# R default
 > model.matrix(~x1*x2, DF)
   (Intercept) x1B x2d x1B:x2d
1           1   0   0       0
2           1   1   0       0
3           1   0   0       0
4           1   1   1       1
5           1   0   1       0
6           1   1   1       1
attr(,"assign")
[1] 0 1 2 3
attr(,"contrasts")
attr(,"contrasts")$x1
[1] "contr.treatment"

attr(,"contrasts")$x2
[1] "contr.treatment"

	  S-Plus uses contr.helmert rather than contr.treatment by default, so 
it's "model.matrix" looks different":

 > options(contrasts=c("contr.helmert", "contr.poly"))# S-Plus default
 > model.matrix(~x1*x2, DF, contrasts=contr.helmert)
   (Intercept) x11 x21 x11:x21
1           1  -1  -1       1
2           1   1  -1      -1
3           1  -1  -1       1
4           1   1   1       1
5           1  -1   1      -1
6           1   1   1       1
attr(,"assign")
[1] 0 1 2 3
attr(,"contrasts")
attr(,"contrasts")$x1
[1] "contr.helmert"

attr(,"contrasts")$x2
[1] "contr.helmert"

	  CONCLUSION:  The definition of "interaction" depends on the choice of 
contrast coding.

	  Fortunately, however, if we consider main effects before 
interactions, we get the same analysis of variance regardless of 
contrast coding:

 > options(contrasts=c("contr.treatment", "contr.poly"))# R default
 > fit.t <- lm(y~x1*x2, DF)
 > anova(fit.t)
Analysis of Variance Table

Response: y
           Df  Sum Sq Mean Sq F value Pr(>F)
x1         1 0.72873 0.72873  0.4958 0.5543
x2         1 0.53283 0.53283  0.3625 0.6083
x1:x2      1 0.24469 0.24469  0.1665 0.7228
Residuals  2 2.93980 1.46990

 > options(contrasts=c("contr.helmert", "contr.poly"))# S-Plus default
 > fit.h <- lm(y~x1*x2, DF)
 > anova(fit.h)
Analysis of Variance Table

Response: y
           Df  Sum Sq Mean Sq F value Pr(>F)
x1         1 0.72873 0.72873  0.4958 0.5543
x2         1 0.53283 0.53283  0.3625 0.6083
x1:x2      1 0.24469 0.24469  0.1665 0.7228
Residuals  2 2.93980 1.46990
 >
	  In brief, the analysis of variance is unchanged, because regression 
is equivalent to projecting "y" considered as a 6 dimensional vector 
onto different subspaces.  We first project y onto the subspace spanned 
by the column of all 1's and x1.  Then we add x2 to that subspace and 
project y onto that larger dimensional subspace.  Then we add the x1:x2 
interaction.  Changing the contrast coding changes the choice of 
coordinate axes but does NOT change the subspaces involved UNLESS you 
extract the terms in a different order.  For more detail, consult any 
good book on regression.

	  spencer graves

Christoph Scherber wrote:

> Dear R users,
> 
> When fitting a lme() object (from the nlme library), is it possible to 
> test interactions *before* main effects? As I understand, R 
> conventionally re-orders all terms such that highest-order interactions 
> come last - but I´d like to know if it´s possible (and sensible) to 
> change this ordering of terms.
> 
> I´ve tried the terms() command (from aov) but I don´t know if something 
> similar exists for lme() objects.
> 
> Thanks a lot for your help!
> 
> Best wishes
> Christoph
> 
> ______________________________________________
> 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