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

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

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

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