[BioC] Using limma for omnibus F tests

[Gordon K Smyth]

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