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

Hi Josef,

Thanks for your contributions to SciPy and statsmodels, and I appreciate
you weighing in on this issue!

The motivation for adding multiple comparison corrections was skimming all
the March 2022 open-access articles in a few top biomedical journals. The
Bonferroni method of FWER control and the Benjamini, Hochberg, and
Yekutieli methods of FDR control appeared frequently. After researching the
topic, I realized that it is somewhat dangerous (scientifically speaking)
for all of SciPy's hypothesis tests to be vectorized for performing
multiple tests simultaneously without also providing a way to limit false
positive rates. So the intent was not to cherry-pick features from other
libraries, but to add these basic features where there was a need.

I wanted to acknowledge that other libraries have related features and,
since all of them include SciPy as a dependency, I suggested that bringing
some of that into SciPy might be helpful. This has worked well in the past;
e.g. statsmodels and pingouin have used SciPy's implementation of the
studentized range distribution since we also added Tukey's HSD test. I had
hoped that this case would be similar, but I understand that there are
differences.

We will still refer users to statsmodels for the more advanced
functionality it offers, but for the ndarray hypothesis tests we have in
SciPy, we need some basic ndarray multiple comparison correction
functionality in SciPy.

I thought I should mention now that a PR for some survival analysis
essentials is on the way. The impetus for the PR is the same review of
biomedical articles mentioned above, but it will also address our
longstanding roadmap item "Add tools for survival analysis" and gh-17431 (
https://github.com/scipy/scipy/issues/17431), which asked for a related
feature in SciPy so that they would not need to add another dependency to
their code.

On Mon, Dec 12, 2022 at 2:39 PM <[email protected]> wrote:

> ```
> - These procedures are so general and influential (~10k citations per
> paper) that they belong closer to the base of the scientific computing
> stack.
> - The fact that all of these statistical packages listed above implement
> the procedure further supports the need to have an implementation toward
> the base of the stack. All the packages listed above have SciPy as a
> dependency, and replacing custom implementations with an upstream
> implementation would reduce their maintenance burden.
> ```
> (partial quote)
>
> While most of this can be considered correct, it leads to scipy
> cherry-picking topics from statsmodels.
> I consider statsmodels as "close to the base of the scientific computing
> stack" as far as it concerns statistics and econometrics.
>
> If scipy starts to cherry-pick our (statsmodels) popular topics, then
> there is not much incentive left for me to spend weeks or months to figure
> out how to do Statistics in Python.
> It's already difficult enough to keep up the motivation after 14 years of
> unpaid work.
>
>
> http://jpktd.blogspot.com/2013/04/multiple-testing-p-value-corrections-in.html
> http://jpktd.blogspot.com/2012/03/10-lines-of-code-and-it-took-you-3.html
>
> statsmodels needs "hot" features to attract users, and not just boring old
> fashioned statistics.
> For example a huge fraction of blog articles are OLS and Logit. But that's
> the entry point for users to move to other features.
>
> Josef
>
>
> On Thu, Dec 8, 2022 at 4:03 PM Matt Haberland <[email protected]>
> wrote:
>
>> Hi team,
>>
>> When many single-inference tests are performed at a given significance
>> level, the rate of Type I errors (incorrect rejection of the null
>> hypothesis) exceeds the significance level. One common approach to dealing
>> with this problem is to adjust p-values to control the family-wise error
>> rate (e.g. Bonferroni correction), but an alternative that tends to be more
>> powerful is to control the false discovery rate: the expected proportion of
>> incorrectly rejected hypotheses.
>>
>> gh-17402 proposes a `false_discovery_control` function that accepts an
>> array of p-values and adjusts them so that if rejection of null hypotheses
>> is based on the adjusted p-values, the false discovery rate is controlled
>> at the desired level. There are similar procedures in other statistics
>> libraries, but given that the papers of the underlying methods are so
>> general and influential (~10k citations per paper) and many of
>> SciPy's hypothesis tests are vectorized to act along slices of
>> N-dimensional arrays, it seems important to have such a function closer to
>> the base of the scientific Python stack.
>>
>> If you're interested, please join the discussion at
>> https://github.com/scipy/scipy/pull/17402
>>
>> Thanks!
>> Matt
>>
>> --
>> Matt Haberland
>> Assistant Professor
>> BioResource and Agricultural Engineering
>> 08A-3K, Cal Poly
>> _______________________________________________
>> 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]
>


-- 
Matt Haberland
Assistant Professor
BioResource and Agricultural Engineering
08A-3K, Cal Poly

_______________________________________________
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]