Using limma for omnibus F tests
0
0
Entering edit mode
@gordon-smyth
Last seen 32 minutes ago
WEHI, Melbourne, Australia
Dear Gavin, Thanks for the complete code example. I think you're perceiving a problem which isn't really there, because your limma code already does exactly what you want to do. It finds genes which are different between a, b and c. The reason your F-statistics for genes 26:50 aren't quite the same as for 1:25 is that your data isn't the same. Ordinary univariate F-statistics wouldn't repeat either. You can do omnibus F-tests even more easily by > mdl <- model.matrix(~treatment, pdat) > fit <- eBayes(lmFit(dat,mdl)) > topTable(fit[,-1]) Best wishes Gordon > Date: Wed, 3 Sep 2008 13:35:48 +0100 > From: "Gavin Kelly" <gavin.kelly at="" cancer.org.uk=""> > Subject: [BioC] Using limma for omnibus F tests > To: <bioconductor at="" stat.math.ethz.ch=""> > > Dear all, > I'm trying to find genes that fail the null hypothesis that they are all > the same across >2 treatment groups using the limma package (which I > find exceptionally useful, by the way). > > If I were doing this for a single gene, I'd use a traditional F test; > limma's eBayes generates an F value (via classifyTestsF) from the > individual t statistics, but if I use a model without an intercept, I > believe the F reported will favour genes that are zero in all treatment > groups. And if I use an intercept I'm struggling to make any contrasts > or manipulations of the design matrix that respect the required symmetry > between all treatment groups > > I'd like a way to more generally report an F that reflects genes that > are the same (not necessarily zero) in all groups; I suspect there's a > clever way to do this using contrasts.fit (I thought maybe looking at > all pairwise comparisons which doesn't work in the snippet belo) but > finding that is predicated on my being clever! Anyone got any > suggestions? > > library(limma) > set.seed(1) > pdat <- expand.grid(treatment=c("a","b","c"), rep=c("r1","r2","r3")) > sig <- matrix(0,nrow=50, ncol=nrow(pdat)) > sig[1:25, pdat$treatment=="a"] <- 3 # 25 are up in a > sig[26:50, pdat$treatment=="b"] <- 3 # next 25 are up in b > err <- matrix(rnorm(25*nrow(pdat)), nrow=25, ncol=nrow(pdat)) > dat <- sig + rbind(err, err) #add the same noise to both > dat <- rbind(dat, err, err, err, err) # add some non differentials > > mdl <- model.matrix(~treatment-1, pdat) > ctr <- makeContrasts(treatmentb-treatmenta, treatmentc-treatmenta, > treatmentc-treatmentb,levels=mdl) > fit <- eBayes(contrasts.fit(lmFit(dat, mdl), ctr)) > plot(fit$F) # I'd like this to repeat 1:25 == 26:50 > > Regards - Gavin > -- > Gavin Kelly > Bioinformatics & Biostatistics > Cancer Research UK
Cancer limma Cancer limma • 923 views
ADD COMMENT

Login before adding your answer.

Traffic: 635 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