[R] using the design matrix to correctly configure contrasts

[Karl Brand]

Esteemed R-forum subscribers,

I'm having a tough time configuring contrasts for my 3-way ANOVA. In short:

I don't know how to configure (all) my contrasts correctly in order to 
specify (all) my comparisons of interest.

I succeeded getting my contrasts of interest set up for a simpler 2-way 
ANOVA based on the fairly intuitive logic of the design col.names.

But i'm not able to extend this logic to a more complex three way ANOVA. 
Obviously there *is* a logic to the "0's" and "1's" of R's design 
'treatment contrast' matrix- i see that the "1's" correspond to each of 
the factors that an observation is associated with. But exactly how this 
information can be used to ascribe contrasts of interest remains beyond 
my understanding, especially where interaction is concerned. Googling 
for days didn't yield the help i was looking for.

Can some one point me to a detailed explanation (suitable for R 
dilettantes) how to specify any contrast of interest given a design matrix?

For anyone really needing some diversion, i'd *really* appreciate a 
personal explanation using my design as an example. Warning- its not 
small...

#load limma package to analyse affy gene expression data
library(limma)

targets <- read.delim("targets.txt")

#factors & levels
Time <- factor(targets$Time, levels = c("t1", "t2", "t3", "t4",
                                         "t5", "t6", "t7", "t8",
                                         "t9", "t10", "t11", "t12",
                                         "t13", "t14", "t15", "t16"))
Tissue <- factor(targets$Tissue, levels = c("R", "C"))
Pperiod <- factor(targets$Pperiod, levels = c("E", "L", "S"))

#model excluding the 3-way interaction (written for convenience to
# build on previous 2-way ANOVA
design <- model.matrix(~Time + Tissue + Pperiod + Time:Tissue + 
Time:Pperiod + Tissue:Pperiod)

fit <- lmFit(rma.pp, design)

#how does this design look:

 > colnames(design)
  [1] "(Intercept)"      "Timet2"
  [3] "Timet3"           "Timet4"
  [5] "Timet5"           "Timet6"
  [7] "Timet7"           "Timet8"
  [9] "Timet9"           "Timet10"
[11] "Timet11"          "Timet12"
[13] "Timet13"          "Timet14"
[15] "Timet15"          "Timet16"
[17] "TissueC"          "PperiodL"
[19] "PperiodS"         "Timet2:TissueC"
[21] "Timet3:TissueC"   "Timet4:TissueC"
[23] "Timet5:TissueC"   "Timet6:TissueC"
[25] "Timet7:TissueC"   "Timet8:TissueC"
[27] "Timet9:TissueC"   "Timet10:TissueC"
[29] "Timet11:TissueC"  "Timet12:TissueC"
[31] "Timet13:TissueC"  "Timet14:TissueC"
[33] "Timet15:TissueC"  "Timet16:TissueC"
[35] "Timet2:PperiodL"  "Timet3:PperiodL"
[37] "Timet4:PperiodL"  "Timet5:PperiodL"
[39] "Timet6:PperiodL"  "Timet7:PperiodL"
[41] "Timet8:PperiodL"  "Timet9:PperiodL"
[43] "Timet10:PperiodL" "Timet11:PperiodL"
[45] "Timet12:PperiodL" "Timet13:PperiodL"
[47] "Timet14:PperiodL" "Timet15:PperiodL"
[49] "Timet16:PperiodL" "Timet2:PperiodS"
[51] "Timet3:PperiodS"  "Timet4:PperiodS"
[53] "Timet5:PperiodS"  "Timet6:PperiodS"
[55] "Timet7:PperiodS"  "Timet8:PperiodS"
[57] "Timet9:PperiodS"  "Timet10:PperiodS"
[59] "Timet11:PperiodS" "Timet12:PperiodS"
[61] "Timet13:PperiodS" "Timet14:PperiodS"
[63] "Timet15:PperiodS" "Timet16:PperiodS"
[65] "TissueC:PperiodL" "TissueC:PperiodS"

#the contrasts that i believe i *can* work out correctly:

cont.matrix <- cbind(

t1.S.RvsAll.S.R.Times =
   c( 0,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,
      0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
      0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
      0,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,
      0,  0),

t1.E.RvsAll.E.R.Times =
   c( 0,	-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
      0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
      0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
      0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
      0,  0),

t1.L.RvsAll.L.R.Times =
   c( 0,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,
      0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
      0,  0,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,  1,
      1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
      0,  0)
                       )

fit2 <- contrasts.fit(fit, cont.matrix)
fit3 <- eBayes(fit2)

#sensical out put follows suggesting i'm on the right track

BUT! Most basically, the contrasts i can NOT work out are:

t1.S.CvsAll.Times & t1.L.CvsAll.Times

ie., similar to t1.S.RvsAll.Times & t1.S.RvsAll.Times above but for the 
"C" tissue, not the "S".




-- 
Karl Brand <k.brand at erasmusmc.nl>
Department of Genetics
Erasmus MC
Dr Molewaterplein 50
3015 GE Rotterdam
P +31 (0)10 704 3409 | F +31 (0)10 704 4743 | M +31 (0)642 777 268