Repository: statsmodels/statsmodels (Python)

Issue:
Bug in /statsmodels/stats/proportion.py?
#### Describe the bug

In the `_std_2prop_power` function, when [calculating](https://github.com/statsmodels/statsmodels/blob/b722afa1ceb3b90df33491d64dfc319f9ab89a4a/statsmodels/stats/proportion.py#L1894) `std_null`, the parameter `ratio` (sample size ratio, nobs2 = ratio * nobs1) is not correctly passed to `_std_diff_prop`, Instead, default value of one is used. 

This potential bug would impact both `samplesize_proportions_2indep_onetail` and `power_proportions_2indep`, when the two samples are unequally sized. 

#### Code Sample

```python
from statsmodels.stats import proportion

diff, p2, ratio = 0.1, 0.1, 2
p1 = p2 + diff
p_pooled = (p1 + p2 * ratio) / (1 + ratio)
p1_vnull, p2_vnull = p_pooled, p_pooled

# the current behavior
std_null_statsmodels = proportion._std_diff_prop(p1_vnull, p2_vnull)
# the expected behavior
std_null_correct = proportion._std_diff_prop(p1_vnull, p2_vnull, ratio = ratio)
# their difference
print([std_null_statsmodels, std_null_correct])
```
<details>

See formula (4) of this paper for technical details:

Fleiss, Joseph L., Alex Tytun, and Hans K. Ury. A simple approximation for calculating sample sizes for comparing independent proportions. Biometrics (1980): 343-346.

Basically, the current behavior incorrectly assumes `ratio = 1`, for the first term in the numerator. 

</details>


If the issue has not been resolved, please file it in the issue tracker.

#### Expected Output

Change [this](https://github.com/statsmodels/statsmodels/blob/b722afa1ceb3b90df33491d64dfc319f9ab89a4a/statsmodels/stats/proportion.py#L1894) line to `std_null = _std_diff_prop(p1_vnull, p2_vnull, ratio=nobs_ratio)`



Hints:
power_proportions_2indep doesn't propagate any keywords  (although alpha and value are unused arguments in the function)

`p_pooled, std_null, std_alt = _std_2prop_power(diff, prop2, ratio=1,
                                                   alpha=0.05, value=0)`

`std_null = _std_diff_prop(p1_vnull, p2_vnull)`
I need to check the normalization/scaling here (nobs or nobs1)
The variance in both samples are the same in this null case. (but the null could have a nonzero `value`


Thanks for the quick reply. 

In the null case, `p1 = p2 = p_pooled`, so I agree that `p1(1 - p1) = p2(1 - p2) =  p_pooled(1 - p_pooled)`. However, variance of the proportion estimate for sample 2 should have been be scaled by `ratio = nobs2 / nobs1`, which was not done.

If you check the paper I referenced (in details above), or other sample size calculators such as this one (see the `r+1` term for $m^\prime$):

https://www2.ccrb.cuhk.edu.hk/stat/proportion/Casagrande.htm

You'll see what I am talking about.


some more info:

nonzero `value` is not supported, 
Would be better to raise ValueError if value != 0, 
or implement it? (AFAIR, other packages don't allow nonzero value. and we don't have power for tost_proportions_2indep)
```
    value : float
        currently only `value=0`, i.e. equality testing, is supported
```


```
std_null
            standard error of difference under the null hypothesis (without
            sqrt(nobs))
        std_alt
            standard error of difference under the alternative hypothesis
            (without sqrt(nobs))
```
AFAICS, this should be sqrt(nobs1) instead of sqrt(nobs)
i.e. AFAIR everything is normalized by `nobs1` which is also the target of the sample size computation

all unit tests in `test_power_2indep` (against R package ? and Stata use ratio=1)
summ: 2 bugs

- power_proportions_2indep does not propagate ratio to _std_2prop_power, and
- _std_2prop_power does not use ratio in std_null as pointed out in the first message

and clarify docstring for std

@jiannanlu Thanks for finding and reporting the bug

related aside:
#6675 original merged PR
I don't see anything about using the power_proportions_2indep function with the power classes to solve for nobs1 and similar.
Solving for nobs-`ratio`, requires that it is properly propagated during rootfinding. 
Finding effect_size is also "non-standard". We need to make sure std_alt adjusts as we change the alt value (diff) or effect size.

brief search and skimming

It looks like they have the one-tail sample size computation with a non-zero null "value" (either diff or risk ratio)
(I didn't read their discussion for the estimator, e.g. Miettinen and Nurminen MLE under Null diff restriction.)

Farrington, Conor P., and Godfrey Manning. "Test statistics and sample size formulae for comparative binomial trials with null hypothesis of non‐zero risk difference or non‐unity relative risk." Statistics in medicine 9, no. 12 (1990): 1447-1454.

see also https://github.com/statsmodels/statsmodels/issues/4828#issuecomment-592985459

@bashtage This needs to go into the next bugfix release

I can prepare it later today or tomorrow
Please get it in.  It looks like aside from a release note 0.13.2 is ready.

Patch:
diff --git a/statsmodels/stats/proportion.py b/statsmodels/stats/proportion.py
--- a/statsmodels/stats/proportion.py
+++ b/statsmodels/stats/proportion.py
@@ -1891,7 +1891,7 @@ def _std_2prop_power(diff, p2, ratio=1, alpha=0.05, value=0):
     p2_alt = p2
     p1_alt = p2_alt + diff
 
-    std_null = _std_diff_prop(p1_vnull, p2_vnull)
+    std_null = _std_diff_prop(p1_vnull, p2_vnull, ratio=nobs_ratio)
     std_alt = _std_diff_prop(p1_alt, p2_alt, ratio=nobs_ratio)
     return p_pooled, std_null, std_alt
 
@@ -1947,16 +1947,16 @@ def power_proportions_2indep(diff, prop2, nobs1, ratio=1, alpha=0.05,
             pooled proportion, used for std_null
         std_null
             standard error of difference under the null hypothesis (without
-            sqrt(nobs))
+            sqrt(nobs1))
         std_alt
             standard error of difference under the alternative hypothesis
-            (without sqrt(nobs))
+            (without sqrt(nobs1))
     """
     # TODO: avoid possible circular import, check if needed
     from statsmodels.stats.power import normal_power_het
 
-    p_pooled, std_null, std_alt = _std_2prop_power(diff, prop2, ratio=1,
-                                                   alpha=0.05, value=0)
+    p_pooled, std_null, std_alt = _std_2prop_power(diff, prop2, ratio=ratio,
+                                                   alpha=alpha, value=value)
 
     pow_ = normal_power_het(diff, nobs1, alpha, std_null=std_null,
                             std_alternative=std_alt,


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
@@ -915,6 +915,16 @@ def test_power_2indep():
                                                alternative='two-sided')
     assert_allclose(n2, n, rtol=1e-13)
 
+    # with nobs ratio != 1
+    # note Stata has reversed ratio compared to ours, see #8049
+    pwr_st = 0.7995659211532175
+    n = 154
+    res = power_proportions_2indep(-0.1, 0.2, n, ratio=2.)
+    assert_allclose(res.power, pwr_st, atol=1e-7)
+
+    n2 = samplesize_proportions_2indep_onetail(-0.1, 0.2, pwr_st, ratio=2)
+    assert_allclose(n2, n, rtol=1e-4)
+
 
 @pytest.mark.parametrize("count", np.arange(10, 90, 5))
 @pytest.mark.parametrize("method", list(probci_methods.keys()) + ["binom_test"])