[BioC] Interaction contrasts with RCBD with replicates

[Gordon K Smyth]

Dear Daniel,

You ask, "is there a way to get simultaneously the correct contrasts for 
both main effects and interactions using limma", but I don't understand 
what you are asking.  Correct them for what?  For example, it would make 
no sense to correct interactions for a block effect, when it is the 
interaction with the block that you are interested in.

It seems to me that the first analysis you give below is satisfactory and 
allows you to extract any contrasts you want.  This includes both 
interactions and "main effects".

You have a balanced two way factorial design.  Even if it makes sense to 
consider Block as random from a scientific point of view (and I'm not 
convinced that it does, although of course I don't know the background to 
your experiment), there is nothing to be gained from a statistical point 
of view from fitting Block as random in the linear model.  For your 
balanced design, there is no information about the treatments in the 
between block error stratum.  It's all in the within block strata.  So the 
fixed effect model gives you a full information analysis of the treatment 
effects.

There is no surprise that the random effect model below gives different 
results.  Not only does it treat the block as random, but it also assumes 
an additive rather than an interaction model.  In addition, the within 
block correlation estimate is pooled over genes, so it is not equivalent 
to a classical univariate mixed model.

In summary, it still seems to me that the obvious analysis (your first 
analysis below) is the right one, and I am not understanding what else you 
are trying to achieve.

Best wishes
Gordon

PS. BTW, I don't personally think there is much merit in the classical 
statistical definition of "main effects" for a factorial experiment in a 
genomic context.  Who wants treatment effects averaged over an important 
covariate when an interaction exists?  It would usually be more 
informative to know the separarate effects for each parity.  I would 
prefer to test for interaction; if interaction exists, then give separate 
results for the two Parity levels; if not, fit the additive model and 
report the resulting treatment effects, which are adjusted for baseline 
block differences.

---------------------------------------------
Professor Gordon K Smyth,
Bioinformatics Division,
Walter and Eliza Hall Institute of Medical Research,
1G Royal Parade, Parkville, Vic 3052, Australia.
http://www.wehi.edu.au
http://www.statsci.org/smyth

On Tue, 3 Jul 2012, dantayrod at stat.ufl.edu wrote:

> Dear Gordon,
>
> First of all, thanks for your response.  Yes, this model was utterly
> incorrect, what I aimed to do was to specify the correct denominator for
> the treatment contrasts.
>
> Using the factorial form for the interaction will get me the right
> contrasts for the interaction terms only, but I am concerned since the
> contrasts change for the main effects of the treatment when using the
> TRT*Block form and when using just TRT in the formula and specifying the
> Block as a random effect (they are also interested in the main effects of
> TRT).
>
> In my code TRT (levels: SFA, EFA, Control) denotes the treatment and
> Parity (Levels: Primip, Multip) is the block, so when I run the model and
> contrasts using (in my original code I have more contrasts to test for the
> interactions, not included here):
>
> #Set up design matrix
> TS <- paste(pheno.Liver$Parity, pheno.Liver$TRT, sep=".")
> TS <- factor(TS, levels=c("Primip.Control",
> "Primip.SFA","Primip.EFA","Multip.Control","Multip.SFA", "Multip.EFA"))
> design <- model.matrix(~0+TS)
> colnames(design) <- c(levels(TS))
>
> #Fit model and contrasts
> fit <- lmFit(eset,design)
> MatContrast=makeContrasts(FAT=(Primip.SFA+Primip.EFA+Multip.SFA+Multip.EFA)/4-(Primip.Control+Multip.Control)/2,FA=(Primip.EFA+Multip.EFA)/2-(Primip.SFA+Multip.SFA)/2,levels=design)
> fitMat<-contrasts.fit(fit,MatContrast)
> Contrast.eBayes=eBayes(fitMat)
> TopContrast1=topTable(Contrast.eBayes,coef=1)
> TopContrast2=topTable(Contrast.eBayes,coef=2)
>
> I get different contrast results for the main effects than when running:
>
> #Set up design matrix
> TS2 <- pheno.Liver$TRT
> TS2 <- factor(TS, levels=c("Control", "SFA","EFA"))
> blk <- pheno.Liver$Parity
> design2 <- model.matrix(~0+TS2)
> colnames(design) <- c(levels(TS2))
> correl=duplicateCorrelation(eset, design2,block=blk)
>
> #Fit model and contrasts
> fit2 <- lmFit(eset,design2,block=blk,cor=correl$consensus)
> MatContrast2=makeContrasts(FAT=(SFA+EFA)/2-Control,FA=EFA-SFA,levels=design2)
> fitMat2<-contrasts.fit(fit2,MatContrast2)
> Contrast.eBayes2=eBayes(fitMat2)
> TopContrast2.1=topTable(Contrast.eBayes2,coef=1)
> TopContrast2.2=topTable(Contrast.eBayes2,coef=2)
>
> Is there any way to get simultaneously the correct contrasts for both main
> effects and interactions using limma?
>
> Thanks again for all of your help,
>
> Daniel Taylor
>
>

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