Dear all, I wanted to adjust for continuous covariates using limma. I tried 2 approaches which so far work for 2 group designs as explained in limma user guide on page 40. However, if I add a single covariate/blocker into the equation, I get slightly different statistics results in spite of the same coefficient values.
First example:
design <- model.matrix(~ group+covar)
design
(Intercept) group1 covar
1 1 0 0.2352934
2 1 0 0.2497340
3 1 1 0.3405839
4 1 0 0.3807273
5 1 1 0.3847266
6 1 1 0.3096032
(Intercept) group1 covar
653635 2.0700705 -0.04815855 2.031578
729737 4.7529897 1.02896914 -6.946158
100133331 2.2457020 0.73567908 -3.811348
100288069 1.8392895 0.19191507 -2.054743
100287934 -0.8755875 0.19874871 3.724989
400728 1.7266959 0.38078881 -1.048634
y <- voom(dgeObj,design)
prefit <- lmFit(y,design)
fit <- eBayes(prefit)
logFC AveExpr t P.Value adj.P.Val B
55007 -1.5740320 5.953331 -6.612608 0.0002165173 0.9999677 -4.014277
Second example:
design1 <- model.matrix(~ 0 + covar)
design1
group0 group1 covar
1 1 0 0.2352934
2 1 0 0.2497340
3 0 1 0.3405839
4 1 0 0.3807273
5 0 1 0.3847266
6 0 1 0.3096032
CvsCON
653635 -0.04815855
729737 1.02896914
100133331 0.73567908
100288069 0.19191507
100287934 0.19874871
400728 0.38078881
contrasts <- makeContrasts(CvsCON = group1 - group0, levels = design1)
y1 <- voom(dgeObj,design1)
prefit1 <- lmFit(y1,design1)
fit.cont <- contrasts.fit(prefit1, contrasts)
fit1 <- eBayes(fit.cont)
logFC AveExpr t P.Value adj.P.Val B
55007 -1.5740320 5.953331 -6.756469 0.0001879744 0.9999673 -3.986743
Can anyone help me to figure out where the difference in topTable stat calculation is? Thank you very much in advance!
1 answer
Already figured it out. Thank you!
If the design matrix is not orthogonal, contrasts.fit() uses an approximation (The approximation only arises when there is an additive term). In the case that design matrix is orthogonal, contrasts.fit() is going to produce exact values as they are with the included intercept (~group).
Log in to answer this question.