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@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@gmail.com>
>> To: bioconductor@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 inte...{{dropped:10}}

