[BioC] unexpected voom! / limma contrast inconsistency

[Gordon K Smyth]

Dear Steve,

You are right.  contrast.fit() is not exact for contrasts in the context 
of an additive model with observational weights.

Gordon

On Wed, 25 Sep 2013, Stephen Hoang wrote:

> Dear Gordon,
>
> Thank you very much for the helpful clarification.
>
> In the case that a blocking variable is given to lmFit (in my case, to
> specify duplicate correlation), along with a design of the form
> model.matrix(~0+a),
> will contrasts.fit() produce approximate or exact results? Intuitively, it
> seems that part of the design is captured in the blocking variable,
> possibly meaning that the overall design is not of the orthogonal class
> that ensures exact results.
>
> Thanks,
> Steve
>
>
> On Sat, Sep 21, 2013 at 7:24 PM, Gordon K Smyth <smyth at wehi.edu.au> wrote:
>
>> 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 intended solely for the
>> addressee.
>> You must not disclose, forward, print or use it without the permission of
>> the sender.
>> ______________________________**______________________________**__________
>>
>

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