This example from limma is a bit confusing, because what they do is: they first create a contrast matrix, then they build a design matrix, and then they build another contrast matrix that they actually use. Also, I'm a bit confused, because the contrasts that they end up using do not seem the correct contrasts to me, but this may be due to my limited understanding.
What matters for you are only the two latter, design
and cont.matrix
:
design <- model.matrix(~treatments)
colnames(design) <- c("Intercept","Time","E10","E48")
The design
matrix describes the dummy variables that are used to code the linear regression model:
> design
Intercept Time E10 E48
1 1 0 0 0
2 1 0 0 0
3 1 0 1 0
4 1 0 1 0
5 1 1 0 0
6 1 1 0 0
7 1 1 0 1
8 1 1 0 1
Variable one is the Intercept; all samples are multiplied by 1
, so an average of all samples is the intercept. Second variable is time, and only the late timepoints will go into this. E10 is the mean of the samples with estrogen present at the early timepoint, E48 likewise at the late timepoint.
Then, the authors create another contrast matrix:
> cont.matrix <- cbind(E10=c(0,0,1,0),E48=c(0,0,0,1))
> cont.matrix
E10 E48
[1,] 0 0
[2,] 0 0
[3,] 1 0
[4,] 0 1
Exactly the same contrasts will be achieved with makeContrasts( "E10", "E48", levels= design )
. I have no idea what the purpose of this contrast fitting is, since the two contrasts should be exactly equivalent to the two coefficients in design
. And guess what? They are, check it for yourself:
fit <- lmFit( eset, design )
fit2 <- contrasts.fit( fit, cont.matrix )
fit <- eBayes( fit ) ; fit2 <- eBayes( fit )
cor( fit$coefficients[,"E10"], fit2$coefficients[,"E10"] )
# result: 1; you can also check the p-values, they are identical
Now, to see what is going on I would use an equivalent, but conceptually much simpler encoding.
d <- model.matrix( ~ 0 + treatments ) # no intercept!
colnames( d ) <- levels( treatments )
Now the d looks like this:
e10 E10 e48 E48
1 1 0 0 0
2 1 0 0 0
3 0 1 0 0
4 0 1 0 0
5 0 0 1 0
6 0 0 1 0
7 0 0 0 1
8 0 0 0 1
We encode each group separately; each coefficient represents simply the average of log-signals from microarrays corresponding to the samples from one group only. Likewise, the contrasts are easier and simple to understand:
cmtx <- makeContrasts( "E10-e10", "E48-e48", levels= d )
Here, we see that we are directly testing the difference between E10 and e10 in one contrast, and E48 and e48 in the other. cmtx
looks like that:
> cmtx
Contrasts
Levels E10-e10 E48-e48
e10 -1 0
E10 1 0
e48 0 -1
E48 0 1
If you calculate these contrasts with
fit3 <- eBayes( contrasts.fit( lmFit( eset, d ), cmtx ) )
you will see that the results are exactly the same as before.. but now the setup makes much more sense :-)