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