Repository: statsmodels/statsmodels (Python)

Issue:
ENH: proportion_confint only has two-sided confidence intervals, no "alternative" option
see
https://stats.stackexchange.com/questions/646445/multinomial-one-sided-confidence-intervals

same for multinomial_proportions_confint

I think for proportion_confint it should be just one interval at 2*alpha, except for "beta-test" which is minlike, and the one-sided will correspond to one side of equal-tail intervals.

I don't have a guess on one-sided interval for multinomial_proportions_confint. (too long ago). But it might be based on an equal-tail derivation.



Hints:
@josef-pkt 
Hi, Josef.
I read the thread and agree that the one-sided interval is the same as one side of the two-sided interval for doubled alpha in a binomial distribution.
Can I take this issue on propotion_confint to add an option to explicitly get one-sided ci?
Yes, we want to add the `alternative` option and include one sided intervals.

Note, "binomtest" will have different one-sided intervals which should be easy to compute.
The will be the same as for the central, equal-tail binom test.

pull request would be welcome.


Patch:
diff --git a/statsmodels/stats/proportion.py b/statsmodels/stats/proportion.py
--- a/statsmodels/stats/proportion.py
+++ b/statsmodels/stats/proportion.py
@@ -108,7 +108,7 @@ def _bisection_search_conservative(
     return best_pt, best
 
 
-def proportion_confint(count, nobs, alpha:float=0.05, method="normal"):
+def proportion_confint(count, nobs, alpha:float=0.05, method="normal", alternative:str='two-sided'):
     """
     Confidence interval for a binomial proportion
 
@@ -133,11 +133,16 @@ def proportion_confint(count, nobs, alpha:float=0.05, method="normal"):
          - `jeffreys` : Jeffreys Bayesian Interval
          - `binom_test` : Numerical inversion of binom_test
 
+    alternative : {"two-sided", "larger", "smaller"}
+        default: "two-sided"
+        specifies whether to calculate a two-sided or one-sided confidence interval.
+
     Returns
     -------
     ci_low, ci_upp : {float, ndarray, Series DataFrame}
-        lower and upper confidence level with coverage (approximately) 1-alpha.
+        larger and smaller confidence level with coverage (approximately) 1-alpha.
         When a pandas object is returned, then the index is taken from `count`.
+        When side is not "two-sided", lower or upper bound is set to 0 or 1 respectively.
 
     Notes
     -----
@@ -188,21 +193,24 @@ def _check(x: np.ndarray, name: str) -> np.ndarray:
         nobs_a = _check(np.asarray(nobs_a), "count")
 
     q_ = count_a / nobs_a
-    alpha_2 = 0.5 * alpha
+
+
+    if alternative == 'two-sided':
+        if method != "binom_test":
+            alpha = alpha / 2.0
+    elif alternative not in ['larger', 'smaller']:
+        raise NotImplementedError(f"alternative {alternative} is not available")
 
     if method == "normal":
         std_ = np.sqrt(q_ * (1 - q_) / nobs_a)
-        dist = stats.norm.isf(alpha / 2.0) * std_
+        dist = stats.norm.isf(alpha) * std_
         ci_low = q_ - dist
         ci_upp = q_ + dist
-    elif method == "binom_test":
-        # inverting the binomial test
+    elif method == "binom_test" and alternative == 'two-sided':
         def func_factory(count: int, nobs: int) -> Callable[[float], float]:
             if hasattr(stats, "binomtest"):
-
                 def func(qi):
                     return stats.binomtest(count, nobs, p=qi).pvalue - alpha
-
             else:
                 # Remove after min SciPy >= 1.7
                 def func(qi):
@@ -258,9 +266,9 @@ def func(qi):
                 ci_upp.flat[index] = 1 - ci_low.flat[index]
                 ci_low.flat[index] = 1 - temp
             index = bcast.index
-    elif method == "beta":
-        ci_low = stats.beta.ppf(alpha_2, count_a, nobs_a - count_a + 1)
-        ci_upp = stats.beta.isf(alpha_2, count_a + 1, nobs_a - count_a)
+    elif method == "beta" or (method == "binom_test" and alternative != 'two-sided'):
+        ci_low = stats.beta.ppf(alpha, count_a, nobs_a - count_a + 1)
+        ci_upp = stats.beta.isf(alpha, count_a + 1, nobs_a - count_a)
 
         if np.ndim(ci_low) > 0:
             ci_low.flat[q_.flat == 0] = 0
@@ -269,7 +277,7 @@ def func(qi):
             ci_low = 0 if q_ == 0 else ci_low
             ci_upp = 1 if q_ == 1 else ci_upp
     elif method == "agresti_coull":
-        crit = stats.norm.isf(alpha / 2.0)
+        crit = stats.norm.isf(alpha)
         nobs_c = nobs_a + crit**2
         q_c = (count_a + crit**2 / 2.0) / nobs_c
         std_c = np.sqrt(q_c * (1.0 - q_c) / nobs_c)
@@ -277,7 +285,7 @@ def func(qi):
         ci_low = q_c - dist
         ci_upp = q_c + dist
     elif method == "wilson":
-        crit = stats.norm.isf(alpha / 2.0)
+        crit = stats.norm.isf(alpha)
         crit2 = crit**2
         denom = 1 + crit2 / nobs_a
         center = (q_ + crit2 / (2 * nobs_a)) / denom
@@ -289,9 +297,8 @@ def func(qi):
         ci_upp = center + dist
     # method adjusted to be more forgiving of misspellings or incorrect option name
     elif method[:4] == "jeff":
-        ci_low, ci_upp = stats.beta.interval(
-            1 - alpha, count_a + 0.5, nobs_a - count_a + 0.5
-        )
+        ci_low = stats.beta.ppf(alpha, count_a + 0.5, nobs_a - count_a + 0.5)
+        ci_upp = stats.beta.isf(alpha, count_a + 0.5, nobs_a - count_a + 0.5)
     else:
         raise NotImplementedError(f"method {method} is not available")
     if method in ["normal", "agresti_coull"]:
@@ -301,6 +308,10 @@ def func(qi):
         container = pd.Series if isinstance(count, pd.Series) else pd.DataFrame
         ci_low = container(ci_low, index=count.index)
         ci_upp = container(ci_upp, index=count.index)
+    if alternative == 'larger':
+        ci_low = 0
+    elif alternative == 'smaller':
+        ci_upp = 1
     if is_scalar:
         return float(ci_low), float(ci_upp)
     return ci_low, ci_upp


Test patch:
diff --git a/statsmodels/stats/tests/test_proportion.py b/statsmodels/stats/tests/test_proportion.py
--- a/statsmodels/stats/tests/test_proportion.py
+++ b/statsmodels/stats/tests/test_proportion.py
@@ -995,3 +995,67 @@ def test_ci_symmetry_array(count, method):
     a = proportion_confint([count, count], n, method=method)
     b = proportion_confint([n - count, n - count], n, method=method)
     assert_allclose(np.array(a), 1.0 - np.array(b[::-1]))
+
+
+@pytest.mark.parametrize("count, nobs, alternative, expected", [
+    (1, 1, 'two-sided', [1.0, 1.0]),
+    (0, 1, 'two-sided', [0.0, 0.0]),
+    (1, 1, 'larger',   [0.0, 1.0]),
+    (0, 1, 'larger',   [0.0, 0.0]),
+    (1, 1, 'smaller',    [1.0, 1.0]),
+    (0, 1, 'smaller',    [0.0, 1.0]),
+])
+def test_proportion_confint_edge(count, nobs, alternative, expected):
+    # test against R 4.3.3 / DescTools / BinomCI and binom.test
+    ci_low, ci_upp = proportion_confint(count, nobs, alternative=alternative)
+    assert_allclose(np.array([ci_low, ci_upp]), np.array(expected))
+
+
+@pytest.mark.parametrize("count", [100])
+@pytest.mark.parametrize("nobs", [200])
+@pytest.mark.parametrize("method, alpha, alternative, expected", [
+    # two-sided
+    ("normal",          0.05, "two-sided",  [0.4307048, 0.5692952]),
+    ("normal",          0.01, "two-sided",  [0.4089307, 0.5910693]),
+    ("normal",          0.10, "two-sided",  [0.4418456, 0.5581544]),
+    ("agresti_coull",   0.05, "two-sided",  [0.4313609, 0.5686391]),
+    ("beta",            0.05, "two-sided",  [0.4286584, 0.5713416]),
+    ("wilson",          0.05, "two-sided",  [0.4313609, 0.5686391]),
+    ("jeffreys",        0.05, "two-sided",  [0.4311217, 0.5688783]),
+    # larger
+    ("normal",          0.05, "larger",    [0.0, 0.5581544]),
+    ("agresti_coull",   0.05, "larger",    [0.0, 0.557765]),
+    ("beta",            0.05, "larger",    [0.0, 0.560359]),
+    ("wilson",          0.05, "larger",    [0.0, 0.557765]),
+    ("jeffreys",        0.05, "larger",    [0.0, 0.5578862]),
+    # smaller
+    ("normal",          0.05, "smaller",     [0.4418456, 1.0]),
+    ("agresti_coull",   0.05, "smaller",     [0.442235, 1.0]),
+    ("beta",            0.05, "smaller",     [0.439641, 1.0]),
+    ("wilson",          0.05, "smaller",     [0.442235, 1.0]),
+    ("jeffreys",        0.05, "smaller",     [0.4421138, 1.0]),
+])
+def test_proportion_confint(count, nobs, method, alpha, alternative, expected):
+    # test against R 4.3.3 / DescTools / BinomCI and binom.test
+    ci_low, ci_upp = proportion_confint(count, nobs, alpha=alpha, method=method, alternative=alternative)
+    assert_allclose(np.array([ci_low, ci_upp]), np.array(expected))
+
+
+@pytest.mark.parametrize("count", [100])
+@pytest.mark.parametrize("nobs", [200])
+@pytest.mark.parametrize("alpha, alternative, expected, rtol", [
+    # binom.exact with method="two.sided" returns values with a precision of up to 4 decimal places.
+    (0.01, "two-sided",  [0.4094, 0.5906],  1e-4),
+    (0.05, "two-sided",  [0.4299, 0.5701],  1e-4),
+    (0.10, "two-sided",  [0.44, 0.56],      1e-4),
+    (0.01, "larger",    [0.0, 0.5840449],  1e-7),
+    (0.05, "larger",    [0.0, 0.560359],   1e-7),
+    (0.10, "larger",    [0.0, 0.5476422],  1e-7),
+    (0.01, "smaller",     [0.4159551, 1.0],  1e-7),
+    (0.05, "smaller",     [0.439641, 1.0],   1e-7),
+    (0.10, "smaller",     [0.4523578, 1.0],  1e-7),
+])
+def test_proportion_confint_binom_test(count, nobs, alpha, alternative, expected, rtol):
+    # test against R 4.3.3 / DescTools / exactci
+    ci_low, ci_upp = proportion_confint(count, nobs, alpha=alpha, method="binom_test", alternative=alternative)
+    assert_allclose(np.array([ci_low, ci_upp]), np.array(expected), rtol=rtol)