[BioC] Limma: difficulties making multiple comparisons within one multiple-group model

[Gordon K Smyth]

Dear Scott,

The code below seems to be from Section 8.3 of the User's Guide (rather 
than 8.6).

Yes, the contrast RNA2 for the second model is the same as RNA2-RNA1 in 
the first model.

The first model doesn't set the intercept to anything, because it has no 
intercept.

Yes, the batch effect will be adjusted for if you include it in the model, 
and this is true for the first model or the second model.

Yes, you could use either of these models to compare RNA3 to the average 
of RNA1 and RNA2.

As far as I can see, there is no difficulty or problem at all.  limma is 
doing doing the right thing for you regardless of how you parametrize the 
model.  You just seem to need some assurance that this is so.

Best wishes
Gordon

> Date: Mon,  4 Mar 2013 14:26:10 -0800 (PST)
> From: "Scott Robinson [guest]" <guest at bioconductor.org>
> To: bioconductor at r-project.org, Scott.Robinson at glasgow.ac.uk
> Subject: [BioC] Limma: difficulties making multiple comparisons within
> 	one	multiple-group model
>
> Dear List,
>
> I am looking at differential methylation (from Illumina 450K) with Limma 
> and have a situation in which I have several groups and want to make 
> comparisons between particular groups. It seems that section 8.6 of the 
> manual is relevant to my situation so I will use that as an example:
>
>> f <- factor(targets$Target, levels=c("RNA1","RNA2","RNA3"))
>> design <- model.matrix(~0+f)
>> colnames(design) <- c("RNA1","RNA2","RNA3")
> To make all pair-wise comparisons between the three groups one could proceed
>> fit <- lmFit(eset, design)
>> contrast.matrix <- makeContrasts(RNA2-RNA1, RNA3-RNA2, RNA3-RNA1,
> + levels=design)
>> fit2 <- contrasts.fit(fit, contrast.matrix)
>> fit2 <- eBayes(fit2)
>
> I am unsure of a couple of things here. One thing I am worried about is 
> whether the comparison RNA2-RNA1 here would be equivelant to doing:
>
>> f <- factor(targets$Target, levels=c("RNA1","RNA2","RNA3"))
>> design <- model.matrix(~f)
>> colnames(design) <- c("(intercept)","RNA2","RNA3")
>> fit <- lmFit(eset, design)
>> fit <- eBayes(fit)
>> topTable(fit, coef="RNA2")
>
> i.e. I don't understand what is going on with the intercept here. If in 
> this second instance we are taking the RNA1 samples to calculate the 
> intercept, then how is it calculated in the previous example? Or is it 
> simply set to (0,0) in the previous example? Or am I way off the mark on 
> how the intercept functions in Limma altogether?
>
> My second problem is that I was wondering if my model were '~0+f+batch' 
> and I then followed through using the first chunk of code would the 
> batch effect be taken into account or not?
>
> Also if I wanted to pool RNA1 and RNA2 and compare these against RNA3 
> can I use this same model?
>
> Thanks in advance, and apologies for ignorance about the intercept 
> issue,
>
> Scott
>
>
> -- output of sessionInfo():
>
> R version 2.15.2 (2012-10-26)
> Platform: x86_64-w64-mingw32/x64 (64-bit)
>
> locale:
> [1] LC_COLLATE=English_United Kingdom.1252
> [2] LC_CTYPE=English_United Kingdom.1252
> [3] LC_MONETARY=English_United Kingdom.1252
> [4] LC_NUMERIC=C
> [5] LC_TIME=English_United Kingdom.1252
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
> [1] WriteXLS_2.3.0 limma_3.14.3
>

______________________________________________________________________
The information in this email is confidential and intend...{{dropped:4}}