Hello, I have run into a new issue (for me) in fitting a limma model.

The model I am fitting includes a number of covariates/additional factors, such as ages and weights of individuals associated with a given "array", several random factors, such as race/ethnicity, region, and the factor defining my primary groups of interest.

E.g. model.matrix(~ 0 + groups + race + region + scale(age) + scale(weight), data = full_dataset)

I am making contrasts from pairwise interactions of levels of the groups factor, e.g. - high - medium - medium - low - high - low

mat_data <- # defined earlier
design <- # defined earlier
contrasts <- # defined earlier

When fitting a limma model, is there any difference between the following, and if so, what exactly:

fit <- lmFit(mat_data, design)
fit <- contrasts.fit(fit, contrasts)
fit <- eBayes(fit)

and

fit <- lmFit(mat_data, design)
fit <- eBayes(fit)
fit <- contrasts.fit(fit, contrasts)

Clearly there is some difference, since I see large differences in topTable for each coefficient, particularly in the t-statistic and P-values. However, logFC and AveExpr stay the same, so I believe that the coefficients being estimated are the same.

The only reason I stumbled up this scenario is that I wanted to check the residuals of my fit and found that after calling contrasts.fit I could no longer use residuals or fitted with the MArrayLM object. It seems that after adjusting the model with contrasts.fit, since the coefficients no longer match the design, you can't compute with the original coefficients. So, rather than computing the fit twice, I wanted to see if I could use the fit after the eBayes call to extract the residuals, and then call contrasts.fit to get the contrasts of interest.

Thank you for any help you can provide!

I recommend saving the result of contrasts.fit as a separate object, so that you still have access to the original lmFit result as well. They both have their uses, as you have noted. Run eBayes on both of them. My typical workflow looks something like this:

fit <- lmFit(mat_data, design)
fit <- eBayes(fit)
cfit <- contrasts.fit(fit, contrasts)
cfit <- eBayes(cfit)
topTable(cfit)

I think you already know which order is correct. The order with contrasts.fit before eBayes is standard and correct. You've observed that the second order gives different results so you know it can't be right.

It's easy to see what's going on:

1. lmFit computes coefficients, residual variances and standard errors.

2. contrasts.fit converts the coefficients and standard errors to reflect the contrasts rather than the original design matrix, but does not compute t-statistics or p-values.

3. eBayes computes t-statistics and p-values from the coefficients and standard errors.

So you can see that, if you run contrasts.fit after eBayes without running eBayes again, then the object contains coefficients and standard errors matching the contrasts but t-statistics and p-values matching the original coefficients. Obviously not sensible.

If you run the code as Ryan suggested:

fit <- lmFit(mat_data, design)
fit2 <- contrasts.fit(fit, contrasts)
fit2 <- eBayes(fit2)

then you preserve both the original coefficients and the contrasts. It's no coincidence that the example code is written exactly like this on the help page for contrasts.fit and every time the contrasts.fit function is used in the limma User's Guide.

It is true that you can't compute fitted values or residuals from fit2, but I find that I hardly ever want to do that. I could have programmed constrasts.fit so that it would store both the original coefficients and the contrasts but I didn't want to make the output object any bigger than it needed to be.