[BioC] unexpected voom! / limma contrast inconsistency

[Gordon K Smyth]

Dear Stephen,

Just for clarification, the class of orthogonal design matrices includes 
matrices of the form model.matrix(~0+a).  contrasts.fit() gives exact 
results for these design matrices even with observation weights.  So 
contrasts.fit() do give exact results in conjunction with voom() for some 
of the most common design matrices.

If you do have a non-orthogonal design matrix then the solution right now 
is to use 0~a or, as Ryan Thompson points out, to re-parametrize the 
design matrix so that the contrast becomes a coefficient.

contrasts.fit() was written as it was for speed, ten years ago when 
computers were much slower than they are now.  I have been asked several 
times to provide an always-exact contrast computation, and voom is the 
prompt for me to finally do this.

Best wishes
Gordon

> Date: Fri, 20 Sep 2013 14:16:26 -0400
> From: Stephen Hoang <stephen.a.hoang at gmail.com>
> To: bioconductor at r-project.org
> Subject: [BioC] unexpected voom! / limma contrast inconsistency
>
> 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()?


______________________________________________________________________
The information in this email is confidential and intend...{{dropped:4}}