Fix to p-value calculation for permutation and exact

src/bmdscore/trend.cpp

@@ -79,7 +79,7 @@ double cwilcox_new(int k, int m, int n) {
   return (w[i][j][k]);
 }
 
-int cwilcox(int k, int m, int n) {
+double cwilcox(int k, int m, int n) {
   // reserve memory for w
   w_init(m, n);
 
@@ -90,5 +90,5 @@ int cwilcox(int k, int m, int n) {
   // if (m > WILCOX_MAX || n > WILCOX_MAX) w_free(m, n);
   w_free(m, n);
 
-  return ((int)ret);
+  return ret;
 }

src/bmdscore/trend.h

@@ -4,6 +4,6 @@
 const int WILCOX_MAX = 100;
 // int imax(int i, int j);
 // void w_init(int m, int n);
-int cwilcox(int k, int m, int n);
+double cwilcox(int k, int m, int n);
 
 #endif

src/pybmds/stats/jonckheere.py

@@ -45,7 +45,7 @@ def calculate_pvalue(self, decreasing: float, increasing: float) -> float:
 class Approach(StrEnum):
     exact = "exact"
     approximate = "approximate"
-    permtuation = "permutation"
+    permutation = "permutation"
 
 
 class TestResult(BaseModel):
@@ -126,16 +126,16 @@ def jonckheere(
 
     if nperm:
         # calculate P-value via permutation
-        approach = Approach.permtuation
+        approach = Approach.permutation
         rng = np.random.default_rng(seed)
         observed_stat = statistic
         perm_stats = np.empty(nperm, dtype=np.float64)
         x_perm = x.copy()
         for i in range(nperm):
             perm_stats[i] = _calc_stat(x_perm, group_indices, group_count)
             x_perm = rng.permutation(x)
-        decreasing = np.sum(perm_stats >= observed_stat) / nperm
-        increasing = np.sum(perm_stats <= observed_stat) / nperm
+        decreasing = np.sum(perm_stats <= observed_stat) / nperm
+        increasing = np.sum(perm_stats >= observed_stat) / nperm
 
     else:
         # calculate P-value via normal approximation
@@ -151,8 +151,12 @@ def jonckheere(
         # calculate p-value with exact method
         else:
             approach = Approach.exact
-            decreasing = sum(_conv_pdf(group_size)[: int(statistic) + 1])
-            increasing = 1 - sum(_conv_pdf(group_size)[: int(statistic)])
+            pdf = _conv_pdf(group_size)
+
+            stat = int(statistic)
+
+            decreasing = np.sum(pdf[: stat + 1])  # P(J <= observed)
+            increasing = 1 - np.sum(pdf[:stat])  # P(J >= observed)
 
     pval = hypothesis.calculate_pvalue(decreasing, increasing)
 

tests/test_pybmds/stats/test_jonckheere.py

@@ -1,7 +1,7 @@
 import numpy as np
 import pytest
 
-from pybmds.stats.jonckheere import jonckheere
+from pybmds.stats.jonckheere import Approach, jonckheere
 
 
 @pytest.fixture
@@ -49,12 +49,16 @@ def test_hypothesis_paths(self, hypothesis, valid_x, valid_group):
         result = jonckheere(valid_x, valid_group, hypothesis=hypothesis)
         assert 0 <= result.p_value <= 1
 
-        # non-unique x
-        result = jonckheere(
-            np.repeat(np.linspace(1, 10, 10), 2),
-            np.repeat(np.linspace(1, 5, 5), 4),
-            hypothesis=hypothesis,
-        )
+        with pytest.warns(
+            UserWarning,
+            match="P-value estimated using normal distribution; total observations < 30",
+        ):
+            result = jonckheere(
+                np.repeat(np.linspace(1, 10, 10), 2),
+                np.repeat(np.linspace(1, 5, 5), 4),
+                hypothesis=hypothesis,
+            )
+
         assert 0 <= result.p_value <= 1
 
         # permutations
@@ -99,3 +103,163 @@ def test_one_group(self, valid_x):
         group = np.ones(8)
         with pytest.raises(ValueError, match="Only one group has non-missing data"):
             jonckheere(valid_x, group, hypothesis="two-sided")
+
+    def test_exact_increasing_trend_small_pvalue(self):
+        """Exact method should give a small p-value for a strong increasing trend."""
+        x = np.array(
+            [
+                1,
+                2,
+                3,
+                10,
+                11,
+                12,
+                20,
+                21,
+                22,
+            ],
+            dtype=float,
+        )
+
+        group = np.array(
+            [
+                0,
+                0,
+                0,
+                1,
+                1,
+                1,
+                2,
+                2,
+                2,
+            ],
+            dtype=float,
+        )
+
+        result = jonckheere(x, group, hypothesis="increasing")
+
+        assert result.approach == Approach.exact
+        assert result.p_value < 0.05
+
+    def test_exact_decreasing_trend_small_pvalue(self):
+        """Exact method should give a small p-value for a strong decreasing trend."""
+        x = np.array(
+            [
+                20,
+                21,
+                22,
+                10,
+                11,
+                12,
+                1,
+                2,
+                3,
+            ],
+            dtype=float,
+        )
+
+        group = np.array(
+            [
+                0,
+                0,
+                0,
+                1,
+                1,
+                1,
+                2,
+                2,
+                2,
+            ],
+            dtype=float,
+        )
+
+        result = jonckheere(x, group, hypothesis="decreasing")
+
+        assert result.approach == Approach.exact
+        assert result.p_value < 0.05
+
+    def test_permutation_increasing_trend_small_pvalue(self):
+        """Permutation method should use the upper tail for an increasing trend."""
+        x = np.array(
+            [
+                1,
+                2,
+                3,
+                10,
+                11,
+                12,
+                20,
+                21,
+                22,
+            ],
+            dtype=float,
+        )
+
+        group = np.array(
+            [
+                0,
+                0,
+                0,
+                1,
+                1,
+                1,
+                2,
+                2,
+                2,
+            ],
+            dtype=float,
+        )
+
+        result = jonckheere(
+            x,
+            group,
+            hypothesis="increasing",
+            nperm=10000,
+            seed=123,
+        )
+
+        assert result.approach == Approach.permutation
+        assert result.p_value < 0.05
+
+    def test_permutation_decreasing_trend_small_pvalue(self):
+        """Permutation method should use the lower tail for a decreasing trend."""
+        x = np.array(
+            [
+                20,
+                21,
+                22,
+                10,
+                11,
+                12,
+                1,
+                2,
+                3,
+            ],
+            dtype=float,
+        )
+
+        group = np.array(
+            [
+                0,
+                0,
+                0,
+                1,
+                1,
+                1,
+                2,
+                2,
+                2,
+            ],
+            dtype=float,
+        )
+
+        result = jonckheere(
+            x,
+            group,
+            hypothesis="decreasing",
+            nperm=10000,
+            seed=123,
+        )
+
+        assert result.approach == Approach.permutation
+        assert result.p_value < 0.05