Repository: statsmodels/statsmodels

Issue:
BUG: tost uses min pvalue instead of max pvalue
(mostly copied from https://github.com/statsmodels/statsmodels/pull/8166#issuecomment-1086704606)

serious bug in tost_poisson_2indep

it looks like it uses min pvalue and not max pvalue
`idx_max = 0 if tt1.pvalue < tt2.pvalue else 1`

the unit test checking against the confidence interval didn't catch this, because AFAIU it has both pvalue equal, so min and max are the same.

```
def test_tost_poisson():

    count1, n1, count2, n2 = 60, 51477.5, 30, 54308.7
    # # central conf_int from R exactci
    low, upp = 1.339735721772650, 3.388365573616252

    res = smr.tost_poisson_2indep(count1, n1, count2, n2, low, upp,
                                  method="exact-cond")

    assert_allclose(res.pvalue, 0.025, rtol=1e-12)
```

discovered while reading and adjusting the code for nonequivalence test https://github.com/statsmodels/statsmodels/issues/3418
see also https://github.com/statsmodels/statsmodels/issues/8212
We need to make sure that consistency checks check all relevant cases

tost_prob_superior might have the same problem, copying the same pattern of test function and unit test

example:

t0 is the unit test in test_tost_poisson

```
t0 = tost_poisson_2indep(count1, n1, count2, n2, low, upp,
                              method="exact-cond")
t1 = tost_poisson_2indep(count1, n1, count2, n2, low-0.1, upp,
                              method="exact-cond")
t2 = tost_poisson_2indep(count1, n1, count2, n2, low, upp+0.1,
                              method="exact-cond")
for t in [t0, t1, t2]:
    print(t.pvalue, t.results_larger.pvalue, t.results_smaller.pvalue)

0.02499999999999971 (0.02499999999999971, 0.024999999999999818)
0.010095852181662136 (0.010095852181662136, 0.024999999999999818)
0.018672947738695933 (0.02499999999999971, 0.018672947738695933)
```

widening one of the equivalence margins should keep tost pvalue at level 0.025,
i.e. confint is strictly inside equivalence margin at one margin and is equal to other equivalence margin


Discussion:

Patch:
diff --git a/statsmodels/stats/nonparametric.py b/statsmodels/stats/nonparametric.py
index 7d98c12c1..09f4cf1fd 100644
--- a/statsmodels/stats/nonparametric.py
+++ b/statsmodels/stats/nonparametric.py
@@ -220,7 +220,7 @@ class RankCompareResult(HolderTuple):
         t1 = self.test_prob_superior(low, alternative='larger')
         t2 = self.test_prob_superior(upp, alternative='smaller')
 
-        # idx_max = 0 if t1.pvalue < t2.pvalue else 1
+        # idx_max = 1 if t1.pvalue < t2.pvalue else 0
         idx_max = np.asarray(t1.pvalue < t2.pvalue, int)
         title = "Equivalence test for Prob(x1 > x2) + 0.5 Prob(x1 = x2) "
         res = HolderTuple(statistic=np.choose(idx_max,
diff --git a/statsmodels/stats/proportion.py b/statsmodels/stats/proportion.py
index 79af109e9..25528204e 100644
--- a/statsmodels/stats/proportion.py
+++ b/statsmodels/stats/proportion.py
@@ -1442,6 +1442,9 @@ def _shrink_prob(count1, nobs1, count2, nobs2, shrink_factor=2,
         false.
 
     """
+    vectorized = any(np.size(i) > 1 for i in [count1, nobs1, count2, nobs2])
+    if vectorized:
+        raise ValueError("function is not vectorized")
     nobs_col = np.array([count1 + count2, nobs1 - count1 + nobs2 - count2])
     nobs_row = np.array([nobs1, nobs2])
     nobs = nobs1 + nobs2
@@ -1746,6 +1749,9 @@ def test_proportions_2indep(count1, nobs1, count2, nobs2, value=None,
         # TODO: odds ratio does not work if value=1 for score test
         value = 0 if compare == 'diff' else 1
 
+    count1, nobs1, count2, nobs2 = map(np.asarray,
+                                       [count1, nobs1, count2, nobs2])
+
     p1 = count1 / nobs1
     p2 = count2 / nobs2
     diff = p1 - p2
@@ -1980,9 +1986,13 @@ def tost_proportions_2indep(count1, nobs1, count2, nobs2, low, upp,
                                   correction=correction,
                                   return_results=True)
 
-    idx_max = 0 if tt1.pvalue < tt2.pvalue else 1
-    res = HolderTuple(statistic=[tt1.statistic, tt2.statistic][idx_max],
-                      pvalue=[tt1.pvalue, tt2.pvalue][idx_max],
+    # idx_max = 1 if t1.pvalue < t2.pvalue else 0
+    idx_max = np.asarray(tt1.pvalue < tt2.pvalue, int)
+    statistic = np.choose(idx_max, [tt1.statistic, tt2.statistic])
+    pvalue = np.choose(idx_max, [tt1.pvalue, tt2.pvalue])
+
+    res = HolderTuple(statistic=statistic,
+                      pvalue=pvalue,
                       compare=compare,
                       method=method,
                       results_larger=tt1,
diff --git a/statsmodels/stats/rates.py b/statsmodels/stats/rates.py
index 0da5ed1a6..8f7cf5d27 100644
--- a/statsmodels/stats/rates.py
+++ b/statsmodels/stats/rates.py
@@ -1152,9 +1152,13 @@ def tost_poisson_2indep(count1, exposure1, count2, exposure2, low, upp,
                               compare=compare,
                               alternative='smaller')
 
-    idx_max = 0 if tt1.pvalue > tt2.pvalue else 1
-    res = HolderTuple(statistic=[tt1.statistic, tt2.statistic][idx_max],
-                      pvalue=[tt1.pvalue, tt2.pvalue][idx_max],
+    # idx_max = 1 if t1.pvalue < t2.pvalue else 0
+    idx_max = np.asarray(tt1.pvalue < tt2.pvalue, int)
+    statistic = np.choose(idx_max, [tt1.statistic, tt2.statistic])
+    pvalue = np.choose(idx_max, [tt1.pvalue, tt2.pvalue])
+
+    res = HolderTuple(statistic=statistic,
+                      pvalue=pvalue,
                       method=method,
                       compare=compare,
                       equiv_limits=(low, upp),


Test patch:
diff --git a/statsmodels/stats/tests/test_nonparametric.py b/statsmodels/stats/tests/test_nonparametric.py
index c5bac26ae..62ea81083 100644
--- a/statsmodels/stats/tests/test_nonparametric.py
+++ b/statsmodels/stats/tests/test_nonparametric.py
@@ -391,9 +391,15 @@ def test_rank_compare_2indep1():
     assert_allclose(res_ub[1], 0.05, rtol=1e-13)
 
     # test consistency of tost and confint
-    res_tost = res.tost_prob_superior(*ci)
-    assert_allclose(res_tost.results_smaller.pvalue, 0.025, rtol=1e-13)
+    # lower margin is binding, alternative larger
+    res_tost = res.tost_prob_superior(ci[0], ci[1] * 1.05)
     assert_allclose(res_tost.results_larger.pvalue, 0.025, rtol=1e-13)
+    assert_allclose(res_tost.pvalue, 0.025, rtol=1e-13)
+
+    # upper margin is binding, alternative smaller
+    res_tost = res.tost_prob_superior(ci[0] * 0.85, ci[1])
+    assert_allclose(res_tost.results_smaller.pvalue, 0.025, rtol=1e-13)
+    assert_allclose(res_tost.pvalue, 0.025, rtol=1e-13)
 
     # use t-distribution
     # our ranking is defined as reversed from lawstat, and BM article
@@ -419,9 +425,15 @@ def test_rank_compare_2indep1():
     assert_allclose(res_ub[1], 0.05, rtol=1e-11)
 
     # test consistency of tost and confint
-    res_tost = res.tost_prob_superior(*ci)
-    assert_allclose(res_tost.results_smaller.pvalue, 0.025, rtol=1e-11)
-    assert_allclose(res_tost.results_larger.pvalue, 0.025, rtol=1e-11)
+    # lower margin is binding, alternative larger
+    res_tost = res.tost_prob_superior(ci[0], ci[1] * 1.05)
+    assert_allclose(res_tost.results_larger.pvalue, 0.025, rtol=1e-10)
+    assert_allclose(res_tost.pvalue, 0.025, rtol=1e-10)
+
+    # upper margin is binding, alternative smaller
+    res_tost = res.tost_prob_superior(ci[0] * 0.85, ci[1])
+    assert_allclose(res_tost.results_smaller.pvalue, 0.025, rtol=1e-10)
+    assert_allclose(res_tost.pvalue, 0.025, rtol=1e-10)
 
     # extras
     # cohen's d
diff --git a/statsmodels/stats/tests/test_proportion.py b/statsmodels/stats/tests/test_proportion.py
index 6ef65cdd8..fbbe69013 100644
--- a/statsmodels/stats/tests/test_proportion.py
+++ b/statsmodels/stats/tests/test_proportion.py
@@ -424,7 +424,7 @@ def test_binom_tost():
     # vectorized, TODO: observed proportion = 0 returns nan
     ci = smprop.proportion_confint(np.arange(1, 20), 20, method='beta',
                                    alpha=0.05)
-    bt = smprop.binom_tost(np.arange(1, 20), 20, *ci)
+    bt = smprop.binom_tost(np.arange(1, 20), 20, ci[0], ci[1])
     bt = np.asarray(bt)
     assert_almost_equal(bt, 0.025 * np.ones(bt.shape), decimal=12)
 
@@ -816,6 +816,8 @@ def test_equivalence_2indep():
     alpha = 0.05
     count1, nobs1 = 7, 34
     count2, nobs2 = 1, 34
+    count1v, nobs1v = [7, 1], 34
+    count2v, nobs2v = [1, 7], 34
 
     methods_both = [
                     ('diff', 'agresti-caffo'),
@@ -837,11 +839,27 @@ def test_equivalence_2indep():
                                               alpha=2 * alpha,
                                               correction=False)
 
+        # Note: test should have only one margin at confint
         res = smprop.tost_proportions_2indep(
-                count1, nobs1, count2, nobs2, low, upp, compare=co,
+                count1, nobs1, count2, nobs2, low, upp * 1.05, compare=co,
                 method=method, correction=False)
         assert_allclose(res.pvalue, alpha, atol=1e-10)
 
+        res = smprop.tost_proportions_2indep(
+                count1, nobs1, count2, nobs2, low * 0.95, upp, compare=co,
+                method=method, correction=False)
+        assert_allclose(res.pvalue, alpha, atol=1e-10)
+
+        # vectorized
+        if method == 'logit-smoothed':
+            # not correctly vectorized
+            return
+        res1 = res  # for debugging  # noqa
+        res = smprop.tost_proportions_2indep(
+                count1v, nobs1v, count2v, nobs2v, low * 0.95, upp, compare=co,
+                method=method, correction=False)
+        assert_allclose(res.pvalue[0], alpha, atol=1e-10)
+
 
 def test_score_confint_koopman_nam():
 
diff --git a/statsmodels/stats/tests/test_rates_poisson.py b/statsmodels/stats/tests/test_rates_poisson.py
index 0ad227f78..c930bee3a 100644
--- a/statsmodels/stats/tests/test_rates_poisson.py
+++ b/statsmodels/stats/tests/test_rates_poisson.py
@@ -783,6 +783,34 @@ class TestMethodsCompare2indep():
 
             assert_allclose(tst2.pvalue, tst.pvalue, rtol=1e-12)
 
+        # check vectorized equivalence test
+        if compare == "ratio":
+            f = 1.5
+            low, upp = 1 / f, f
+        else:
+            v = 0.5
+            low, upp = -v, v
+
+        tst0 = smr.tost_poisson_2indep(count1[0], n1[0], count2[0], n2[0],
+                                       low, upp,
+                                       method=meth,
+                                       compare=compare)
+        tst1 = smr.tost_poisson_2indep(count1[1], n1[1], count2[1], n2[1],
+                                       low, upp,
+                                       method=meth,
+                                       compare=compare)
+
+        tst2 = smr.tost_poisson_2indep(count1, n1, count2, n2, low, upp,
+                                       method=meth,
+                                       compare=compare)
+        assert tst2.statistic.shape == (2,)
+        assert tst2.pvalue.shape == (2,)
+
+        assert_allclose(tst2.statistic[0], tst0.statistic, rtol=1e-12)
+        assert_allclose(tst2.pvalue[0], tst0.pvalue, rtol=1e-12)
+        assert_allclose(tst2.statistic[1], tst1.statistic, rtol=1e-12)
+        assert_allclose(tst2.pvalue[1], tst1.pvalue, rtol=1e-12)
+
 
 def test_y_grid_regression():
     y_grid = arange(1000)