Question: Re: F-tests for factorial effects - limma
0
14.5 years ago by
Ramon Diaz1.1k
Ramon Diaz1.1k wrote:
Dear Naomi and Gordon, Just in case it helps someone else, I am attaching a more verbose version of the code snippet that Gordon sent, because I found the "diag(p) [,attr(X,"assign")==3]" non-intuitive, so I called makeContrasts explicitly. This is the example: ***************************************************************** ## treat is a factor with three levels; age is a continuous variable. X <- model.matrix( ~ treat*age.centered) fit <- lmFit(d.clean, X) p <- ncol(X) cont.ia <- diag(p)[,attr(d.trt.BY.age3,"assign")==3] fit.ia <- eBayes(contrasts.fit(fit, cont.ia)) ############# Using my more verbose apporach X2 <- model.matrix( ~ treat*age.centered) colnames(X2) <- c("Intercept", "Colon", "Mama", "age", "Colon.by.age", "Mama.by.age") contrasts.trt.BY.age <- makeContrasts(Colon.by.age, Mama.by.age, levels = X2) fit2 <- lmFit(d.clean, X2) fit.ia2 <- eBayes(contrasts.fit(fit2, contrasts.trt.BY.age)) ### We can of course do instead (which I like better) X2 <- model.matrix(~ -1 + treat + age.centered + treat*age.centered) colnames(X2) <- c("Colon", "Mama", "Normal", "age", "Colon.by.age", "Mama.by.age") contrasts.trt.BY.age <- makeContrasts(Colon.by.age, Mama.by.age, levels = X2) ************************************************* Best, R. On Wednesday 22 December 2004 14:23, Gordon K Smyth wrote: > > Date: Tue, 21 Dec 2004 17:21:19 -0500 > > From: Naomi Altman <naomi@stat.psu.edu> > > Subject: [BioC] F-tests for factorial effects - limma > > To: bioconductor@stat.math.ethz.ch > > > > I am analyzing a 2-factor factorial Affy experiment, with 3 d.f. for each > > factor. > > > > I would like to get the F-tests for the main effects and interactions > > using limma. > > > > I have computed all the contrasts, and got the t-tests (both unadjusted > > and eBayes). I do know how to combine these into F-tests "by hand" but I > > could not figure out if there was a simple way to do this using limma. > > limma doesn't have any easy way to deal with main effects and interactions, > at least not with main effects, interactions are actually simpler. I > haven't implemented this because I've never been able to figure out what > one would do with these things in a microarray context. > > To compute F-tests for main effects and interaction, the easiest way would > probably be to compute the SS for main effects and interactions by > non-limma means, then use shrinkVar() to adjust the residual mean squares, > i.e., the F-statistic denominators. > > If you only want F-tests for interactions, the following code would work: > > X <- model.matrix(~a*b) > fit <- lmFit(eset, X) > p <- ncol(X) > cont.ia <- diag(p)[,attr(X,"assign")==3] > fit.ia <- eBayes(contrasts.fit(fit, cont.ia)) > > Now fit.ia contains the F-statistic and p-values for the interaction in > fit.ia$F and fit.ia$F.p.value. > > > I had a look at FStat (classifyTestsF). There seems to be a problem, in > > that the matrix tstat is not premultiplied by the contrast matrix when > > the F-statistics are computed. So, if the contrasts are not full- rank, > > an error is generated (instead of the F-statistics) because nrow(Q) != > > ncol(tstat).. (See the lines below). > > No, the code is correct. FStat is quite happy with non full rank contrasts > but the contrast matrix must be applied using contrasts.fit() before > entering FStat(). You should not expect to see a contrast matrix inside > the classifyTestsF() code. > > Gordon > > > if (fstat.only) { > > fstat <- drop((tstat%*% Q)^2 %*% array(1, c(r, 1))) > > attr(fstat, "df1") <- r > > attr(fstat, "df2") <- df > > return(fstat) > > } > > > > I figured that before I fiddled with the code, I would check to make sure > > that I didn't miss an existing routine to do this. > > > > Thanks in advance. > > > > Naomi S. Altman 814-865-3791 (voice) > > Associate Professor > > Bioinformatics Consulting Center > > Dept. of Statistics 814-863-7114 (fax) > > Penn State University 814-865-1348 (Statistics) > > University Park, PA 16802-2111 > > _______________________________________________ > Bioconductor mailing list > Bioconductor@stat.math.ethz.ch > https://stat.ethz.ch/mailman/listinfo/bioconductor -- Ram?n D?az-Uriarte Bioinformatics Unit Centro Nacional de Investigaciones Oncol?gicas (CNIO) (Spanish National Cancer Center) Melchor Fern?ndez Almagro, 3 28029 Madrid (Spain) Fax: +-34-91-224-6972 Phone: +-34-91-224-6900 http://ligarto.org/rdiaz PGP KeyID: 0xE89B3462 (http://ligarto.org/rdiaz/0xE89B3462.asc) Este correo electronico y, en su caso, cualquier fichero anexo al mismo, contiene informacion exclusivamente dirigida a su destinatario o destinatarios. Si Vd. ha recibido este mensaje por error, se ruega notificar esta circunstancia al remitente. Las ideas y opiniones manifestadas en este mensaje corresponden unicamente a su autor y no representan necesariamente a las del Centro Nacional de Investigaciones Oncologicas (CNIO). The information contained in this message is intended for the addressee only. If you have received this message in error or there are any problems please notify the originator. Please note that the Spanish National Cancer Centre (CNIO), does not accept liability for any statements or opinions made which are clearly the sender's own and not expressly made on behalf of the Centre.
microarray cancer affy limma • 607 views
modified 14.5 years ago by Naomi Altman6.0k • written 14.5 years ago by Ramon Diaz1.1k
Answer: Re: F-tests for factorial effects - limma
0
14.5 years ago by
Naomi Altman6.0k
Naomi Altman6.0k wrote:
Thanks, Ramon. What we did was more straight forward: We applied aov to each row, and extracted all of the mean squares. We used the prior error variance and d.f. as computed by eBayes from the original model. We then computed an eBayes F-test by computing the posterior variance. --Naomi At 02:03 PM 1/5/2005 +0100, Ramon Diaz-Uriarte wrote: >Dear Naomi and Gordon, > >Just in case it helps someone else, I am attaching a more verbose version of >the code snippet that Gordon sent, because I found the "diag(p) >[,attr(X,"assign")==3]" non-intuitive, so I called makeContrasts explicitly. > > >This is the example: >***************************************************************** >## treat is a factor with three levels; age is a continuous variable. > >X <- model.matrix( ~ treat*age.centered) >fit <- lmFit(d.clean, X) >p <- ncol(X) >cont.ia <- diag(p)[,attr(d.trt.BY.age3,"assign")==3] >fit.ia <- eBayes(contrasts.fit(fit, cont.ia)) > > >############# Using my more verbose apporach > >X2 <- model.matrix( ~ treat*age.centered) >colnames(X2) <- c("Intercept", "Colon", "Mama", > "age", "Colon.by.age", "Mama.by.age") >contrasts.trt.BY.age <- makeContrasts(Colon.by.age, > Mama.by.age, > levels = X2) >fit2 <- lmFit(d.clean, X2) >fit.ia2 <- eBayes(contrasts.fit(fit2, > contrasts.trt.BY.age)) > >### We can of course do instead (which I like better) >X2 <- model.matrix(~ -1 + treat + age.centered + > treat*age.centered) >colnames(X2) <- c("Colon", "Mama", "Normal", > "age", "Colon.by.age", "Mama.by.age") >contrasts.trt.BY.age <- makeContrasts(Colon.by.age, > Mama.by.age, > levels = X2) > >************************************************* > >Best, > > >R. > > > > > > >On Wednesday 22 December 2004 14:23, Gordon K Smyth wrote: > > > Date: Tue, 21 Dec 2004 17:21:19 -0500 > > > From: Naomi Altman <naomi@stat.psu.edu> > > > Subject: [BioC] F-tests for factorial effects - limma > > > To: bioconductor@stat.math.ethz.ch > > > > > > I am analyzing a 2-factor factorial Affy experiment, with 3 d.f. for each > > > factor. > > > > > > I would like to get the F-tests for the main effects and interactions > > > using limma. > > > > > > I have computed all the contrasts, and got the t-tests (both unadjusted > > > and eBayes). I do know how to combine these into F-tests "by hand" but I > > > could not figure out if there was a simple way to do this using limma. > > > > limma doesn't have any easy way to deal with main effects and interactions, > > at least not with main effects, interactions are actually simpler. I > > haven't implemented this because I've never been able to figure out what > > one would do with these things in a microarray context. > > > > To compute F-tests for main effects and interaction, the easiest way would > > probably be to compute the SS for main effects and interactions by > > non-limma means, then use shrinkVar() to adjust the residual mean squares, > > i.e., the F-statistic denominators. > > > > If you only want F-tests for interactions, the following code would work: > > > > X <- model.matrix(~a*b) > > fit <- lmFit(eset, X) > > p <- ncol(X) > > cont.ia <- diag(p)[,attr(X,"assign")==3] > > fit.ia <- eBayes(contrasts.fit(fit, cont.ia)) > > > > Now fit.ia contains the F-statistic and p-values for the interaction in > > fit.ia$F and fit.ia$F.p.value. > > > > > I had a look at FStat (classifyTestsF). There seems to be a problem, in > > > that the matrix tstat is not premultiplied by the contrast matrix when > > > the F-statistics are computed. So, if the contrasts are not full-rank, > > > an error is generated (instead of the F-statistics) because nrow(Q) != > > > ncol(tstat).. (See the lines below). > > > > No, the code is correct. FStat is quite happy with non full rank contrasts > > but the contrast matrix must be applied using contrasts.fit() before > > entering FStat(). You should not expect to see a contrast matrix inside > > the classifyTestsF() code. > > > > Gordon > > > > > if (fstat.only) { > > > fstat <- drop((tstat%*% Q)^2 %*% array(1, c(r, 1))) > > > attr(fstat, "df1") <- r > > > attr(fstat, "df2") <- df > > > return(fstat) > > > } > > > > > > I figured that before I fiddled with the code, I would check to make sure > > > that I didn't miss an existing routine to do this. > > > > > > Thanks in advance. > > > > > > Naomi S. Altman 814-865-3791 (voice) > > > Associate Professor > > > Bioinformatics Consulting Center > > > Dept. of Statistics 814-863-7114 (fax) > > > Penn State University 814-865-1348 (Statistics) > > > University Park, PA 16802-2111 > > > > _______________________________________________ > > Bioconductor mailing list > > Bioconductor@stat.math.ethz.ch > > https://stat.ethz.ch/mailman/listinfo/bioconductor > >-- >Ram??n D??az-Uriarte >Bioinformatics Unit >Centro Nacional de Investigaciones Oncol??gicas (CNIO) >(Spanish National Cancer Center) >Melchor Fern??ndez Almagro, 3 >28029 Madrid (Spain) >Fax: +-34-91-224-6972 >Phone: +-34-91-224-6900 > >http://ligarto.org/rdiaz >PGP KeyID: 0xE89B3462 >(http://ligarto.org/rdiaz/0xE89B3462.asc) > > >Este correo electronico y, en su caso, cualquier fichero anexo al mismo, >contiene informacion exclusivamente dirigida a su destinatario o >destinatarios. Si Vd. ha recibido este mensaje por error, se ruega >notificar esta circunstancia al remitente. Las ideas y opiniones >manifestadas en este mensaje corresponden unicamente a su autor y no >representan necesariamente a las del Centro Nacional de Investigaciones >Oncologicas (CNIO). > > >The information contained in this message is intended for the addressee >only. If you have received this message in error or there are any problems >please notify the originator. Please note that the Spanish National Cancer >Centre (CNIO), does not accept liability for any statements or opinions >made which are clearly the sender's own and not expressly made on behalf >of the Centre. Naomi S. Altman 814-865-3791 (voice) Associate Professor Bioinformatics Consulting Center Dept. of Statistics 814-863-7114 (fax) Penn State University 814-865-1348 (Statistics) University Park, PA 16802-2111