[R] using the design matrix to correctly configure contrasts
Karl Brand
k.brand at erasmusmc.nl
Tue Jun 1 17:27:09 CEST 2010
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
More information about the R-help
mailing list