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

[Scott Robinson]

Dear Gordon,

I suppose part of my confusion came from the intercept term and that the zero didn't mean assuming the fit passed through the origin, similar to this post I just came across: https://stat.ethz.ch/pipermail/bioconductor/2005-April/008580.html

Thanks very much for clearing things up for me.

All the best,

Scott


-----Original Message-----
From: Gordon K Smyth [mailto:smyth at wehi.EDU.AU] 
Sent: 05 March 2013 23:52
To: Scott Robinson
Cc: Bioconductor mailing list
Subject: Limma: difficulties making multiple comparisons within one multiple-group model

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:6}}