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+Mult ip.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}}