[BioC] Interaction contrasts with RCBD with replicates

[dantayrod at stat.ufl.edu]

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