Gordon Smyth Thank you for developing such an amazing tool.

I have noticed that when I generate a model with covariates I see differences in p-values between a treatment-contrast parameterization vs group-means parameterization. The fold-changes are the same. I do not notice this with models that do not have covariates.

My code is the following for treatment-contrasts:

design<- model.matrix(~Status + Age + Sex + PMI, metadata)
fitA <- lmFit(batchregressed, design)
fitB <- eBayes(fitA, trend = FALSE, robust = TRUE)
tabA<- topTable(fitB, number = 2500, coef = 2)
coefA<- fitB$coefficients

My code is the following for group-means:

design2<- model.matrix(~0+Status + Age + Sex + PMI, metadata)
contrasts<- makeContrasts(StatusControl-StatusMS, levels=colnames(design2))
fitC <- lmFit(batchregressed, design2)
fitD <- contrasts.fit(fitC, contrasts)
fitE <- eBayes(fitD, trend = FALSE, robust = TRUE)
tabC<- topTable(fitE, number = 2500)
coefC<- fitE$coefficients

Is there anything that I am misspecifying? Would it be predicted that p-values would differ with these two model parametrization types?

Thanks in advance for your assistance.

If your data includes precision weights in the batchregressed object, then the p-values will differ slightly, if it doesn't then they will be identical. For an explanation of when p-values will differ, type ?contrasts.fit and see the section called Note.