On Mon, Aug 04, 2008 at 10:17:38AM +0200, Peter Dalgaard wrote:
> 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.)
That is a neat idea, thanks, Peter, but it doesn't quite fit the bill.
The summary provides t-tests but I am hoping to find F-tests,
otherwise I'm not sure how to efficiently test A (3 levels) at the two
levels of C.
The anova.lme function doesn't help, sadly:
> anova(lme(Y ~ M2 - 1, random = ~ 1 | Block, data = example))
numDF denDF F-value p-value
M2 6 25 23.0198 <.0001
so I'm still flummoxed!
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/
______________________________________________
R-help@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.
