Hi Natasha, It turns out it depends on whether you set the contr.sum() on your factor before calling model.matrix(). It's actually in the limmaUsersGuide() in section 8.7 - you have to look closely at the treatment-contrast parametrization (2nd example) and the sum to zero parametrization (3rd example) and the Comparison section of the tables. Here's an example: #1st do treatment-contrast parametrization > Sex <- factor(rep(c("Female","Male"),each=3),levels=c("Male","Female")) > Sex [1] Female Female Female Male Male Male Levels: Male Female > model.matrix(~Sex) (Intercept) SexFemale 1 1 1 2 1 1 3 1 1 4 1 0 5 1 0 6 1 0 attr(,"assign") [1] 0 1 attr(,"contrasts") attr(,"contrasts")$Sex [1] "contr.treatment" #So yes, in this design matrix the 2nd coef is Female - Male #Now switch to the sum to zero parametrization by setting contr.sum() on the Sex factor: > contrasts(Sex) <- contr.sum(2) > Sex [1] Female Female Female Male Male Male attr(,"contrasts") [,1] Male 1 Female -1 Levels: Male Female > model.matrix(~Sex) (Intercept) Sex1 1 1 -1 2 1 -1 3 1 -1 4 1 1 5 1 1 6 1 1 attr(,"assign") [1] 0 1 attr(,"contrasts") attr(,"contrasts")$Sex [,1] Male 1 Female -1 # In this design matrix, the 2nd coef is now (Male - Female)/2. I don't why it works like this, but it messed me up many times before I figured it out! Cheers, Jenny > At 04:26 AM 6/30/2011, Natasha Sahgal wrote: >Hi Jenny and the BioC list, > >Please correct me if I am wrong here, but isn't the contrast "Female >- Male" rather than "Male - Female" as you say? > >I was of the understanding that the contrasts, when not specified by >contrasts.matrix, by default is the latter - former (e.g. B - A), in >a 2 group comparison (say when levels = c(A,B)). > > >BW, >Natasha > >On 29/06/2011 15:31, Jenny Drnevich wrote: >>Actually, you can do ANY sort of factorial design in limma. It took >>me awhile to figure out how to expand the sum to zero >>parametrization (3rd example, Section 8.7 limmaUsersGuide() ) when >>you have more than 2 levels of any factor. Here's an example of a 2x3 design: >> >>#First set up the 2 group factor, just like in the limma vignette: >> >> > Sex <- factor(rep(c("Female","Male"),each=6),levels=c("Male","Female")) >># Note that the order you specify the groups in the levels argument >>determines the direction of the comparison. See below. >> >> > contrasts(Sex) <- contr.sum(2) >> > Sex >> [1] Female Female Female Female Female Female Male Male Male >>Male Male Male >>attr(,"contrasts") >> [,1] >>Male 1 >>Female -1 >>Levels: Male Female >># Note that the contrast is male - female because female was listed >>last in the levels argument above >> >>#Now set up the 3 group factor: >> >> > Time <- factor(rep(1:3,4),levels=3:1) >> # Again, the group order specified in the levels >> determines the direction of comparison >> >> > contrasts(Time) <- contr.sum(3) >> # You use contr.sum(3) here because you have 3 groups. If >> you have 4 groups, you'd use 4, etc. >> >> > Time >> [1] 1 2 3 1 2 3 1 2 3 1 2 3 >>attr(,"contrasts") >> [,1] [,2] >>3 1 0 >>2 0 1 >>1 -1 -1 >>Levels: 3 2 1 >># Note the contrasts are 3 - 1 and 2 - 1, because group 1 was >>listed last in the levels argument above >> >> > design <- model.matrix(~Sex*Time) >> > design >> (Intercept) Sex1 Time1 Time2 Sex1:Time1 Sex1:Time2 >>1 1 -1 -1 -1 1 1 >>2 1 -1 0 1 0 -1 >>3 1 -1 1 0 -1 0 >>4 1 -1 -1 -1 1 1 >>5 1 -1 0 1 0 -1 >>6 1 -1 1 0 -1 0 >>7 1 1 -1 -1 -1 -1 >>8 1 1 0 1 0 1 >>9 1 1 1 0 1 0 >>10 1 1 -1 -1 -1 -1 >>11 1 1 0 1 0 1 >>12 1 1 1 0 1 0 >>attr(,"assign") >>[1] 0 1 2 2 3 3 >>attr(,"contrasts") >>attr(,"contrasts")$Sex >> [,1] >>Male 1 >>Female -1 >> >>attr(,"contrasts")$Time >> [,1] [,2] >>3 1 0 >>2 0 1 >>1 -1 -1 >> >># In this design matrix, the (Intercept) coefficient is the grand >>mean, the Sex1 coef is the main effect of Sex, >># the Time1 and Time2 coef taken together will give you the main >>effect of Time, and the Sex1:Time1 and >># Sex1:Time2 coef taken together will give you the interaction term. >> >># Now fit the model with your data >> >> > fit.2x3 <- eBayes(lmFit(YourData, design)) >> >># To get the overall F test for the 2x3 take all coef except the Intercept: >> >> > overall.2x3 <- topTable(fit, coef=2:6, number=Inf) >> >># To get the F test (equivalent to t test) main effect of Sex, do: >> >> > mainSex <- topTable(fit, coef=2, number=Inf) >># Note that the logFC values need to be multiplied by 2 to get the >>actual Male:Female logFC value! >> >># To get the F test for the main effect of Time, do: >> >> > mainTime <- topTable(fit, coef=3:4, number=Inf) >># I usually don't use the actual logFC values listed for each coef, >>so I usually don't care which group is listed last >> >># To get the F test for the interaction term, do: >> >> > Interact <- topTable(fit, coef=5:6, number=Inf) >> >> >>If you have 4 groups in a level, you'll have 3 coef for the main >>effect of the level and 3 coef for the interaction term, and so on >>up the line. You can also do n-way factorial designs in the same >>way. Just make sure to keep track of the coefficients! >> >>Cheers, >>Jenny >> >>At 03:39 AM 6/29/2011, axel.klenk at actelion.com wrote: >>>Hi Dana, >>> >>>ooops, obviously I shouldn't be posting immediately after lunch... >>>re-reading >>>your message after some cups of coffee I realize that my brain has skipped >>>everything after male/female and treatment/control... so yours is clearly >>>NOT >>>a 2x2 design... my apologies... >>> >>>Cheers, >>> >>> - axel >>> >>> >>>Axel Klenk >>>Research Informatician >>>Actelion Pharmaceuticals Ltd / Gewerbestrasse 16 / CH-4123 Allschwil / >>>Switzerland >>> >>> >>> >>> >>>From: >>>axel.klenk at actelion.com >>>To: >>><dana.stanley at="" csiro.au=""> >>>Cc: >>>bioconductor at r-project.org, bioconductor-bounces at r-project.org >>>Date: >>>28.06.2011 15:13 >>>Subject: >>>Re: [BioC] 2 way anova in Bioconductor >>>Sent by: >>>bioconductor-bounces at r-project.org >>> >>> >>> >>>Hi Dana, >>> >>>limma can do that easily. The current version of the user's guide >>>contains one chapter and two case studies on 2x2 factorial designs. >>> >>>Cheers, >>> >>> - axel >>> >>> >>>Axel Klenk >>>Research Informatician >>>Actelion Pharmaceuticals Ltd / Gewerbestrasse 16 / CH-4123 Allschwil / >>>Switzerland >>> >>> >>> >>> >>> >>>From: >>><dana.stanley at="" csiro.au=""> >>>To: >>><bioconductor at="" r-project.org=""> >>>Date: >>>28.06.2011 03:27 >>>Subject: >>>[BioC] 2 way anova in Bioconductor >>>Sent by: >>>bioconductor-bounces at r-project.org >>> >>> >>> >>>Hi Everyone >>> >>>I use custom designed chicken high-density multiplex one channel >>>microarray (NimbleGen 12x135K). My usual pipeline includes mostly oligo, >>>arrayQualityMetrics, preprocessCore, genefilter, GeneSelector and limma. >>>So far I only worked with simple design, 2 conditions; control/treatment >>>style. Now I have a dataset with embryonic development timeline for male >>>and female chicks. I still compare different timepoints using limma but I >>>want to do 2 way anova to identify interaction significant genes, ie genes >>> >>>where gender influences development timeline. I looked at many packages >>>and could not find 2 way anova, though I am sure it is there somewhere. I >>>tested package ABarray that seemed ideal for my me but after days of >>>trying I contacted the author to find out that expression sets (mine is >>>coming from limma) can no longer be used by ABarray as the package has >>>not been updated for a while. Can anyone suggest a package (and an >>>example code if possible) for 2 way anova? I am so temp! >>> ted to go and do it in GeneSpring... >>> >>>Thanks all >>> >>>Dana