[BioC] Limma - small bug in decideTests?

[Gordon Smyth]

Dear Misha,

>Date: Sun, 17 Jun 2007 11:17:06 +0100 (BST)
>From: Misha Kapushesky <ostolop at ebi.ac.uk>
>Subject: [BioC] Limma - small bug in decideTests?
>To: bioconductor at stat.math.ethz.ch
>Message-ID: <Pine.LNX.4.64.0706171104290.17169 at pigeon.ebi.ac.uk>
>Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed
>
>Dear List,
>
>I looked through the archives and didn't notice anyone pointing this out,
>so I thought I'd ask:
>
>If decideTests() is called with the method argument set to hierarchical
>(by the way, that word is misspelled in the implementation as
>'heirarchical' - if someone should call it with correct spelling, they get
>an error), then internally classifyTestsP() is called to adjust across
>contrasts (after adjusting across genes). So far so good (exc the
>spelling).

Thanks for pointing out the spelling. Mind you, the software is 
correct as it stands, because the mispelling is consistent through 
the code and documentation. Nevertheless, I'll correct the spelling.

>However, it seems the multiple adjustment method is not propagated to
>classifyTestsP - so if the user requests adjust.method="BH",
>classifyTestsP will perform the default adjustment, "holm", which is more
>conservative. Is that behavior by design? Or is there just a
>"method=adjust.method" missing on the classifyTestsP() call?

Your question is about hierarchical testing. At the moment, the 
top-level (F-test) multiple adjustment method can be specified by the 
user, but the second level (t-test) multiple testing method is always 
"holm". I honestly can't remember whether this is by design or not. I 
may have been put off by the fact that there is no existing theory 
relating to FDR control for hierarchical tests. Benjamini's group has 
started to publish in the past year on hierarchical tests, but 
nothing yet which applies to the limma situation. So the 
method="hierarchical" should be viewed as experimental. The method 
I've implemented, passing a penalized p-value down to the second 
level, seems reasonable, but it is running ahead of theory to support 
it. I realize that it's high time I wrote a review of the existing 
theory to explain the motivation behind the different decideTests() 
methods. Perhaps you have some ideas of your own?

As a result of your query, I've added "method=adjust.method" to the 
classifyTestsP() call made from decideTests(method="heirarchical"). 
This means that the same multiple testing adjustment will be used for 
both levels of testing. This is as user would probably expect, but 
please view it as experimental for the time being.

>Also one could do fewer computations in the code, perhaps, if one did
>something like
>
>    padj = p.adjust(object$F.p.value, method=adjust.method)
>    sel = padj < 0.05
>    pmax=min(padj[!sel], na.rm=TRUE)
>    results=classifyTestsP(object[sel,], p.value=pmax, method=adjust.method)
>
>for the hierarchical and nestedF branches of decideTests() - this avoids
>computing i, n, and a there. Or are these computed there for a reason,
>again?

Your suggestion is not actually equivalent to the existing the code, 
your 'pmax' being (stochastically) larger than my 'pvalue*a'. The 
existing code is any case very fast so that computation speed isn't an issue.

Best wishes
Gordon

>Thanks,
>--Misha K.