Re: ENH: stats: add false discovery rate control function

I have some more thoughts on the specifics of how Statsmodels and SciPy
could relate here. First, perhaps SciPy should stick to including only
Benjamini-Hochberg
and Benjamini-Yekutieli. Currently Statsmodels has three additional
methods. The two stage step-up method Benjamini, Krieger and Yekutieli, a
two stage step-up Benjamini-Hochberg method I'm not familiar with, and an
adaptive step-down method Gavrilov, Benjamini, Sarkar, that I'd never heard
of. There was talk in gh-17402 about adding Storey-Tibshirani. As far as I
can tell, this is currently not included in Statsmodels, and I propose
adding it to Statsmodels. I'd be happy to either submit a PR or help review
one that someone else submits. Do you have any thoughts on this Josef?

Also, I noticed that at the moment the documentation
<https://www.statsmodels.org/dev/generated/statsmodels.stats.multitest.multipletests.html#generated-statsmodels-stats-multitest-multipletests--page-root>
for statsmodels.stats.multitest.multipletests is sparse and at times
slightly inconsistent. Before adding the more detailed documentation page
suggested earlier, it seems like it would be a good idea to fill out and
clean the existing docs up. This could be done on a shorter timescale. I'd
be happy to either submit a PR or help review one. I can submit an issue to
Statsmodels with more details of what I mean here.

Ralf said

Agreed. The question is whether it should be that new API that may
> materialize (xref https://github.com/statsmodels/statsmodels/issues/8627),
> or the old one.
>

Josef, for a new API in Statsmodels, what do you think of the idea of
splitting the methods between two functions,
statsmodels.stats.multitest.false_discovery_control and
statsmodels.stats.multitest.fwer_control? I think it could help better
demarcate the differences between these tasks, allow the documentation for
each to only focus on the task at hand, allow SciPy and Statsmodels to
share flags for the "bh" and "by" methods common between them, and overall
just decrease the surface area of where we would try to maintain
API/behavior compatibility. Since we'd hope to have some level of
compatibility, Josef I think your input on SciPy's false_discovery_control
API/behavior would be really appreciated.

To summarize and lay things out clearly, I think following through on these
things

   1. Cleaning up and filling out
   statsmodels.stats.multitest.multipletests' current documentation
   2. Adding a page to Statmodels' documentation which concisely and
   clearly states the problem of multiple comparisons, describes the
   difference between FWER and FDR control, gives an overview of each of the
   methods describing the pros and cons of each and the situations in which
   one would or would not use them.
   3. Referring to Statsmodels in Scipy's documentation for
   false_discovery_control.
   4. Pointing out the importance of adjusting for multiple comparisons
   with info on how to do it in the documentation for each of SciPy's
   vectorized hypothesis tests.
   5. Updating the API for Statsmodels, making a mostly compatible
   statsmodels.stats.multitest.false_discovery_control which differs from
   SciPy mainly in having more methods available.
   6. Adding Storey-Tibshirani to Statsmodels.

would leave both SciPy and Statsmodels in a stronger position and help
spread awareness of the need to correct for multiple comparisons. I think
all of this work should be done in any case, and it's just a matter of
dividing up the benefits fairly.

Best,
Albert





On Fri, Jan 27, 2023 at 2:57 PM Pamphile Roy <[email protected]> wrote:

> Hi everyone,
>
> Thank you all for your inputs. Based on that, I propose we do the
> following:
>
> We add a recommendation for the corresponding method of statsmodels in the
> docstring of the method.
> Besides the accompanying text, there are 2 ways we can reference the
> method(s): using a plain text ref or url, or we make a real external link.
>
> Currently we do not have an inter-sphinx mapping for statsmodels, we could
> add one. This would strengthen our connection with statsmodels
> and open the possibility to add more API links in the documentation. This
> can either be done here, or considered in a follow up to make it stand out
> and communicated widely even.
>
> Cheers,
> Pamphile
>
>
>
>
> On 27.01.2023, at 20:12, Ralf Gommers <[email protected]> wrote:
>
>
>
> On Thu, Jan 26, 2023 at 3:00 PM <[email protected]> wrote:
>
>> > I noticed the docstring for statsmodels.stats.multitest.multipletests
>> mentions that the API could change in the future.
>>
>> It has been stable for almost 10 years.
>> However, I recently decided on a rewrite to modernize it into a new
>> function (for backwards compatibility)
>> https://github.com/statsmodels/statsmodels/issues/8627
>> (blame the authors of the original papers who did not always provide
>> formulas for the pvalues :)
>> To some extent the policy for backwards compatibility and reproducibility
>> across versions is stricter in statsmodels than in scipy.stats.
>>
>> I don't mind (much) some overlap between scipy.stats and statsmodels, as
>> long as it doesn't become excessive.
>>
>
> I feel the same way. A bit of overlap is sometimes unavoidable, and
> probably even a good thing as long as that overlap is compatible in
> API/behavior. But there shouldn't be too much duplication, we should just
> cross-link and help each other out.
>
>
>> Multiple testing correction is a favorite of mine because the ratio
>> effort per line of code is among the highest that I have ever had.
>>
>> to the docs:
>> A long time ago, I still had a blog
>> FWER and FDR:
>> http://jpktd.blogspot.com/2013/04/multiple-testing-p-value-corrections-in.html
>>
>
> I still remember this one and the other post you linked. "10 lines of code
> and it took you 3 weeks?!" is in the same vein as "this issue has been open
> for 2 years and you still haven't fixed it?" that I hadn't learned to
> ignore back then yet:)
>
>
>>
>> Josef
>>
>> On Thu, Jan 26, 2023 at 9:27 AM Albert Steppi <[email protected]>
>> wrote:
>>
>>> Hi Matt and Josef,
>>>
>>> Some thoughts:
>>>
>>> I agree that offering vectorized hypothesis tests presents a danger that
>>> people may start running multiple tests without realizing they need to
>>> correct for them. I think it's important to point out the need to correct
>>> for multiple-comparisons in the docs and also point to utilities that make
>>> it easy to do so. For the latter, we could of course just refer to
>>> statsmodels.stats.multitest.multipletests, but I there are good reasons we
>>> may want to offer some multiple comparison correction methods directly in
>>> SciPy.
>>>
>>> I noticed the docstring for statsmodels.stats.multitest.multipletests
>>> mentions that the API could change in the future. I think it's important
>>> that we can point to something which tries to guarantee a stable API.
>>> There's a lot of scientific code that's just kind of treated as set it and
>>> forget it, scripts passed down through grad students, shared utilities
>>> within particular labs, that sort of thing. Such code is often not DRY, and
>>> it could be frustrating to have to clean things up after an API change,
>>> especially for a grad student or post-doc who is already under a lot of
>>> stress and time pressure. Also, in my opinion, correcting for multiple
>>> comparisons is of fundamental importance. It would be good if we could make
>>> the experience as seamless as possible on the happy path.
>>>
>>> I have some thoughts for a compromise though. I checked through the docs
>>> and found that for scipy.stats.probplot it says
>>>
>>> *    probplot generates a probability plot, which should not be confused
>>> with a Q-Q or a P-P plot. Statsmodels has more extensive functionality of
>>> this type, see statsmodels.api.ProbPlot.*
>>>
>>> Statsmodels has 10 different multiple comparisons methods available.
>>> Perhaps it would be best to only implement a core set in SciPy that we
>>> consider to be the most important, and then make a similar referral to
>>> Statsmodels for more methods. It would be nice if Statsmodels'
>>> multipletests API could be settled down first though. I think it would be
>>> best if the APIs agreed.
>>>
>>
> Agreed. The question is whether it should be that new API that may
> materialize (xref https://github.com/statsmodels/statsmodels/issues/8627),
> or the old one.
>
>
>> Also, there's an issue that there are subtleties involved when deciding
>>> when and how to correct for multiple comparisons, with some disagreement
>>> among statisticians. It would be good to be able to point to an accessible
>>> document which gives a clear overview of the basics, some rules of thumb,
>>> and refers the user to more in depth literature. The Wikipedia article for
>>> the Multiple comparisons problem
>>> <https://en.wikipedia.org/wiki/Multiple_comparisons_problem> isn't very
>>> satisfactory for this purpose in my opinion. If we write such a document,
>>> perhaps we could make it part of Statsmodels' documentation, and point
>>> SciPy users there for more information. This document could go into the
>>> reasons why someone may want to use the methods that would be available in
>>> Statsmodels but not SciPy.
>>>
>>
> This sounds like not only a good compromise on functionality and
> cross-linking, but a good idea. This document may fit better in Statsmodels
> docs than in SciPy's docs, and it'll be useful to users of both packages.
> Perhaps as part of getting the SciPy feature landed, Matt & other
> scipy.stats` regulars can help with this in Statsmodels? That would
> hopefully be a win-win for everyone.
>
>
>
>> I'd hope the referrals to Statsmodels could help make up for any lost
>>> users due to the loss of a monopoly on multiple comparison correction in
>>> Python.
>>>
>>> Statsmodels is important because it offers a friendly API for the kind
>>> of statistical modeling common in clinical trials and the social sciences
>>> which is often done with R or SAS. These kinds of things are painful to do
>>> with only NumPy and SciPy. For something as fundamental as multiple
>>> comparison correction, I think it's good to offer some functionality in
>>> SciPy, but I agree that it's important to think carefully about where the
>>> boundaries should be. The goal shouldn't be to poach Statsmodels' users,
>>> but to strengthen the Python statistical computing ecosystem as a whole. In
>>> my opinion, adding some multiple comparisons methods to SciPy would fall
>>> into the latter category.
>>>
>>
> +1
>
> Cheers,
> Ralf
>
>
>
>>> Best,
>>> Albert
>>>
>> _______________________________________________
> SciPy-Dev mailing list -- [email protected]
> To unsubscribe send an email to [email protected]
> https://mail.python.org/mailman3/lists/scipy-dev.python.org/
> Member address: [email protected]
>
>
> _______________________________________________
> SciPy-Dev mailing list -- [email protected]
> To unsubscribe send an email to [email protected]
> https://mail.python.org/mailman3/lists/scipy-dev.python.org/
> Member address: [email protected]
>

_______________________________________________
SciPy-Dev mailing list -- [email protected]
To unsubscribe send an email to [email protected]
https://mail.python.org/mailman3/lists/scipy-dev.python.org/
Member address: [email protected]