R help - Aug 2008 - Decomposing tests of interaction terms in mixed-effects models

Dear R colleagues,

a friend and I are trying to develop a modest workflow for the problem
of decomposing tests of higher-order terms into interpretable sets of
tests of lower order terms with conditioning.

For example, if the interaction between A (3 levels) and C (2 levels)
is significant, it may be of interest to ask whether or not A is
significant at level 1 of C and level 2 of C.

The textbook response seems to be to subset the data and perform the
tests on the two subsets, but using the RSS and DF from the original
fit.  

We're trying to answer the question using new explanatory variables.
This approach (seems to) work just fine in aov, but not lme.  

For example,

##########################################################################

# Build the example dataset

set.seed(101)

Block <- gl(6, 6, 36, lab = paste("B", 1:6, sep = ""))
A <- gl(3, 2, 36, lab = paste("A", 1:3, sep = ""))
C <- gl(2, 1, 36, lab = paste("C", 1:2, sep = ""))
example <- data.frame(Block, A, C) 
example$Y <- rnorm(36)*2 + as.numeric(example$A)*as.numeric(example$C)^2 + 
        3 * rep(rnorm(6), each=6)

with(example, interaction.plot(A, C, Y)) 

# lme 

require(nlme) 
anova(lme(Y~A*C, random = ~1|Block, data = example)) 

summary(aov(Y ~ Error(Block) + A*C, data = example))

# Significant 2-way interaction.  Now we would like to test the effect
# of A at C1 and the effect of A at C2.  Construct two new variables
# that decompose the interaction, so an LRT is possible.

example$AC <- example$AC1 <- example$AC2 <- interaction(example$A,
example$C)

example$AC1[example$C == "C1"] <- "A1.C1"  # A is
constant at C1
example$AC2[example$C == "C2"] <- "A1.C2"  # A is
constant at C2

# lme fails (as does lmer)

lme(Y ~ AC1 + AC, random = ~ 1 | Block, data = example)

# but aov seems just fine.

summary(aov(Y ~ Error(Block) + AC1 + AC, data = example)) 

## AC was not significant, so A is not significant at C1

summary(aov(Y ~ Error(Block) + AC2 + AC, data = example)) 

## AC was significant, so A is significant at C2

## Some experimentation suggests that lme doesn't like the 'partial
## confounding' approach that we are using, rather than the variables
## that we have constructed.

lme(Y ~ AC, random = ~ 1 | Block, data = example)
lme(Y ~ A + AC, random = ~ 1 | Block, data = example)

##########################################################################

Are we doing anything obviously wrong?   Is there another approach to
the goal that we are trying to achieve?

Many thanks,

Andrew

-- 
Andrew Robinson  
Department of Mathematics and Statistics            Tel: +61-3-8344-6410
University of Melbourne, VIC 3010 Australia         Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/

Andrew Robinson wrote:
> Dear R colleagues,
>
> a friend and I are trying to develop a modest workflow for the problem
> of decomposing tests of higher-order terms into interpretable sets of
> tests of lower order terms with conditioning.
>
> For example, if the interaction between A (3 levels) and C (2 levels)
> is significant, it may be of interest to ask whether or not A is
> significant at level 1 of C and level 2 of C.
>
> The textbook response seems to be to subset the data and perform the
> tests on the two subsets, but using the RSS and DF from the original
> fit.  
>
> We're trying to answer the question using new explanatory variables.
> This approach (seems to) work just fine in aov, but not lme.  
>
> For example,
>
> ##########################################################################
>
> # Build the example dataset
>
> set.seed(101)
>
> Block <- gl(6, 6, 36, lab = paste("B", 1:6, sep =
""))
> A <- gl(3, 2, 36, lab = paste("A", 1:3, sep = ""))
> C <- gl(2, 1, 36, lab = paste("C", 1:2, sep = ""))
> example <- data.frame(Block, A, C) 
> example$Y <- rnorm(36)*2 + as.numeric(example$A)*as.numeric(example$C)^2
+
>         3 * rep(rnorm(6), each=6)
>
> with(example, interaction.plot(A, C, Y)) 
>
> # lme 
>
> require(nlme) 
> anova(lme(Y~A*C, random = ~1|Block, data = example)) 
>
> summary(aov(Y ~ Error(Block) + A*C, data = example))
>
> # Significant 2-way interaction.  Now we would like to test the effect
> # of A at C1 and the effect of A at C2.  Construct two new variables
> # that decompose the interaction, so an LRT is possible.
>
> example$AC <- example$AC1 <- example$AC2 <- interaction(example$A,
example$C)
>
> example$AC1[example$C == "C1"] <- "A1.C1"  # A is
constant at C1
> example$AC2[example$C == "C2"] <- "A1.C2"  # A is
constant at C2
>
> # lme fails (as does lmer)
>
> lme(Y ~ AC1 + AC, random = ~ 1 | Block, data = example)
>
> # but aov seems just fine.
>
> summary(aov(Y ~ Error(Block) + AC1 + AC, data = example)) 
>
> ## AC was not significant, so A is not significant at C1
>
> summary(aov(Y ~ Error(Block) + AC2 + AC, data = example)) 
>
> ## AC was significant, so A is significant at C2
>
> ## Some experimentation suggests that lme doesn't like the 'partial
> ## confounding' approach that we are using, rather than the variables
> ## that we have constructed.
>
> lme(Y ~ AC, random = ~ 1 | Block, data = example)
> lme(Y ~ A + AC, random = ~ 1 | Block, data = example)
>
> ##########################################################################
>
> Are we doing anything obviously wrong?   Is there another approach to
> the goal that we are trying to achieve?
>   
This looks like a generic issue with lme/lmer, in that they are not 
happy with rank deficiencies in the design matrix.

Here's a "dirty" trick which you are not really supposed to know
about
because it is hidden inside the "stats" namespace:

M <- model.matrix(~AC1+AC, data=example)
M2 <- stats:::Thin.col(M)
summary(lme(Y ~ M2 - 1, random = ~ 1 | Block, data = example)

(Thin.col(), Thin.row(), and Rank() are support functions for 
anova.mlm() et al., but maybe we should document them and put them out 
in the open.)

-- 
   O__  ---- Peter Dalgaard             ?ster Farimagsgade 5, Entr.B
  c/ /'_ --- Dept. of Biostatistics     PO Box 2099, 1014 Cph. K
 (*) \(*) -- University of Copenhagen   Denmark      Ph:  (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)              FAX: (+45) 35327907