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

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.

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.

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.

Best,
Albert



On Wed, Jan 25, 2023 at 5:35 PM <[email protected]> wrote:

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

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