[Draft] pseudo-p significance calculation

[JosiahParry]

This PR drafts a function calculate_significance() to provide a consistent way to calculate pseudo-p values from a reference distribution.

It is based on the discussion at https://github.com/pysal/esda/issues/199

[codecov[bot]]

Codecov Report

:x: Patch coverage is 93.42105% with 5 lines in your changes missing coverage. Please review. :white_check_mark: Project coverage is 82.0%. Comparing base (300f8e8) to head (8dc453d). :warning: Report is 25 commits behind head on main.

Files with missing lines	Patch %	Lines
esda/moran.py	75.0%	2 Missing :warning:
esda/significance.py	95.7%	2 Missing :warning:
esda/crand.py	95.5%	1 Missing :warning:

Additional details and impacted files

@@          Coverage Diff          @@
##            main    #281   +/-   ##
=====================================
  Coverage   82.0%   82.0%           
=====================================
  Files         24      25    +1     
  Lines       3489    3538   +49     
=====================================
+ Hits        2861    2902   +41     
- Misses       628     636    +8     

Files with missing lines	Coverage Δ
esda/crand.py	93.7% <95.5%> (-0.7%)	:arrow_down:
esda/moran.py	84.9% <75.0%> (-0.1%)	:arrow_down:
esda/significance.py	95.7% <95.7%> (ø)

... and 1 file with indirect coverage changes

:rocket: New features to boost your workflow:

• :snowflake: Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[ljwolf]

OK, done here on the logic & implementation. Thank you @JosiahParry for getting the ball rolling here 😄 Very much appreciated!

I've re-implemented the percentile-based two-sided test from scratch using scipy.stats.scoreatpercentile. This approach finds the percentile for the test statistic in the reference distribution and counts how many replicates are outside of (p, 1-p). Over simulations, these are always 2*directed. Second, I modified your folding approach to fold around the mean of the replicates, rather than zero (since the expected value of local stats generally isn't zero) and kept it as a folded option for testing.

I don't think we should expose the folded variant to the user in the classes, since the power in each direction is dependent on the symmetry of the distribution. For example, in the illustration below, the smallest replicate, when folded, is not "extreme," but this is accounted for in the percentile-based method.

The percentile will always equal the folded version for symmetric distributions, but the folded version becomes a directed test as skew increases. I think that the (over+under)/all is also the intended estimand of the directed approach, after re-reading the Hope paper referred to in #199

If other maintainers approve these four options (greater, lesser, two-sided, and directed) for the user classes and a folded option for this function only (for replication/testing purposes) I can start propagating this across the user classes.

[ljwolf]

I am still working on this, but I recall now why the implementation of an "alternative" argument was a bit trickier than I expected... because we allow for the user to "discard" the random replicates, rather than store them, we have to push the significance logic all the way down to the conditional randomization numba function. This may have significant performance implications, since we're currently only counting the number of larger random stats in each innermost loop.

It seems clear to me that

1. if the test statistic is as large as k realizations,
2. then there are always 2k simulations outside of the (k/n, (1-(k/n)) interval, plus the test statistic itself.
3. So, the proper p-value for the two-sided test is (2*n_greater+1)/(n_samples + 1),
4. which is off from two times our current p-value by 1/(n_samples+1).

If 1-4 are correct, this means we don't need to change any of the numba code. The correction can be calculated as 2*directed - (1/(n_samples+1)) after the numba calculation. Do I have this right @sjsrey @knaaptime @martinfleis @jGaboardi?

So, if we implement our current test for local stats without flipping (as greater), generate 1-p_sim (as lesser), and implement the above correction for the two-sided test (2*p_sim - (1/(n_samples+1))), none of the numba code needs to change.

Is that OK w/ other maintainers?

[martinfleis]

This is too much stats for me to say anything useful.

[jGaboardi]

This is too much stats for me to say anything useful.

Same for me.

[ljwolf]

One further wrinkle as well: some global Moran tests support directed testing a binary option already. Notably, if two_tailed=False, they pick the test direction based on whether the global I is positive or negative. It's also useful to note: this means we currently pick the smallest one-tailed p-value and, if the test is two-tailed, multiply this by two.

For us to roll-out the testing across all the classes, we need to consider if this option should be deprecated in favor of an explicit "alternative" option? Right now, there's no way to force a direction on these tests.

[sjsrey]

I am still working on this, but I recall now why the implementation of an "alternative" argument was a bit trickier than I expected... because we allow for the user to "discard" the random replicates, rather than store them, we have to push the significance logic all the way down to the conditional randomization numba function. This may have significant performance implications, since we're currently only counting the number of larger random stats in each innermost loop.

It seems clear to me that

1. if the test statistic is as large as k realizations,
2. then there are always 2k simulations outside of the (k/n, (1-(k/n)) interval, plus the test statistic itself.
3. So, the proper p-value for the two-sided test is (2*n_greater+1)/(n_samples + 1),
4. which is off from two times our current p-value by 1/(n_samples+1).

If 1-4 are correct, this means we don't need to change any of the numba code. The correction can be calculated as 2*directed - (1/(n_samples+1)) after the numba calculation. Do I have this right @sjsrey @knaaptime @martinfleis @jGaboardi?

So, if we implement our current test for local stats without flipping (as greater), generate 1-p_sim (as lesser), and implement the above correction for the two-sided test (2*p_sim - (1/(n_samples+1))), none of the numba code needs to change.

Is that OK w/ other maintainers?

I think this is OK.

One thing to check is if:

3. So, the proper p-value for the two-sided test is (2*n_greater+1)/(n_samples + 1), Should be 2(n_greater+1)/(n_samples+1)

[ljwolf]

One thing to check is if:

Sure, that is what I initially thought & what @JosiahParry suggested.

The reason why I'm thinking it's actually 2*p_sim - (1/(n_permutations + 1)) is because using 2*p_sim amounts to counting the test stat twice: 2*p_directed = 2 * (n_outside + 1)/(n+1) = (2*outside + 2)/(n+1). The difference will be vanishingly small as the number of permutations increases, but it's the principle...

Thinking another way, in the percentile-based version of the test, you compute the percentile p for the test statistic, count how many null statistics are outside of (p,1-p) and add one, since the test stat is always at least at its own percentile. This p-value is smaller than p_sim by 1/(n_samples-1), which is what would happen in the percentile test if you counted the test statistic twice.

the simulation code at the end of esda/simulation.py should illustrate?

[weikang9009]

@ljwolf I was looking at the discussions in this PR and the other related issue. The correction for the two-sided test 2*p_sim - (1/(n_samples+1)) looks correct to me.

[JosiahParry]

Thank you for the explanation @ljwolf. I think I'm almost there/onboard!

It's worth calling out explicitly this formula can result in a p-value > 1.0 which should also be handled e.g.

p_sim = 0.65
nsim = 999

(p_corrected = (2*p_sim - (1/(nsim + 1))))
#> [1] 1.299

if (p_corrected > 1) {
  1.0
} else {
  p_corrected
}
#> [1] 1

Additionally, would you mind elaborating why it is - (1/(nsim + 1)) as opposed to + (1/(nsim + 1))? To me, it makes more sense to penalize smaller numbers of simulations rather than larger number of simulations. For example subtracting the second term results in smaller p values for smaller numbers of simulations and larger ones for larger numbers of simulations.

calc_p_sim <- function(p_sim, nsim) {
  (p_corrected = (2*p_sim - (1/(nsim + 1))))

  if (p_corrected > 1) {
    1.0
  } else {
    p_corrected
  }

}

calc_p_sim(0.05, 49)
#> [1] 0.08
calc_p_sim(0.05, 99)
#> [1] 0.09
calc_p_sim(0.05, 999)
#> [1] 0.099

[weikang9009]

It's worth calling out explicitly this formula can result in a p-value > 1.0 which should also be handled e.g.

The formula 2*p_sim - (1/(n_samples+1)) will not result in a p-value larger than 1 since p_sim is equivalent to the directed case which always picks the smaller side, meaning that p_sim will be always no larger than 0.5.

Additionally, would you mind elaborating why it is - (1/(nsim + 1)) as opposed to + (1/(nsim + 1))? To me, it makes more sense to penalize smaller numbers of simulations rather than larger number of simulations. For example subtracting the second term results in smaller p values for smaller numbers of simulations and larger ones for larger numbers of simulations.

This is based on the equation for the two-sided case p_two_sided = (2*n_greater+1)/(n_samples + 1) and the current implementation of p_sim: p_sim = (n_greater+1)/(n_samples + 1)

@ljwolf please correct me if I'm wrong

[ljwolf]

Yep, exactly right!

@JosiahParry, think about it in terms of the percentiles, maybe that will be clearer.

The correction term is to adjust for double counting the "observed" stat- if you just did 2*p_sim, that computes the number of "outside" simulated statistics, plus the observed stat twice (as if we see it at both percentile p and 1-p) However, we actually only see the observed stat once at percentile p---the two tailed boundary in the "other" tail is derived from 1-p.

Like, imagine the number of simulations is 100. Sort them from smallest to largest. If 20 simulations are as big as the test statistic, then the test statistic is at the 80th percentile. At the (100-80)th percentile, there are n*(1-.8) smaller observations.

So, we see 20 more extreme large observations, 1 test statistic, and 20 more extreme smaller observations as if the test statistic were in the other tail.

2p_sim(n+1) would count 20 more extreme large stats, 20 more extreme small stats, and two test statistics-one too many. So, we need to deflate that count by 1, which deflates the p value by 1/(n+1).

[knaaptime]

for some reason this keeps reminding me of Downey's Euro Problem example, probably just because it reminds me "the statistician"'s default test is always two-sided

[JosiahParry]

Going through my stale PRs. Has this been implemented and documented in the most recent releases? I think some long-ish form documentation on the "improved" calculation of pseudo-p values would be quite great.

[ljwolf]

Hi! back again :) this has not been forgotten---it is the highest priority for me when I have development time.

I needed to move to a pure numpy version of the two-sided percentile test in order to push it down into the numba.njit() inner loop of esda.crand().

~~I will push that code up shortly, bandwidth permitting. My intention is then to write an numba.njit() compatible esda.significance._permutation_test() function, and send both esda.crand.parallel_crand() and esda.significance.permutation_test() to that when calculating the permutation test.~~ done.

We now need to do some profiling, but hopefully there is not much gain from doing it inline vs. calling another jitted function on singly-typed input.

[ljwolf]

@weikang9009 notes correctly that we will also need to update the notebooks where local/global statistics are used before merging this.

[ljwolf]

the affected notebook has been fixed and #376 has been addressed!

[ljwolf]

I think this is ready to merge.

[JosiahParry]

Great work, @ljwolf !

[ljwolf]

Thanks! it seems there's a broadcasting issue that's numba-version dependent. I will squash this issue, and then it's ready.

[ljwolf]

OK, all tests are passing except windows, which look like a build issue. Could we merge this? Or, can I get some help to identify the issue with the windows build? This touches numba code, but no code is specific to windows, and all tests on linux/macos pass.

[ljwolf]

Thanks @martinfleis! @sjsrey can you merge?