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

On Wed, Jan 25, 2023 at 4:42 PM Matt Haberland <[email protected]> wrote:

> 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.
>

and many of those cite statsmodels

https://www.statsmodels.org/dev/generated/statsmodels.stats.multitest.multipletests.html
and a few others
https://www.statsmodels.org/dev/stats.html#multiple-tests-and-multiple-comparison-procedures

many months of work for a few lines of code
(I was pretty happy when I figured out the difference between FWER and FDR
control, and step-up and step-down procedures.)

Josef



> 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.
>

I'm happy about using the improved implementation of the distribution
compared to what we have in statsmodels (it was contributed, taken from
another package).
I'm still waiting for multivariate-t cdf in scipy for more general cases of
multiple comparisons.

I'm not that committed to multiple comparison distributions, because the
usual pvalue corrections work better in many cases (according to some of
the literature)
There is also https://pypi.org/project/scikit-posthocs/ which looks good,
so I lowered my priority for this in statsmodels unless it fits into a
larger topic.


>
> 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.
>

ECDF is trivial and nothing that gets citations or required a lot of work
(various confidence intervals for it are still WIP).
I don't think in the survival case there will be much overlap. The focus in
statsmodels is on regression, and parametric survival/lifetime/reliability
models will arrive once we support multi-link models. (last year numfocus
wasn't happy with my application)

Cheers,
Josef


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

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