Interaction contrasts with RCBD with replicates
1
0
Entering edit mode
Guest User ★ 13k
@guest-user-4897
Last seen 9.6 years ago
Dear all, I have a RCBD with three treatment levels and two blocks. Within each block, I have 5 replicates for each treatment (that is, my block size is 15, 5 for each trt level). I am aware that since I have replicates within the blocks, the denominator for the treatment is supposed to be the treatment by block interaction. However, the person who I am analyzing the data for is interested in comparing the interaction (block*trt) and I am not quite certain how to setup the model statement in lmFit in this case, so that I can eventually test the interaction. Does the code below fit the model and tests using the correct residuals (i.e. TRT*Block for TRT and Block, and sqrt(MSE) for TRT*Block)? design <- model.matrix(~ 0+TRT*Block, pheno.Data) correl=duplicateCorrelation(eset, design,block=Block) fit <- lmFit(eset,design,block=Block,cor = correl$consensus) I sincerely appreciate any help I can get. Thanks, Daniel -- output of sessionInfo(): > sessionInfo() R version 2.15.0 (2012-03-30) Platform: x86_64-pc-mingw32/x64 (64-bit) locale: [1] LC_COLLATE=English_United States.1252 [2] LC_CTYPE=English_United States.1252 [3] LC_MONETARY=English_United States.1252 [4] LC_NUMERIC=C [5] LC_TIME=English_United States.1252 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] gcrma_2.28.0 BiocInstaller_1.4.6 car_2.0-12 [4] nnet_7.3-1 MASS_7.3-18 affy_1.34.0 [7] Biobase_2.16.0 BiocGenerics_0.2.0 limma_3.12.1 loaded via a namespace (and not attached): [1] affyio_1.24.0 Biostrings_2.24.1 IRanges_1.14.3 [4] preprocessCore_1.18.0 splines_2.15.0 stats4_2.15.0 [7] zlibbioc_1.2.0 > -- Sent via the guest posting facility at bioconductor.org.
• 1.8k views
ADD COMMENT
0
Entering edit mode
@gordon-smyth
Last seen 4 hours ago
WEHI, Melbourne, Australia
Hi Daniel, You're making it more complicated than it needs to be. You are setting Block to be a fixed factor (in the design matrix) and to be a random factor (using the block argument) at the same time, but this is not permitted. Just use model.matrix(~TRT*Block), forget about "0+" in the formula and forget about duplicate correlation, and then it's all straightforward and correct: design <- model.matrix(~TRT*Block, pheno.Data) fit <- lmFit(eset,design) fit <- eBayes(fit) topTable(fit, coef=5:6) Best wishes Gordon On Mon, 2 Jul 2012, Daniel [guest] wrote: > > Dear all, > > I have a RCBD with three treatment levels and two blocks. Within each > block, I have 5 replicates for each treatment (that is, my block size is > 15, 5 for each trt level). > > I am aware that since I have replicates within the blocks, the > denominator for the treatment is supposed to be the treatment by block > interaction. > > However, the person who I am analyzing the data for is interested in > comparing the interaction (block*trt) and I am not quite certain how to > setup the model statement in lmFit in this case, so that I can > eventually test the interaction. > > Does the code below fit the model and tests using the correct residuals > (i.e. TRT*Block for TRT and Block, and sqrt(MSE) for TRT*Block)? > > design <- model.matrix(~ 0+TRT*Block, pheno.Data) > correl=duplicateCorrelation(eset, design,block=Block) > fit <- lmFit(eset,design,block=Block,cor = correl$consensus) > > I sincerely appreciate any help I can get. > > Thanks, > > Daniel > > -- output of sessionInfo(): > >> sessionInfo() > R version 2.15.0 (2012-03-30) > Platform: x86_64-pc-mingw32/x64 (64-bit) > > locale: > [1] LC_COLLATE=English_United States.1252 > [2] LC_CTYPE=English_United States.1252 > [3] LC_MONETARY=English_United States.1252 > [4] LC_NUMERIC=C > [5] LC_TIME=English_United States.1252 > > attached base packages: > [1] stats graphics grDevices utils datasets methods base > > other attached packages: > [1] gcrma_2.28.0 BiocInstaller_1.4.6 car_2.0-12 > [4] nnet_7.3-1 MASS_7.3-18 affy_1.34.0 > [7] Biobase_2.16.0 BiocGenerics_0.2.0 limma_3.12.1 > > loaded via a namespace (and not attached): > [1] affyio_1.24.0 Biostrings_2.24.1 IRanges_1.14.3 > [4] preprocessCore_1.18.0 splines_2.15.0 stats4_2.15.0 > [7] zlibbioc_1.2.0 >> > > -- > Sent via the guest posting facility at bioconductor.org. > ______________________________________________________________________ The information in this email is confidential and intend...{{dropped:4}}
ADD COMMENT
0
Entering edit mode
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
ADD REPLY
0
Entering edit mode
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}}
ADD REPLY

Login before adding your answer.

Traffic: 805 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6