Repository: statsmodels/statsmodels

Issue:
BUG: FTestPower always Fails to converge on a solution, wrong names, unclear usage
**edit** (josef)
FTestPower mixes up df_num, df_denom names of keywords, meaning is reversed. 
Code itself works correctly but not obvious how to use it.


#### Describe the bug

I am trying to do a power analysis for regression - compute the number of necessary observations.

No matter what parameters I supply, I get 
```
python3.10/site-packages/statsmodels/stats/power.py:415: ConvergenceWarning: 
Failed to converge on a solution.
```
and the return value is `array([3.])`.

1. Why does it return an array?
2. Why does it always return the same value?

#### Code Sample, a copy-pastable example if possible

```python
import statsmodels.api as sm
analysis = sm.stats.FTestPower()
analysis.solve_power(nobs=None, effect_size=1, alpha=0.1, power=0.9, df_num=1)
```

I assume that the effect size is `f^2=r^2/(1-r^2)` ([Cohen's f^2](https://en.wikipedia.org/wiki/Effect_size#Cohen's_f2)).
This is a simple regression, so `df=1`.

#### Expected Output

I expect a single number returned - at least for _some_ inputs, no warnings.

#### Output of ``import statsmodels.api as sm; sm.show_versions()``

<details>

INSTALLED VERSIONS
------------------
Python: 3.10.6.final.0
OS: Linux 6.0.6-76060006-generic #202210290932~1669062050~22.04~d94609a SMP PREEMPT_DYNAMIC Mon N x86_64
byteorder: little
LC_ALL: None
LANG: C.UTF-8

statsmodels
===========

Installed: 0.13.5 (/home/sds/.virtualenvs/myself/lib/python3.10/site-packages/statsmodels)

Required Dependencies
=====================

cython: Not installed
numpy: 1.24.1 (/home/sds/.virtualenvs/myself/lib/python3.10/site-packages/numpy)
scipy: 1.10.0 (/home/sds/.virtualenvs/myself/lib/python3.10/site-packages/scipy)
pandas: 1.5.3 (/home/sds/.virtualenvs/myself/lib/python3.10/site-packages/pandas)
    dateutil: 2.8.2 (/home/sds/.virtualenvs/myself/lib/python3.10/site-packages/dateutil)
patsy: 0.5.3 (/home/sds/.virtualenvs/myself/lib/python3.10/site-packages/patsy)

Optional Dependencies
=====================

matplotlib: 3.6.3 (/home/sds/.virtualenvs/myself/lib/python3.10/site-packages/matplotlib)
    backend: module://matplotlib_inline.backend_inline 
cvxopt: Not installed
joblib: 1.2.0 (/home/sds/.virtualenvs/myself/lib/python3.10/site-packages/joblib)

Developer Tools
================

IPython: 8.9.0 (/home/sds/.virtualenvs/myself/lib/python3.10/site-packages/IPython)
    jinja2: 3.1.2 (/home/sds/.virtualenvs/myself/lib/python3.10/site-packages/jinja2)
sphinx: Not installed
    pygments: 2.14.0 (/home/sds/.virtualenvs/myself/lib/python3.10/site-packages/pygments)
pytest: 7.2.1 (/home/sds/.virtualenvs/myself/lib/python3.10/site-packages/pytest)
virtualenv: 20.17.1 (/home/sds/.virtualenvs/myself/lib/python3.10/site-packages/virtualenv)

</details>


Discussion:

Patch:
diff --git a/statsmodels/stats/oneway.py b/statsmodels/stats/oneway.py
index 42d9dd8ed..04dbd1d97 100644
--- a/statsmodels/stats/oneway.py
+++ b/statsmodels/stats/oneway.py
@@ -11,6 +11,9 @@ import numpy as np
 from scipy import stats
 from scipy.special import ncfdtrinc
 
+# functions that use scipy.special instead of boost based function in stats
+from statsmodels.stats.power import ncf_cdf, ncf_ppf
+
 from statsmodels.stats.robust_compare import TrimmedMean, scale_transform
 from statsmodels.tools.testing import Holder
 from statsmodels.stats.base import HolderTuple
@@ -774,7 +777,7 @@ def equivalence_oneway_generic(f_stat, n_groups, nobs, equiv_margin, df,
         type_effectsize = "Cohen's f_squared"
     else:
         raise ValueError('`margin_type` should be "f2" or "wellek"')
-    crit_f = stats.ncf.ppf(alpha, df[0], df[1], nc_null)
+    crit_f = ncf_ppf(alpha, df[0], df[1], nc_null)
 
     if margin_type == "wellek":
         # TODO: do we need a sqrt
@@ -784,8 +787,8 @@ def equivalence_oneway_generic(f_stat, n_groups, nobs, equiv_margin, df,
 
     reject = (es < crit_es)
 
-    pv = stats.ncf.cdf(f_stat, df[0], df[1], nc_null)
-    pwr = stats.ncf.cdf(crit_f, df[0], df[1], 1e-13)  # scipy, cannot be 0
+    pv = ncf_cdf(f_stat, df[0], df[1], nc_null)
+    pwr = ncf_cdf(crit_f, df[0], df[1], 1e-13)  # scipy, cannot be 0
     res = HolderTuple(statistic=f_stat,
                       pvalue=pv,
                       effectsize=es,  # match es type to margin_type
@@ -916,7 +919,7 @@ def _power_equivalence_oneway_emp(f_stat, n_groups, nobs, eps, df, alpha=0.05):
     nobs_mean = nobs.sum() / n_groups
     fn = f_stat  # post-hoc power, empirical power at estimate
     esn = fn * (n_groups - 1) / nobs_mean  # Wellek psi
-    pow_ = stats.ncf.cdf(res.crit_f, df[0], df[1], nobs_mean * esn)
+    pow_ = ncf_cdf(res.crit_f, df[0], df[1], nobs_mean * esn)
 
     return pow_
 
@@ -980,8 +983,8 @@ def power_equivalence_oneway(f2_alt, equiv_margin, nobs_t, n_groups=None,
     else:
         raise ValueError('`margin_type` should be "f2" or "wellek"')
 
-    crit_f_margin = stats.ncf.ppf(alpha, df[0], df[1], nobs_t * f2_null)
-    pwr_alt = stats.ncf.cdf(crit_f_margin, df[0], df[1], nobs_t * f2_alt)
+    crit_f_margin = ncf_ppf(alpha, df[0], df[1], nobs_t * f2_null)
+    pwr_alt = ncf_cdf(crit_f_margin, df[0], df[1], nobs_t * f2_alt)
     return pwr_alt
 
 
diff --git a/statsmodels/stats/power.py b/statsmodels/stats/power.py
index bfa1cd508..c50d832ab 100644
--- a/statsmodels/stats/power.py
+++ b/statsmodels/stats/power.py
@@ -33,9 +33,30 @@ refactoring
 import warnings
 
 import numpy as np
-from scipy import stats, optimize
+from scipy import stats, optimize, special
 from statsmodels.tools.rootfinding import brentq_expanding
 
+
+def nct_cdf(x, df, nc):
+    return special.nctdtr(df, nc, x)
+
+
+def nct_sf(x, df, nc):
+    return 1 - special.nctdtr(df, nc, x)
+
+
+def ncf_cdf(x, dfn, dfd, nc):
+    return special.ncfdtr(dfn, dfd, nc, x)
+
+
+def ncf_sf(x, dfn, dfd, nc):
+    return 1 - special.ncfdtr(dfn, dfd, nc, x)
+
+
+def ncf_ppf(q, dfn, dfd, nc):
+    return special.ncfdtri(dfn, dfd, nc, q)
+
+
 def ttest_power(effect_size, nobs, alpha, df=None, alternative='two-sided'):
     '''Calculate power of a ttest
     '''
@@ -60,16 +81,20 @@ def ttest_power(effect_size, nobs, alpha, df=None, alternative='two-sided'):
             # avoid endless loop, https://github.com/scipy/scipy/issues/2667
             pow_ = np.nan
         else:
-            pow_ = stats.nct._sf(crit_upp, df, d*np.sqrt(nobs))
+            # pow_ = stats.nct._sf(crit_upp, df, d*np.sqrt(nobs))
+            # use scipy.special
+            pow_ = nct_sf(crit_upp, df, d*np.sqrt(nobs))
     if alternative in ['two-sided', '2s', 'smaller']:
         crit_low = stats.t.ppf(alpha_, df)
         #print crit_low, df, d*np.sqrt(nobs)
         if np.any(np.isnan(crit_low)):
             pow_ = np.nan
         else:
-            pow_ += stats.nct._cdf(crit_low, df, d*np.sqrt(nobs))
+            # pow_ += stats.nct._cdf(crit_low, df, d*np.sqrt(nobs))
+            pow_ += nct_cdf(crit_low, df, d*np.sqrt(nobs))
     return pow_
 
+
 def normal_power(effect_size, nobs, alpha, alternative='two-sided', sigma=1.):
     """Calculate power of a normal distributed test statistic
 
@@ -227,11 +252,11 @@ def ftest_anova_power(effect_size, nobs, alpha, k_groups=2, df=None):
 
     should be general nobs observations, k_groups restrictions ???
     '''
-    df_num = nobs - k_groups
-    df_denom = k_groups - 1
-    crit = stats.f.isf(alpha, df_denom, df_num)
-    pow_ = stats.ncf.sf(crit, df_denom, df_num, effect_size**2 * nobs)
-    return pow_#, crit
+    df_num = k_groups - 1
+    df_denom = nobs - k_groups
+    crit = stats.f.isf(alpha, df_num, df_denom)
+    pow_ = ncf_sf(crit, df_num, df_denom, effect_size**2 * nobs)
+    return pow_
 
 
 def ftest_power(effect_size, df2, df1, alpha, ncc=1):
@@ -240,8 +265,7 @@ def ftest_power(effect_size, df2, df1, alpha, ncc=1):
     Parameters
     ----------
     effect_size : float
-        standardized effect size, mean divided by the standard deviation.
-        effect size has to be positive.
+        The effect size is here Cohen's ``f``, the square root of ``f2``.
     df2 : int or float
         Denominator degrees of freedom.
         This corresponds to the df_resid in Wald tests.
@@ -282,10 +306,73 @@ def ftest_power(effect_size, df2, df1, alpha, ncc=1):
     df_num, df_denom = df1, df2
     nc = effect_size**2 * (df_denom + df_num + ncc)
     crit = stats.f.isf(alpha, df_num, df_denom)
-    pow_ = stats.ncf.sf(crit, df_num, df_denom, nc)
+    # pow_ = stats.ncf.sf(crit, df_num, df_denom, nc)
+    # use scipy.special for ncf
+    pow_ = ncf_sf(crit, df_num, df_denom, nc)
     return pow_ #, crit, nc
 
 
+def ftest_power_f2(effect_size, df_num, df_denom, alpha, ncc=1):
+    '''Calculate the power of a F-test.
+
+    Based on Cohen's `f^2` effect size.
+
+    This assumes
+
+        df_num : numerator degrees of freedom, (number of constraints)
+        df_denom : denominator degrees of freedom (df_resid in regression)
+        nobs = df_denom + df_num + ncc
+        nc = effect_size * nobs  (noncentrality index)
+
+    Power is computed one-sided in the upper tail.
+
+    Parameters
+    ----------
+    effect_size : float
+        Cohen's f2 effect size or noncentrality divided by nobs.
+    df_num : int or float
+        Numerator degrees of freedom.
+        This corresponds to the number of constraints in Wald tests.
+    df_denom : int or float
+        Denominator degrees of freedom.
+        This corresponds to the df_resid in Wald tests.
+    alpha : float in interval (0,1)
+        significance level, e.g. 0.05, is the probability of a type I
+        error, that is wrong rejections if the Null Hypothesis is true.
+    ncc : int
+        degrees of freedom correction for non-centrality parameter.
+        see Notes
+
+    Returns
+    -------
+    power : float
+        Power of the test, e.g. 0.8, is one minus the probability of a
+        type II error. Power is the probability that the test correctly
+        rejects the Null Hypothesis if the Alternative Hypothesis is true.
+
+    Notes
+
+    The sample size is given implicitly by ``df_denom`` with fixed number of
+    constraints given by numerator degrees of freedom ``df_num``:
+
+        nobs = df_denom + df_num + ncc
+
+    Set ncc=0 to match t-test, or f-test in LikelihoodModelResults.
+    ncc=1 matches the non-centrality parameter in R::pwr::pwr.f2.test
+
+    ftest_power with ncc=0 should also be correct for f_test in regression
+    models, with df_num (df1) as number of constraints and d_denom (df2) as
+    df_resid.
+    '''
+
+    nc = effect_size * (df_denom + df_num + ncc)
+    crit = stats.f.isf(alpha, df_num, df_denom)
+    # pow_ = stats.ncf.sf(crit, df_num, df_denom, nc)
+    # use scipy.special for ncf
+    pow_ = ncf_sf(crit, df_num, df_denom, nc)
+    return pow_
+
+
 #class based implementation
 #--------------------------
 
@@ -881,10 +968,17 @@ class NormalIndPower(Power):
 class FTestPower(Power):
     """Statistical Power calculations for generic F-test of a constraint
 
+    This class is not recommended, use `FTestPowerF2` with corrected interface.
+
     This is based on Cohen's f as effect size measure.
 
     Warning: Methods in this class have the names df_num and df_denom reversed.
 
+    See Also
+    --------
+    FTestPowerF2 :
+        Class with Cohen's f-squared as effect size, corrected keyword names.
+
     Examples
     --------
     Sample size and power for multiple regression base on R-squared
@@ -1020,7 +1114,7 @@ class FTestPower(Power):
 
         Notes
         -----
-        The function uses scipy.optimize for finding the value that satisfies
+        The method uses scipy.optimize for finding the value that satisfies
         the power equation. It first uses ``brentq`` with a prior search for
         bounds. If this fails to find a root, ``fsolve`` is used. If ``fsolve``
         also fails, then, for ``alpha``, ``power`` and ``effect_size``,
@@ -1040,6 +1134,151 @@ class FTestPower(Power):
                                                       power=power,
                                                       ncc=ncc)
 
+
+
+class FTestPowerF2(Power):
+    """Statistical Power calculations for generic F-test of a constraint
+
+    This is based on Cohen's f^2 as effect size measure.
+
+    Examples
+    --------
+    Sample size and power for multiple regression base on R-squared
+
+    Compute effect size from R-squared
+
+    >>> r2 = 0.1
+    >>> f2 = r2 / (1 - r2)
+    >>> f = np.sqrt(f2)
+    >>> r2, f2, f
+    (0.1, 0.11111111111111112, 0.33333333333333337)
+
+    Find sample size by solving for denominator degrees of freedom.
+
+    >>> df1 = 1  # number of constraints in hypothesis test
+    >>> df2 = FTestPowerF2().solve_power(effect_size=f2, alpha=0.1, power=0.9,
+                                         df_num=df1)
+    >>> ncc = 1  # default
+    >>> nobs = df2 + df1 + ncc
+    >>> df2, nobs
+    (76.46459758305376, 78.46459758305376)
+
+    verify power at df2
+
+    >>> FTestPowerF2().power(effect_size=f, alpha=0.1, df_num=df1, df_denom=df2)
+    0.8999999972109698
+
+    """
+
+    def power(self, effect_size, df_num, df_denom, alpha, ncc=1):
+        '''Calculate the power of a F-test.
+
+        The effect size is Cohen's ``f^2``.
+
+        The sample size is given by ``nobs = df_denom + df_num + ncc``
+
+        Parameters
+        ----------
+        effect_size : float
+            The effect size is here Cohen's ``f2``. This is equal to
+            the noncentrality of an F-test divided by nobs.
+        df_num : int or float
+            Numerator degrees of freedom,
+            This corresponds to the number of constraints in Wald tests.
+        df_denom : int or float
+            Denominator degrees of freedom.
+            This corresponds to the df_resid in Wald tests.
+        alpha : float in interval (0,1)
+            Significance level, e.g. 0.05, is the probability of a type I
+            error, that is wrong rejections if the Null Hypothesis is true.
+        ncc : int
+            Degrees of freedom correction for non-centrality parameter.
+            see Notes
+
+        Returns
+        -------
+        power : float
+            Power of the test, e.g. 0.8, is one minus the probability of a
+            type II error. Power is the probability that the test correctly
+            rejects the Null Hypothesis if the Alternative Hypothesis is true.
+
+        Notes
+        -----
+        The sample size is given implicitly by df_denom
+
+        set ncc=0 to match t-test, or f-test in LikelihoodModelResults.
+        ncc=1 matches the non-centrality parameter in R::pwr::pwr.f2.test
+
+        ftest_power with ncc=0 should also be correct for f_test in regression
+        models, with df_num and d_denom as defined there. (not verified yet)
+        '''
+
+        pow_ = ftest_power_f2(effect_size, df_num, df_denom, alpha, ncc=ncc)
+        return pow_
+
+    #method is only added to have explicit keywords and docstring
+    def solve_power(self, effect_size=None, df_num=None, df_denom=None,
+                    alpha=None, power=None, ncc=1):
+        '''Solve for any one parameter of the power of a F-test
+
+        for the one sample F-test the keywords are:
+            effect_size, df_num, df_denom, alpha, power
+
+        Exactly one needs to be ``None``, all others need numeric values.
+
+        The effect size is Cohen's ``f2``.
+
+        The sample size is given by ``nobs = df_denom + df_num + ncc``, and
+        can be found by solving for df_denom.
+
+        Parameters
+        ----------
+        effect_size : float
+            The effect size is here Cohen's ``f2``. This is equal to
+            the noncentrality of an F-test divided by nobs.
+        df_num : int or float
+            Numerator degrees of freedom,
+            This corresponds to the number of constraints in Wald tests.
+        df_denom : int or float
+            Denominator degrees of freedom.
+            This corresponds to the df_resid in Wald tests.
+        alpha : float in interval (0,1)
+            significance level, e.g. 0.05, is the probability of a type I
+            error, that is wrong rejections if the Null Hypothesis is true.
+        power : float in interval (0,1)
+            power of the test, e.g. 0.8, is one minus the probability of a
+            type II error. Power is the probability that the test correctly
+            rejects the Null Hypothesis if the Alternative Hypothesis is true.
+        ncc : int
+            degrees of freedom correction for non-centrality parameter.
+            see Notes
+
+        Returns
+        -------
+        value : float
+            The value of the parameter that was set to None in the call. The
+            value solves the power equation given the remaining parameters.
+
+
+        Notes
+        -----
+        The function uses scipy.optimize for finding the value that satisfies
+        the power equation. It first uses ``brentq`` with a prior search for
+        bounds. If this fails to find a root, ``fsolve`` is used. If ``fsolve``
+        also fails, then, for ``alpha``, ``power`` and ``effect_size``,
+        ``brentq`` with fixed bounds is used. However, there can still be cases
+        where this fails.
+
+        '''
+
+        return super(FTestPowerF2, self).solve_power(effect_size=effect_size,
+                                                      df_num=df_num,
+                                                      df_denom=df_denom,
+                                                      alpha=alpha,
+                                                      power=power,
+                                                      ncc=ncc)
+
+
 class FTestAnovaPower(Power):
     '''Statistical Power calculations F-test for one factor balanced ANOVA
 


Test patch:
diff --git a/statsmodels/stats/tests/test_power.py b/statsmodels/stats/tests/test_power.py
index a6e94df9c..154e26c4a 100644
--- a/statsmodels/stats/tests/test_power.py
+++ b/statsmodels/stats/tests/test_power.py
@@ -88,7 +88,7 @@ class CheckPowerMixin:
 
     @pytest.mark.matplotlib
     def test_power_plot(self, close_figures):
-        if self.cls == smp.FTestPower:
+        if self.cls in [smp.FTestPower, smp.FTestPowerF2]:
             pytest.skip('skip FTestPower plot_power')
         fig = plt.figure()
         ax = fig.add_subplot(2,1,1)
@@ -701,6 +701,33 @@ class TestFtestPower(CheckPowerMixin):
                 junk=3)
 
 
+class TestFtestPowerF2(CheckPowerMixin):
+
+    @classmethod
+    def setup_class(cls):
+        res2 = Holder()
+        #> rf = pwr.f2.test(u=5, v=19, f2=0.3**2, sig.level=0.1)
+        #> cat_items(rf, "res2.")
+        res2.u = 5
+        res2.v = 19
+        res2.f2 = 0.09
+        res2.sig_level = 0.1
+        res2.power = 0.235454222377575
+        res2.method = 'Multiple regression power calculation'
+
+        cls.res2 = res2
+        cls.kwds = {'effect_size': res2.f2, 'df_num': res2.u,
+                     'df_denom': res2.v, 'alpha': res2.sig_level,
+                     'power': res2.power}
+        # keyword for which we do not look for root:
+        # solving for n_bins does not work, will not be used in regular usage
+        cls.kwds_extra = {}
+        cls.args_names = ['effect_size', 'df_num', 'df_denom', 'alpha']
+        cls.cls = smp.FTestPowerF2
+        # precision for test_power
+        cls.decimal = 5
+
+
 
 def test_power_solver():
     # messing up the solver to trigger backup