DOC: stats: add realistic examples to variance tests

[tupui]

What does this implement/fix?

Adds a realistic example to scipy.stats.levene. The example is from:

C.I. BLISS (1952), The Statistics of Bioassay: With Special
Reference to the Vitamins, pp 499-503,
doi: 10.1016/C2013-0-12584-6.

Additional information

When this looks good, I'll add similar examples to the other variance tests.

[mdhaber]

I'm curious if a variance test was performed in the paper?

I'd be interested in seeing the results of the permutation version of the test. If that portion of the example belongs in the spearmanr/kendalltau, might as well follow the rest of that format here.

Then again, I'm having second thoughts about putting these in all the function documentation. The files are growing considerably, and I think they might be too lengthy (too much plotting code for the amount we use the function itself). Maybe we should move them over into a section of tutorials for hypothesis tests like we do for the distributions?

[tupui]

I'm curious if a variance test was performed in the paper?

It was not done, but I've seen some blog posts (of questionable quality) using this examples for that.

I'd be interested in seeing the results of the permutation version of the test. If that portion of the example belongs in the spearmanr/kendalltau, might as well follow the rest of that format here.

Permutation version? Not sure I understand.

Then again, I'm having second thoughts about putting these in all the function documentation. The files are growing considerably, and I think they might be too lengthy (too much plotting code for the amount we use the function itself). Maybe we should move them over into a section of tutorials for hypothesis tests like we do for the distributions?

I think there are 2 separate things here.

• File length: I think we could think about splitting some of the big files into logical grouping. Similarly to your refactoring of the API doc. It's common for the API doc structure to match the file/code structure. Having huge files like more_stats or stats is not very self explanatory.

• Details/plots: I also have this concern. I think we have 2 ways of proceeding: 1. continue like that and consolidate things at the end. 2. think now about the result we want to have and refactor/document functions. We started 1. which is not optimal, but still an improvement. One benefit is that for a user it's self contained. At the moment users don't really navigate a lot in our docs (from the metrics). But it's a biased assessment since our doc is constructed to favor that.

I propose we discuss that next week.

[mdhaber]

Permutation version? Not sure I understand.

It is possible to use permutation_test to approximate the null distribution similarly to how it is done in the correlation tests, and it is useful to do so if the sample size is small. (If that is not desired here, we should be able to answer the question: why is it desirable to show how to do that in the correlation tests but not here? I think should be consistent about the depth of each example, ideally.)

I think there are 2 separate things here.

Yes, you're right, there are really two issues.

1. Yes, I was planning on that. I want to to be careful to preserve git blame and the history of the code, if possible. I'm hoping that if we copy the big files in one commit so that git recognizes that they have the same history, we can then remove things from each copy to make them unique.
2. Yes, I think proceeding as we are and refactoring at the end is OK. That's what I had in mind.

[tupui]

It is possible to use permutation_test to approximate the null distribution similarly to how it is done in the correlation tests, and it is useful to do so if the sample size is small. (If that is not desired here, we should be able to answer the question: why is it desirable to show how to do that in the correlation tests but not here? I think should be consistent about the depth of each example, ideally.)

We would get the following:

def statistic(x, y, z):
    return stats.levene(x, y, z).statistic

ref = stats.permutation_test((small_dose, medium_dose, large_dose), statistic,
                             permutation_type='pairings')

# PermutationTestResult(statistic=0.6457341109631506, pvalue=1.0, null_distribution=array([0.64573411, 0.64573411, 0.64573411, ..., 0.64573411, 0.64573411, 0.64573411]))

Shall I include it then?

[mdhaber]

I think it would be good to include the stuff about the F distribution being an asymptotic approximation and to show how permutation_test can be used to get a more accurate p-value.

The code is not quite right.

1. Based on how the levene statistic works, it looks like we should only consider values of statistic larger than the observed value to be "more extreme". Therefore, we need to use alternative='greater' in the call to permutation_test.
2. Also, the levene statistic doesn't require the sample sizes to be equal, so we immediately know that method='pairings' and method='samples' are inappropriate. method='independent' is appropriate because we are testing the null hypothesis that the samples are independently drawn from distributions with the same variance.

The call should be:

ref = stats.permutation_test((small_dose, medium_dose, large_dose), statistic, 
                             permutation_type='independent', 
                             alternative='greater')

And we would see a pretty substanctial discrepancy between the asymptotic null distribution and the randomized approximation:

[tupui]

I did the update. (I noticed the colour of the hist does not corresponds to the legend for kendalltau etc, I can adjust here if you want.)

For the rest I am honestly getting lost with wording. Could you make suggestions? Thanks!

[mdhaber]

Thanks @tupui! Yes, go ahead and add this example to the other variance tests.

[tupui]

@mdhaber failure is not related. Depending on how close this is to be merged, I could put the fix here which is just a variable name change in a stats test which was introduced yesterday.

[mdhaber]

I could put the fix here which is just a variable name change in a stats test which was introduced yesterday.

Oops, I took care of that in gh-17865.

[mdhaber]

Looks close. If you wouldn't mind digging a bit into the conditions under which the null distributions used in the tests were derived, I'd appreciate it. (Never mind. The computational experiments don't lie.) The statement about the null distributions being an asymptotic approximation is not always true. In any case, it does not appear to be the most important thing to say for bartlett.

[mdhaber]

The name of this test could be changed. In some PRs I removed the prefix and just made it plot so that it would be easier to copy.

[tupui]

Yeah I will do this 👍

[mdhaber]

I'm not sure if this is true here. It might be, but I'm not sure it's the thing that's worth mentioning.

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
rng = np.random.default_rng(1638083107694713882823079058616272161)

ps = []
ss = []
for i in range(10000):
    samples = rng.normal(size=(5, 5))
    s, p = stats.bartlett(*samples)
    ss.append(s)
    ps.append(p)

x = np.linspace(0, np.max(ss))
dist = stats.chi2(df=len(samples)-1)
plt.plot(x, dist.pdf(x))
plt.hist(ss, density=True, bins=50)
plt.show()

The graph matches the chi2 distribution quite well even for just a few observations per sample.

If we just change to the uniform distribution, it's a very different story:

The following might be more important:

Note that the chi-square distribution provides the null distribution when the observations are normally distributed. For samples drawn from non-normal populations, it may be more appropriate to perform a permutation test...

But it would probably be worth doing some research to confirm.

[tupui]

Thanks Matt. Since these are just examples I propose to only write we know to be true. I don't think there is much value to gain, for the purpose of the examples, going deeper.

[mdhaber]

Since these are just examples I propose to only write we know to be true.

But we don't know that this is true. It doesn't appear to be true.

We do know that bartlett is sensitive to non-normality. It is mentioned in the docstring.

So my suggestion was to remove the part we don't know and replace it with something we do know:

Note that the chi-square distribution provides the null distribution when the observations are normally distributed. For small samples drawn from non-normal populations, it may be more appropriate to perform a permutation test...

[tupui]

Thanks Matt!