Question: 2 way anova in Bioconductor
0
7.7 years ago by
Dana.Stanley@csiro.au110 wrote:
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 tempted to go and do it in GeneSpring... Thanks all Dana [[alternative HTML version deleted]]
modified 7.7 years ago by Jenny Drnevich1.9k • written 7.7 years ago by Dana.Stanley@csiro.au110
Answer: 2 way anova in Bioconductor
0
7.7 years ago by
Axel Klenk920
Switzerland
Axel Klenk920 wrote:
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 [[alternative HTML version deleted]] _______________________________________________ Bioconductor mailing list Bioconductor at r-project.org https://stat.ethz.ch/mailman/listinfo/bioconductor Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor The information of this email and in any file transmitted with it is strictly confidential and may be legally privileged. It is intended solely for the addressee. If you are not the intended recipient, any copying, distribution or any other use of this email is prohibited and may be unlawful. In such case, you should please notify the sender immediately and destroy this email. The content of this email is not legally binding unless confirmed by letter. Any views expressed in this message are those of the individual sender, except where the message states otherwise and the sender is authorised to state them to be the views of the sender's company. For further information about Actelion please see our website at http://www.actelion.com
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 [[alternative HTML version deleted]] _______________________________________________ Bioconductor mailing list Bioconductor at r-project.org https://stat.ethz.ch/mailman/listinfo/bioconductor Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor The information of this email and in any file transmitted with it is strictly confidential and may be legally privileged. It is intended solely for the addressee. If you are not the intended recipient, any copying, distribution or any other use of this email is prohibited and may be unlawful. In such case, you should please notify the sender immediately and destroy this email. The content of this email is not legally binding unless confirmed by letter. Any views expressed in this message are those of the individual sender, except where the message states otherwise and the sender is authorised to state them to be the views of the sender's company. For further information about Actelion please see our website at http://www.actelion.com _______________________________________________ Bioconductor mailing list Bioconductor at r-project.org https://stat.ethz.ch/mailman/listinfo/bioconductor Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor The information of this email and in any file transmitted with it is strictly confidential and may be legally privileged. It is intended solely for the addressee. If you are not the intended recipient, any copying, distribution or any other use of this email is prohibited and may be unlawful. In such case, you should please notify the sender immediately and destroy this email. The content of this email is not legally binding unless confirmed by letter. Any views expressed in this message are those of the individual sender, except where the message states otherwise and the sender is authorised to state them to be the views of the sender's company. For further information about Actelion please see our website at http://www.actelion.com
Answer: 2 way anova in Bioconductor
0
7.7 years ago by
Maciej Jończyk720 wrote:
Dear Dana, you can try maanova package. I haven't use it extensively but as I know it can handle 2-way anova as well as more complicated designs. Details and examples could be found on developers' page, given at package page: http://www.bioconductor.org/packages/release/bioc/html/maanova.html. HTH, Maciej > Date: Tue, 28 Jun 2011 11:26:40 +1000 > From: <dana.stanley at="" csiro.au=""> > To: <bioconductor at="" r-project.org=""> > Subject: [BioC] 2 way anova in Bioconductor > Message-ID: > <605D22ED4C355C44AE8703F8787832BFB11F40D276 at EXVIC- MBX02.nexus.csiro.au> > > > 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
Answer: 2 way anova in Bioconductor
0
7.7 years ago by
Jenny Drnevich1.9k
United States
Jenny Drnevich1.9k 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 > > > [[alternative HTML version deleted]] > >_______________________________________________ >Bioconductor mailing list >Bioconductor at r-project.org >https://stat.ethz.ch/mailman/listinfo/bioconductor >Search the archives: >http://news.gmane.org/gmane.science.biology.informatics.conductor > > > > >The information of this email and in any file transmitted with it is >strictly confidential and may be legally privileged. >It is intended solely for the addressee. If you are not the intended >recipient, any copying, distribution or any other use of this email is >prohibited and may be unlawful. In such case, you should please notify the >sender immediately and destroy this email. >The content of this email is not legally binding unless confirmed by >letter. >Any views expressed in this message are those of the individual sender, >except where the message states otherwise and the sender is authorised to >state them to be the views of the sender's company. For further >information about Actelion please see our website at >http://www.actelion.com > >_______________________________________________ >Bioconductor mailing list >Bioconductor at r-project.org >https://stat.ethz.ch/mailman/listinfo/bioconductor >Search the archives: >http://news.gmane.org/gmane.science.biology.informatics.conductor > > > > >The information of this email and in any file transmitted with it is >strictly confidential and may be legally privileged. >It is intended solely for the addressee. If you are not the intended >recipient, any copying, distribution or any other use of this email >is prohibited and may be unlawful. In such case, you should please >notify the sender immediately and destroy this email. >The content of this email is not legally binding unless confirmed by letter. >Any views expressed in this message are those of the individual >sender, except where the message states otherwise and the sender is >authorised to state them to be the views of the sender's company. >For further information about Actelion please see our website at >http://www.actelion.com > >_______________________________________________ >Bioconductor mailing list >Bioconductor at r-project.org >https://stat.ethz.ch/mailman/listinfo/bioconductor >Search the archives: >http://news.gmane.org/gmane.science.biology.informatics.conductor
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
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