Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

Sorry for asking such stupid question, but 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
share|improve this question

migrated from stackoverflow.com Jul 13 '13 at 20:45

This question came from our site for professional and enthusiast programmers.

4  
This question appears to be off-topic because it is about a statistical method and thus belongs on Cross Validated. – Roland Jul 10 '13 at 13:11
up vote 4 down vote accepted

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

share|improve this answer
    
Thanks a lot, that was clear and easy to understand. – Rozita Akrami Jul 10 '13 at 14:29
    
Hi again, I know really a little in math, so I have problem with the notations and terms here in making contrast matrix: I don't know what tensor product is? so can anyone explain what these formulae mean? M_1 : M_2 M_1 * M_2 Maybe some explanation with examples related to above case study. Thanks – Rozita Akrami Jul 16 '13 at 14:29
    
A:B is the interaction term in R's formula notation. A*B is a shorthand for A + B + A:B, so both effects plus the interaction. – January Jul 16 '13 at 14:49
    
Thanks, hopefully my last question: What is intercept mean, and why did you ignored it? Thanks – Rozita Akrami Jul 16 '13 at 15:38
    
Look here: stats.stackexchange.com/questions/63357/… – January Jul 16 '13 at 19:48

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.