[BioC] unexpected voom! / limma contrast inconsistency

[Ryan C. Thompson]

Hi,

Gordon recently answered a request of mine to split out limma's 
design-matrix-reparametrizing logic into a separate function, which he 
added under the name "contrastAsCoef":

https://stat.ethz.ch/pipermail/bioc-devel/2013-September/004652.html

This function takes a design matrix and a set of N contrasts and 
returns an equivalent design matrix where those contrasts are now the 
first N coefficients in the design matrix.

Maybe you could use this to reparametrize your design matrix prior to 
lmFit so that you never have to use contrasts.fit.

-Ryan

On Fri 20 Sep 2013 11:16:26 AM PDT, Stephen Hoang wrote:
> Hi All,
>
> We have come across somewhat unexpected behavior from voom! / limma that is
> worth sharing. We found that using the same data, but different design
> matrices, "equivalent" contrasts can yield different results.
>
> We constructed a design matrix, such that to make our desired contrast
> (treatment - control), we compared two coefficients of the resulting linear
> model using contrasts.fit(). We also constructed a separate design matrix
> in such a way that the control was the baseline of the linear model. In
> principle, the statistics from the contrast.fit() comparison in the first
> design should be equivalent to the statistics for the treatment coefficient
> in the second design. However, we found that this is not the case. Fold
> changes and average expression values for each gene do not differ between
> the two designs; however, the t-statistics, and P-values do.
>
> The reason for this difference is noted in the contrasts.fit() help page in
> the limma vignette, but this caveat is sufficiently subtle that it may fly
> under many radars (as it did with ours). From the help page: "...if the
> design matrix is non-orthogonal and the original fit included quality
> weights or missing values, then the unscaled standard deviations produced
> by this function are approximate rather than exact." Since voom() always
> generates quality weights, the unscaled standard deviations generated by
> contrasts.fit() will always be approximate for non-orthogonal design
> matrices (the vast majority of design matrices). Consistently, the P-values
> are, on balance, worse for contrasts made with contrasts.fit() in our case
> study. As a sanity check, we forced all of the quality weights to 1, which
> indeed caused the two comparisons to generate the same statistics.
>
> Does anyone have any suggestions on how to avoid this problem for
> comparisons that require the use of contrasts.fit()?
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor