This is side-stepping the question, but I am not aware that it ever makes
sense to include the "A:B" term in a design matrix without also including
the main effects of A and B (here I include obvious extensions such as A +
A:B + C where B is a coarser factor than C, so here the main effects of B
are essentially included).  Of course, the A+A:B+C example seems to also
fail in the DESeq2 code you quote, but these examples are rare in comp bio.

If you're just fitting a model like
  A + A:B
you're almost certainly doing something wrong from a statistical point of
view.

Kasper


On Tue, Jul 9, 2013 at 7:00 PM, Steve Lianoglou <lianoglou.steve@gene.com>wrote:

> Greetings DESeq2ers,
>
> I've been playing around with different ways to encode a 2x2 factorial
> design I'm working with: 4 cell types, 4 treatments and exploring
> what's cooking in there.
>
> I thought that the nested interaction formula used in the limma user's
> guide (section 8.5.3) would be an easy way to explore some obvious
> questions since differential expression from different treatments
> within each cell type (among other things) shake out quite cleanly
> from the results of running the wald test since these end up as
> columns in the design matrix, eg:
>
> ~ cell + cell:treatment
>
> as long as the control treatment is the first level of the treatment
> factor, we got columns for all cell:treatment effects, eg.
> cellA:treatment1, cellA:treatment2, ..., cellD:treatment4)
>
> Yay!
>
> ... almost ...
>
> The problem is that the "renaming column" stuff in the fitNbinomGLMs
> function seems to assume that we'll have columns for all main effects
> in the design matrix, so the following line triggers an error:
>
> modelMatrixNames[match(convertNames$from, modelMatrixNames)] <-
> convertNames$to
>
> since you will find elements like "treatment2, ..., treatment4" among
> the elements in `covertNames$from`. The call to match then returns NA
> for these, and *boom*.
>
> An easy fix would be to add the following line after the call to
> `renameModelMatrixColumns`:
>
> convertNames <- subset(convertNames, from %in% modelMatrixNames)
>
> So the "if (renameCols)" block (line 1066 in core.R) now becomes:
>
>   if (renameCols) {
>     names(modelMatrixNames) <- modelMatrixNames
>     convertNames <- renameModelMatrixColumns(modelMatrixNames,
>
>  as.data.frame(colData(object)),
>                                              modelFormula)
>     convertNames <- subset(convertNames, from %in% modelMatrixNames)
>     modelMatrixNames[match(convertNames$from, modelMatrixNames)] <-
> convertNames$to
>   }
>
>
>
> --
> Steve Lianoglou
> Computational Biologist
> Bioinformatics and Computational Biology
> Genentech
>
> _______________________________________________
> Bioc-devel@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/bioc-devel
>

	[[alternative HTML version deleted]]

