Creating contrast matrix for linear regression in R

Question

I cannot understand the criteria behind contrast matrix, should be created for designing linear models in R. I have read limma user guide (P: 101) and there was an example for that. The example is about estrogen(present/absent) and its late (48h) and early (10h) effect on a cancer. I cannot understand in what sense they could take away the effect of time by assigning 1 for both e48 and E48.

Here are the R commands:

treatments <- factor(c(1,1,2,2,3,3,4,4),labels=c("e10","E10","e48","E48"))
contrasts(treatments) <- cbind(Time=c(0,0,1,1),E10=c(0,1,0,0),E48=c(0,0,0,1))
contrasts

> treatments
[1] e10 e10 E10 E10 e48 e48 E48 E48
attr(,"contrasts")
    Time E10 E48
e10    0   0   0
E10    0   1   0
e48    1   0   0
E48    1   0   1
Levels: e10 E10 e48 E48

Answer 1

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 :-)