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]]

