Repository: statsmodels/statsmodels (Python)

Issue:
Incorrect comment character in R code
#### Describe the bug

The R script [CH.r](https://github.com/statsmodels/statsmodels/blob/main/statsmodels/stats/libqsturng/CH.r) appears to use `%` as a comment character rather than `#`.  This invalid R code is problematic when also using [`renv`](https://rstudio.github.io/renv/index.html) (an R virtual environment tool, similar to `venv` for Python) for dependency detection.

Changing the comment character to `#` should fix the issue – happy to open a simple PR addressing this if it makes things easier.

#### Code Sample

Steps to reproduce (assuming [R is installed](https://www.r-project.org/))

```bash
# Create venv and install statsmodels
mkdir renv-issue && cd renv-issue
python3 -m venv venv/
source venv/bin/activate
pip install statsmodels

# Install renv R package
R -e "install.packages('remotes', repos='http://cran.us.r-project.org')"
R -e "remotes::install_github('rstudio/renv@0.16.0')"

# initialize renv library (and discover dependencies)
R -e "renv::init()"
```

That last line gives the following error, with unexpected input at the `%` character:

```
WARNING: One or more problems were discovered while enumerating dependencies.

/Users/mike.garcia/renv-issue/venv/lib/python3.9/site-packages/statsmodels/stats/libqsturng/CH.r
------------------------------------------------------------------------------------------------

ERROR 1: /Users/mike.garcia/renv-issue/venv/lib/python3.9/site-packages/statsmodels/stats/libqsturng/CH.r:1:1: unexpected input
1: % Copyright (c) 2011, Roger Lew BSD [see LICENSE.txt]
    ^
```

<details>

This is not directly an issue for using `statsmodels` (as far as I can tell), but is creating some issues when using `renv` in a multi-lingual project that also contains a Python virtual environment containing `statsmodels`. The `renv` package automatically discovers R dependencies in the project directory: when it encounters the invalid R code in `CH.r` it throws an `unexpected input` warning requiring user y/n prompt to proceed.

I realize that `renv` provides a option to use a [`.renvignore` file](https://rstudio.github.io/renv/reference/dependencies.html#ignoring-files) analogous to `.gitignore` to which the virtual environment could (should?) be added. But thought I'd raise the flag on the R syntax here anyway.

- ✅  I did not find any issue referencing `CH.r` in current issues.
- ✅  Issue is present in current version of `main` branch.

</details>

#### Expected Output

`R -e renv::init()` in the above reprex should not display the error in 


#### Output of ``import statsmodels.api as sm; sm.show_versions()``

<details>

```

INSTALLED VERSIONS
------------------
Python: 3.9.5.final.0
OS: Darwin 21.6.0 Darwin Kernel Version 21.6.0: Sat Jun 18 17:07:25 PDT 2022; root:xnu-8020.140.41~1/RELEASE_X86_64 x86_64
byteorder: little
LC_ALL: None
LANG: en_US.UTF-8

statsmodels
===========

Installed: 0.13.5 (/Users/mike.garcia/renv-issue/venv/lib/python3.9/site-packages/statsmodels)

Required Dependencies
=====================

cython: Not installed
numpy: 1.24.1 (/Users/mike.garcia/renv-issue/venv/lib/python3.9/site-packages/numpy)
scipy: 1.10.0 (/Users/mike.garcia/renv-issue/venv/lib/python3.9/site-packages/scipy)
pandas: 1.5.3 (/Users/mike.garcia/renv-issue/venv/lib/python3.9/site-packages/pandas)
    dateutil: 2.8.2 (/Users/mike.garcia/renv-issue/venv/lib/python3.9/site-packages/dateutil)
patsy: 0.5.3 (/Users/mike.garcia/renv-issue/venv/lib/python3.9/site-packages/patsy)

Optional Dependencies
=====================

matplotlib: Not installed
cvxopt: Not installed
joblib: Not installed

Developer Tools
================

IPython: Not installed
    jinja2: Not installed
sphinx: Not installed
    pygments: Not installed
pytest: Not installed
virtualenv: Not installed
```

</details>



Hints:

Patch:
diff --git a/statsmodels/discrete/count_model.py b/statsmodels/discrete/count_model.py
--- a/statsmodels/discrete/count_model.py
+++ b/statsmodels/discrete/count_model.py
@@ -369,7 +369,7 @@ def hessian(self, params):
     def predict(self, params, exog=None, exog_infl=None, exposure=None,
                 offset=None, which='mean', y_values=None):
         """
-        Predict response variable or other statistic given exogenous variables.
+        Predict expected response or other statistic given exogenous variables.
 
         Parameters
         ----------
@@ -396,8 +396,8 @@ def predict(self, params, exog=None, exog_infl=None, exposure=None,
         which : str (optional)
             Statitistic to predict. Default is 'mean'.
 
-            - 'mean' : the conditional expectation of endog E(y | x),
-              i.e. exp of linear predictor.
+            - 'mean' : the conditional expectation of endog E(y | x). This
+              takes inflated zeros into account.
             - 'linear' : the linear predictor of the mean function.
             - 'var' : returns the estimated variance of endog implied by the
               model.
diff --git a/statsmodels/distributions/copula/transforms.py b/statsmodels/distributions/copula/transforms.py
--- a/statsmodels/distributions/copula/transforms.py
+++ b/statsmodels/distributions/copula/transforms.py
@@ -14,7 +14,6 @@
 from scipy.special import expm1, gamma
 
 
-# not used yet
 class Transforms:
 
     def __init__(self):
diff --git a/statsmodels/genmod/families/family.py b/statsmodels/genmod/families/family.py
--- a/statsmodels/genmod/families/family.py
+++ b/statsmodels/genmod/families/family.py
@@ -34,6 +34,10 @@ class Family:
     variance : a variance function
         Measures the variance as a function of the mean probabilities.
         See the individual families for the default variance function.
+    check_link : bool
+        If True (default), then and exception is raised if the link is invalid
+        for the family.
+        If False, then the link is not checked.
 
     See Also
     --------
@@ -383,6 +387,10 @@ class Poisson(Family):
         The default link for the Poisson family is the log link. Available
         links are log, identity, and sqrt. See statsmodels.families.links for
         more information.
+    check_link : bool
+        If True (default), then and exception is raised if the link is invalid
+        for the family.
+        If False, then the link is not checked.
 
     Attributes
     ----------
@@ -535,6 +543,10 @@ class Gaussian(Family):
         The default link for the Gaussian family is the identity link.
         Available links are log, identity, and inverse.
         See statsmodels.genmod.families.links for more information.
+    check_link : bool
+        If True (default), then and exception is raised if the link is invalid
+        for the family.
+        If False, then the link is not checked.
 
     Attributes
     ----------
@@ -705,6 +717,10 @@ class Gamma(Family):
         The default link for the Gamma family is the inverse link.
         Available links are log, identity, and inverse.
         See statsmodels.genmod.families.links for more information.
+    check_link : bool
+        If True (default), then and exception is raised if the link is invalid
+        for the family.
+        If False, then the link is not checked.
 
     Attributes
     ----------
@@ -867,6 +883,10 @@ class Binomial(Family):
         The default link for the Binomial family is the logit link.
         Available links are logit, probit, cauchy, log, loglog, and cloglog.
         See statsmodels.genmod.families.links for more information.
+    check_link : bool
+        If True (default), then and exception is raised if the link is invalid
+        for the family.
+        If False, then the link is not checked.
 
     Attributes
     ----------
@@ -1136,6 +1156,10 @@ class InverseGaussian(Family):
         inverse squared link.
         Available links are InverseSquared, Inverse, Log, and Identity.
         See statsmodels.genmod.families.links for more information.
+    check_link : bool
+        If True (default), then and exception is raised if the link is invalid
+        for the family.
+        If False, then the link is not checked.
 
     Attributes
     ----------
@@ -1299,6 +1323,10 @@ class NegativeBinomial(Family):
         The ancillary parameter for the negative binomial distribution.
         For now ``alpha`` is assumed to be nonstochastic.  The default value
         is 1.  Permissible values are usually assumed to be between .01 and 2.
+    check_link : bool
+        If True (default), then and exception is raised if the link is invalid
+        for the family.
+        If False, then the link is not checked.
 
     Attributes
     ----------
@@ -1518,6 +1546,10 @@ class Tweedie(Family):
         likelihood is used.
         In both cases, for likelihood computations the var_power
         must be between 1 and 2.
+    check_link : bool
+        If True (default), then and exception is raised if the link is invalid
+        for the family.
+        If False, then the link is not checked.
 
     Attributes
     ----------
diff --git a/statsmodels/genmod/generalized_linear_model.py b/statsmodels/genmod/generalized_linear_model.py
--- a/statsmodels/genmod/generalized_linear_model.py
+++ b/statsmodels/genmod/generalized_linear_model.py
@@ -218,6 +218,10 @@ class GLM(base.LikelihoodModel):
 
     Notes
     -----
+    Note: PerfectSeparationError exception has been converted to a
+    PerfectSeparationWarning and perfect separation or perfect prediction will
+    not raise an exception by default. (changed in version 0.14)
+
     Only the following combinations make sense for family and link:
 
      ============= ===== === ===== ====== ======= === ==== ====== ====== ====
diff --git a/statsmodels/robust/norms.py b/statsmodels/robust/norms.py
--- a/statsmodels/robust/norms.py
+++ b/statsmodels/robust/norms.py
@@ -831,10 +831,12 @@ class MQuantileNorm(RobustNorm):
     in QuantileRegression but replaces the L1 absolute value by a chosen
     base norm.
 
-    rho_q(u) = |q - I(q < 0)| * rho_base(u)
-    or
-    rho_q(u) = q * rho_base(u)  if u >= 0
-    rho_q(u) = (1 - q) * rho_base(u)  if u < 0
+        rho_q(u) = |q - I(q < 0)| * rho_base(u)
+
+    or, equivalently,
+
+        rho_q(u) = q * rho_base(u)  if u >= 0
+        rho_q(u) = (1 - q) * rho_base(u)  if u < 0
 
 
     Parameters
@@ -854,9 +856,22 @@ class MQuantileNorm(RobustNorm):
 
     References
     ----------
-    Newey Powell
-    Jones
-    Bianchi Salvati
+
+    .. [*] Bianchi, Annamaria, and Nicola Salvati. 2015. “Asymptotic Properties
+       and Variance Estimators of the M-Quantile Regression Coefficients
+       Estimators.” Communications in Statistics - Theory and Methods 44 (11):
+       2416–29. doi:10.1080/03610926.2013.791375.
+
+    .. [*] Breckling, Jens, and Ray Chambers. 1988. “M-Quantiles.”
+       Biometrika 75 (4): 761–71. doi:10.2307/2336317.
+
+    .. [*] Jones, M. C. 1994. “Expectiles and M-Quantiles Are Quantiles.”
+       Statistics & Probability Letters 20 (2): 149–53.
+       doi:10.1016/0167-7152(94)90031-0.
+
+    .. [*] Newey, Whitney K., and James L. Powell. 1987. “Asymmetric Least
+       Squares Estimation and Testing.” Econometrica 55 (4): 819–47.
+       doi:10.2307/1911031.
     ...
 
     """
diff --git a/statsmodels/stats/power.py b/statsmodels/stats/power.py
--- a/statsmodels/stats/power.py
+++ b/statsmodels/stats/power.py
@@ -30,6 +30,8 @@
 
 
 """
+import warnings
+
 import numpy as np
 from scipy import stats, optimize
 from statsmodels.tools.rootfinding import brentq_expanding
@@ -232,7 +234,7 @@ def ftest_anova_power(effect_size, nobs, alpha, k_groups=2, df=None):
     return pow_#, crit
 
 
-def ftest_power(effect_size, df_num, df_denom, alpha, ncc=1):
+def ftest_power(effect_size, df2, df1, alpha, ncc=1):
     '''Calculate the power of a F-test.
 
     Parameters
@@ -240,10 +242,12 @@ def ftest_power(effect_size, df_num, df_denom, alpha, ncc=1):
     effect_size : float
         standardized effect size, mean divided by the standard deviation.
         effect size has to be positive.
-    df_num : int or float
-        numerator degrees of freedom.
-    df_denom : int or float
-        denominator degrees of freedom.
+    df2 : int or float
+        Denominator degrees of freedom.
+        This corresponds to the df_resid in Wald tests.
+    df1 : int or float
+        Numerator degrees of freedom.
+        This corresponds to the number of constraints 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.
@@ -260,18 +264,25 @@ def ftest_power(effect_size, df_num, df_denom, alpha, ncc=1):
 
     Notes
     -----
+    changed in 0.14: use df2, df1 instead of df_num, df_denom as arg names.
+    The latter had reversed meaning.
+
+    The sample size is given implicitly by ``df2`` with fixed number of
+    constraints given by numerator degrees of freedom ``df1``:
 
-    sample size is given implicitly by df_num
+        nobs = df2 + df1 + ncc
 
-    set ncc=0 to match t-test, or f-test in LikelihoodModelResults.
+    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)
+    models, with df_num (df1) as number of constraints and d_denom (df2) as
+    df_resid.
     '''
+    df_num, df_denom = df1, df2
     nc = effect_size**2 * (df_denom + df_num + ncc)
-    crit = stats.f.isf(alpha, df_denom, df_num)
-    pow_ = stats.ncf.sf(crit, df_denom, df_num, nc)
+    crit = stats.f.isf(alpha, df_num, df_denom)
+    pow_ = stats.ncf.sf(crit, df_num, df_denom, nc)
     return pow_ #, crit, nc
 
 
@@ -581,8 +592,8 @@ def solve_power(self, effect_size=None, nobs=None, alpha=None, power=None,
         Parameters
         ----------
         effect_size : float
-            standardized effect size, mean divided by the standard deviation.
-            effect size has to be positive.
+            Standardized effect size.The effect size is here Cohen's f, square
+            root of "f2".
         nobs : int or float
             sample size, number of observations.
         alpha : float in interval (0,1)
@@ -868,22 +879,63 @@ def solve_power(self, effect_size=None, nobs1=None, alpha=None, power=None,
 
 
 class FTestPower(Power):
-    '''Statistical Power calculations for generic F-test
+    """Statistical Power calculations for generic F-test of a constraint
 
-    '''
+    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.
+
+    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 df, wrongly named ``df_num``
+
+    >>> df1 = 1  # number of constraints in hypothesis test
+    >>> df2 = FTestPower().solve_power(effect_size=f, alpha=0.1, power=0.9,
+                                       df_denom=df1)
+    >>> ncc = 1  # default
+    >>> nobs = df2 + df1 + ncc
+    >>> df2, nobs
+    (76.46459758305376, 78.46459758305376)
+
+    verify power at df2
+
+    >>> FTestPower().power(effect_size=f, alpha=0.1, df_denom=df1, df_num=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``, square root of ``f2``.
+
+        The sample size is given by ``nobs = df_denom + df_num + ncc``
+
+        Warning: The meaning of df_num and df_denom is reversed.
+
         Parameters
         ----------
         effect_size : float
-            standardized effect size, mean divided by the standard deviation.
-            effect size has to be positive.
+            Standardized effect size. The effect size is here Cohen's ``f``,
+            square root of ``f2``.
         df_num : int or float
-            numerator degrees of freedom.
+            Warning incorrect name
+            denominator degrees of freedom,
+            This corresponds to the number of constraints in Wald tests.
         df_denom : int or float
-            denominator degrees of freedom.
+            Warning incorrect name
+            numerator 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.
@@ -916,7 +968,7 @@ def power(self, effect_size, df_num, df_denom, alpha, ncc=1):
 
     #method is only added to have explicit keywords and docstring
     def solve_power(self, effect_size=None, df_num=None, df_denom=None,
-                    nobs=None, alpha=None, power=None, ncc=1):
+                    alpha=None, power=None, ncc=1, **kwargs):
         '''solve for any one parameter of the power of a F-test
 
         for the one sample F-test the keywords are:
@@ -924,14 +976,26 @@ def solve_power(self, effect_size=None, df_num=None, df_denom=None,
 
         Exactly one needs to be ``None``, all others need numeric values.
 
+        The effect size is Cohen's ``f``, square root of ``f2``.
+
+        The sample size is given by ``nobs = df_denom + df_num + ncc``.
+
+        Warning: The meaning of df_num and df_denom is reversed.
 
         Parameters
         ----------
         effect_size : float
-            standardized effect size, mean divided by the standard deviation.
-            effect size has to be positive.
-        nobs : int or float
-            sample size, number of observations.
+            Standardized effect size. The effect size is here Cohen's ``f``,
+            square root of ``f2``.
+        df_num : int or float
+            Warning incorrect name
+            denominator degrees of freedom,
+            This corresponds to the number of constraints in Wald tests.
+            Sample size is given by ``nobs = df_denom + df_num + ncc``
+        df_denom : int or float
+            Warning incorrect name
+            numerator 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.
@@ -939,10 +1003,13 @@ def solve_power(self, effect_size=None, df_num=None, df_denom=None,
             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.
-        alternative : str, 'two-sided' (default) or 'one-sided'
-            extra argument to choose whether the power is calculated for a
-            two-sided (default) or one sided test.
-            'one-sided' assumes we are in the relevant tail.
+        ncc : int
+            degrees of freedom correction for non-centrality parameter.
+            see Notes
+        kwargs : empty
+            ``kwargs`` are not used and included for backwards compatibility.
+            If ``nobs`` is used as keyword, then a warning is issued. All
+            other keywords in ``kwargs`` raise a ValueError.
 
         Returns
         -------
@@ -961,6 +1028,11 @@ def solve_power(self, effect_size=None, df_num=None, df_denom=None,
         where this fails.
 
         '''
+        if kwargs:
+            if "nobs" in kwargs:
+                warnings.warn("nobs is not used")
+            else:
+                raise ValueError(f"incorrect keyword(s) {kwargs}")
         return super(FTestPower, self).solve_power(effect_size=effect_size,
                                                       df_num=df_num,
                                                       df_denom=df_denom,
@@ -971,6 +1043,12 @@ def solve_power(self, effect_size=None, df_num=None, df_denom=None,
 class FTestAnovaPower(Power):
     '''Statistical Power calculations F-test for one factor balanced ANOVA
 
+    This is based on Cohen's f as effect size measure.
+
+    See Also
+    --------
+    statsmodels.stats.oneway.effectsize_oneway
+
     '''
 
     def power(self, effect_size, nobs, alpha, k_groups=2):
@@ -979,8 +1057,8 @@ def power(self, effect_size, nobs, alpha, k_groups=2):
         Parameters
         ----------
         effect_size : float
-            standardized effect size, mean divided by the standard deviation.
-            effect size has to be positive.
+            standardized effect size. The effect size is here Cohen's f, square
+            root of "f2".
         nobs : int or float
             sample size, number of observations.
         alpha : float in interval (0,1)
diff --git a/statsmodels/tsa/stattools.py b/statsmodels/tsa/stattools.py
--- a/statsmodels/tsa/stattools.py
+++ b/statsmodels/tsa/stattools.py
@@ -1472,9 +1472,11 @@ def grangercausalitytests(x, maxlag, addconst=True, verbose=None):
     >>> data = data.data[["realgdp", "realcons"]].pct_change().dropna()
 
     All lags up to 4
+
     >>> gc_res = grangercausalitytests(data, 4)
 
     Only lag 4
+
     >>> gc_res = grangercausalitytests(data, [4])
     """
     x = array_like(x, "x", ndim=2)


Test patch:
diff --git a/statsmodels/stats/tests/test_power.py b/statsmodels/stats/tests/test_power.py
--- a/statsmodels/stats/tests/test_power.py
+++ b/statsmodels/stats/tests/test_power.py
@@ -688,6 +688,19 @@ def setup_class(cls):
         # precision for test_power
         cls.decimal = 5
 
+    def test_kwargs(self):
+
+        with pytest.warns(UserWarning):
+            smp.FTestPower().solve_power(
+                effect_size=0.3, alpha=0.1, power=0.9, df_denom=2,
+                nobs=None)
+
+        with pytest.raises(ValueError):
+            smp.FTestPower().solve_power(
+                effect_size=0.3, alpha=0.1, power=0.9, df_denom=2,
+                junk=3)
+
+
 
 def test_power_solver():
     # messing up the solver to trigger backup