ENH: stats: add nonparametric one-sample quantile test and CI

[romain-jacob]

Reference issue

Relate to statsmodels/statsmodels#6562
Loosely relate to #10577
(additional statistical methods, including confidence intervals)

What does this implement?

This PR implements a non-parametric approach to compute confidence intervals for quantiles for a given confidence level and input x which is either a set of samples (one-dimensional array_like) or the number of samples available. The confidence intervals are valid if and only if the samples are i.i.d.

Both one-sided and two-sided confidence intervals can be obtained (default is one-sided). The function returns two values: either the bounds for the two one-sided confidence intervals, or the lower and upper bounds of a two-sided confidence interval.
The return values are either the indexes of the bounds (if x is an integer) or sample values (if x is the set of samples). None is returned when there are not enough samples to compute the desired confidence interval.

This approach is well-known, present is statistics textbooks but got somehow overlooked despite its simplicity and usefulness. An implementation in SciPy would help changing that.

Additional information

• The current PR does not include any benchmark for that function, which seemed a bit of an overkill. It can be added it necessary.
• [FIXED] The Sphinx documentation renders fine locally, although I do get some error messages I did not manage to eradicate (Inline interpreted text or phrase reference start-string without end-string)...
• From the code in this PR, it's easy to write a function that returns the minimal number of samples to compute any confidence interval of interest. I have this code lying around but was not sure whether it makes sense to integrate in SciPy as well. Opinions?

[tylerjereddy]

I just added a few drive-by comments--not a substitute for a stats guru review.

[chrisb83]

Thanks for the contribution. At the moment, there are no functions in SciPy that return a CI, #12609 adds a CI to a test. In the same spirit, one could implement a quantile test (see e.g. Conover: Practical nonparam. statistics) and in addition to the pvalue, return the CI. (we already have binom_test to test the probability). It is just a suggestion that came to my mind.

[romain-jacob]

The only CI failure is in sparse.linalg and unrelated to this PR.

Unless there are further comments, I believe this PR is ready to be merged.

[chrisb83]

@romain-jacob have you thought about the following suggestion?

At the moment, there are no functions in SciPy that return a CI, #12609 adds a CI to a test. In the same spirit, one could implement a quantile test (see e.g. Conover: Practical nonparam. statistics) and in addition to the pvalue, return the CI. (we already have binom_test to test the probability). It is just a suggestion that came to my mind.

[romain-jacob]

Hi @chrisb83, yes, sorry I forgot to follow-up on that one.

The method I am proposing here is not related to hypothesis testing, it's an estimation method for quantiles. I do not know the quantile test you are referring to (I did not look it up), but if it is a test of practical relevance, then I guess that could be certainly added too. But I think that should be in a distinct PR, no?

[chrisb83]

The method I am proposing here is not related to hypothesis testing, it's an estimation method for quantiles.

In general, there is a very close relationship between CIs and hypothesis tests. E.g. if you estimate a mean and you can construct confidence intervals, you can define a test and vice versa. I am just wondering about the best way to integrate the functionality into SciPy. At the moment, just returning a CI without an estimate of the quantile somehow looks odd. are there other implementations, e.g. in R, that we can use as a benchmark?

[romain-jacob]

I see your point, and indeed I am aware of the close connection between CI and NHST, but I am not aware of any "test" that explicitly focus on percentiles (except the median). That being said, I know very little about statistics in general, so such tests may well exist, but I could not find any.

I just browsed through the R packages and found QuantileNPCI which seems to be doing about the same thing, although it seems it computes only two-sided CI (not 100% sure yet, I just had a quick look).

Unfortunately I am not proficient in R, so it's not easy for me to quickly compare the two implementations. But I can try computing the CI presented in the QuantileNPCI example; that will already do for a quick check.

[mdhaber]

Yes

My concerns were resolved.

[mdhaber]

Suggested change

	greater than equal likelihood that the population quantile is larger
	greater than or equal likelihood that the population quantile is larger

?

[mdhaber]

Suggested change

	distribution. `T1` is sensitive whether the population quantile is
	distribution. `T1` is sensitive to whether the population quantile is

[mdhaber]

I'm not familiar with this use of "likelihood". Also, what is the likelihood greater than?

[steppi]

Sorry. It’s nonsense. I wrote it on the plane earlier today while dozing. I’ll fix it up.

[steppi]

Naming things is so hard. I'm not happy with the following repr for QuantileTestResult but just wanted to throw something out there..

    >>> rvs = stats.uniform.rvs(size=100, random_state=rng)
    >>> stats.quantile_test(rvs, q=0.6, p=0.75, alternative='greater')
    QuantileTestResult(statistic=64, statistic_sign=1, pvalue=0.00940696592998271)

My thought is that if we're encoding whether the statistic T1 or T2 is used to calculate the p-value with a sign $\pm1$, then ideally the name statistic_sign should make sense for this attribute. If I try to take a step back and forget how wrapped up I am in this though, I don't think it would be immediately clear to me what statistic_sign should mean in this context. Some other names I thought of are statistic_direction and statistic_side.

To me, I think the important thing that should be communicated is that +1 means the test is reporting a p-value quantifying the likelihood of a sample quantile being greater than or equal to the observed value under the null distribution, and a value of -1 means the test is reporting a p-value quantifying the likelihood of a sample quantile being less than or equal to the observed value under the null distribution.

[mdhaber]

Part of the difficulty is that the statistic itself is always a count of values to the same side of the quantile, regardless of the alternative. Could this be refactored such that one is the statistics is changed to its complement? Then statistic_direction could make a lot of sense.

[steppi]

Part of the difficulty is that the statistic itself is always a count of values to the same side of the quantile, regardless of the alternative. Could this be refactored such that one is the statistics is changed to its complement? Then statistic_direction could make a lot of sense.

Part of the difficulty is that the statistic itself is always a count of values to the same side of the quantile, regardless of the alternative. Could this be refactored such that one is the statistics is changed to its complement? Then statistic_direction could make a lot of sense.

I’ll think about it too. Sounds very promising.

[romain-jacob]

Hello there!

Thanks a lot for your help @steppi! I just wanted to say that it's been strange to see the work on this PR unfolding without me raising a finger, but I'm really grateful to see more experienced contributors wrapping things up.

So, many thanks :-)

[steppi]

Hello there!

Thanks a lot for your help @steppi! I just wanted to say that it's been strange to see the work on this PR unfolding without me raising a finger, but I'm really grateful to see more experienced contributors wrapping things up.

So, many thanks :-)

Happy to help @romain-jacob! A group of us stats maintainers got together at SciPy 2023 and started working on wrapping up / reviewing outstanding PRs. Thanks for this valuable contribution!

[steppi]

@mdhaber I decided to go with statistic_type equal to 1 when test statistic T1 is used and 2 when test statistic T2 is used. I started playing with taking the "complement" of one of the statistics so that statistic_direction would become meaningful but it seemed to complicate things.

I think T1 and T2 are nice because for the continuous case each approximates $n$ times CDF evaluated at x=q, the hypothesized value for the $p_{th}$ quantile. Under the null this should just equal $np$. This seems like a simpler way to organize things. Then I think it's intuitive that for the discrete or mixed case you need both statistics since there might be multiple observations with x = q.

[mdhaber]

@steppi isn't this supposed to be statistic, not statistic1? Here and elsewhere.

[steppi]

Yes

[mdhaber]

@romain-jacob Could you check the recent changes for accuracy? When you approve, these corrections have been made, and @steppi takes a final pass through the rendered documentation to look for inconsistencies like the ones I've highlighted, I think it's fine for @steppi to merge.

[romain-jacob]

I just reviewed the current state and double-checked the code and the examples. The change to the statistic_type trick makes sense to me, but I must say I think this makes the function more difficult to understand. I don't think it matters too much though, and I'm confident it is correct.

So all good from my side :-)

[mdhaber]

Looks like most of the stylization issues are fixed. I'm re-reviewing so it's easier to keep track of what remains.

[mdhaber]

Looking pretty good @steppi. I think we should rebase, squashing consecutive commits by the same author. You want to do that, or want me to?

[steppi]

Looking pretty good @steppi. I think we should rebase, squashing all my commits and squashing all your commits. You want to do that, or want me to?

You can do it if you have time. Eating dinner now and won’t have a chance for a while.

[mdhaber]

I've rebased some pretty messy histories before, but this one seems really tough. Would you take a look? If it's too messy, LMK.

[steppi]

I've rebased some pretty messy histories before, but this one seems really tough. Would you take a look? If it's too messy, LMK.

Sure. Will do.

[romain-jacob]

I almost can't believe this finally got merged! Thank you so much @mdhaber for your patience and dedication throughout the process. I've been a messy contributor to work with -__-'

Many thanks also @steppi for helping closing this! I really appreciate it.

[mdhaber]

Congratulations and thanks for seeing it through!