Repository: statsmodels/statsmodels (Python)

Issue:
BUG: GLS fit_regularized, different params if sigma is 1-dim or diag 2 dim
I don't see the problem right now. Needs investigating.

Fixing a unit test for GLS.fit_regularized that used incorrect kwarg "weights" #7751

I think the results, specifically params, should be the same if `sigma = 1/ wghts` and `sigma = np.diag(1 / wgths)`
The second matches the params from the `WLS(... weights=wgths)` estimate.
The former looks like an unpenalized estimate


my current version form test_regression
the last two asserts fail if I use `model_gls = GLS(yvec, xmat, sigma=1 / wgt)

```
def test_regularized_predict():
    n = 100
    p = 5
    np.random.seed(3132)
    xmat = np.random.normal(size=(n, p))
    yvec = xmat.sum(1) + np.random.normal(size=n)
    wgt = np.random.uniform(1, 2, n)
    model_wls = WLS(yvec, xmat, weights=wgt)
    # TODO: params is not the same in GLS if sigma=1 / wgth, i.e 1-dim
    model_gls = GLS(yvec, xmat, sigma=np.diag(1 / wgt))
    res = []
    for model1 in [model_wls, model_gls]:
        result1 = model1.fit_regularized(alpha=20.0, L1_wt=0.5, refit=True)
        res.append(result1)
        params = result1.params
        fittedvalues = np.dot(xmat, params)
        pr = model1.predict(result1.params)
        assert_allclose(fittedvalues, pr)
        assert_allclose(result1.fittedvalues, pr)

        pr = result1.predict()
        assert_allclose(fittedvalues, pr)

    assert_allclose(res[0].model.wendog, res[1].model.wendog, rtol=1e-10)
    assert_allclose(res[0].model.wexog, res[1].model.wexog, rtol=1e-10)
    assert_allclose(res[0].fittedvalues, res[1].fittedvalues, rtol=1e-10)
    assert_allclose(res[0].params, res[1].params, rtol=1e-10)
```


Hints:

Patch:
diff --git a/statsmodels/base/_penalized.py b/statsmodels/base/_penalized.py
--- a/statsmodels/base/_penalized.py
+++ b/statsmodels/base/_penalized.py
@@ -51,7 +51,7 @@ def __init__(self, *args, **kwds):
 
         self._init_keys.extend(['penal', 'pen_weight'])
         self._null_drop_keys = getattr(self, '_null_drop_keys', [])
-        self._null_drop_keys.extend(['penal'])
+        self._null_drop_keys.extend(['penal', 'pen_weight'])
 
     def _handle_scale(self, params, scale=None, **kwds):
 
diff --git a/statsmodels/base/model.py b/statsmodels/base/model.py
--- a/statsmodels/base/model.py
+++ b/statsmodels/base/model.py
@@ -29,6 +29,8 @@
 from statsmodels.tools.tools import nan_dot, recipr
 from statsmodels.tools.validation import bool_like
 
+ERROR_INIT_KWARGS = False
+
 _model_params_doc = """Parameters
     ----------
     endog : array_like
@@ -79,6 +81,10 @@ class Model(object):
     # Default is 1, which is more common. Override in models when needed
     # Set to None to skip check
     _formula_max_endog = 1
+    # kwargs that are generically allowed, maybe not supported in all models
+    _kwargs_allowed = [
+        "missing", 'missing_idx', 'formula', 'design_info', "hasconst",
+        ]
 
     def __init__(self, endog, exog=None, **kwargs):
         missing = kwargs.pop('missing', 'none')
@@ -106,6 +112,22 @@ def _get_init_kwds(self):
 
         return kwds
 
+    def _check_kwargs(self, kwargs, keys_extra=None, error=ERROR_INIT_KWARGS):
+
+        kwargs_allowed = [
+            "missing", 'missing_idx', 'formula', 'design_info', "hasconst",
+            ]
+        if keys_extra:
+            kwargs_allowed.extend(keys_extra)
+
+        kwargs_invalid = [i for i in kwargs if i not in kwargs_allowed]
+        if kwargs_invalid:
+            msg = "unknown kwargs " + repr(kwargs_invalid)
+            if error is False:
+                warnings.warn(msg, ValueWarning)
+            else:
+                raise ValueError(msg)
+
     def _handle_data(self, endog, exog, missing, hasconst, **kwargs):
         data = handle_data(endog, exog, missing, hasconst, **kwargs)
         # kwargs arrays could have changed, easier to just attach here
diff --git a/statsmodels/discrete/discrete_model.py b/statsmodels/discrete/discrete_model.py
--- a/statsmodels/discrete/discrete_model.py
+++ b/statsmodels/discrete/discrete_model.py
@@ -457,8 +457,9 @@ class BinaryModel(DiscreteModel):
     _continuous_ok = False
 
     def __init__(self, endog, exog, check_rank=True, **kwargs):
+        # unconditional check, requires no extra kwargs added by subclasses
+        self._check_kwargs(kwargs)
         super().__init__(endog, exog, check_rank, **kwargs)
-
         if not issubclass(self.__class__, MultinomialModel):
             if not np.all((self.endog >= 0) & (self.endog <= 1)):
                 raise ValueError("endog must be in the unit interval.")
@@ -774,6 +775,7 @@ def _derivative_exog(self, params, exog=None, transform='dydx',
 class CountModel(DiscreteModel):
     def __init__(self, endog, exog, offset=None, exposure=None, missing='none',
                  check_rank=True, **kwargs):
+        self._check_kwargs(kwargs)
         super().__init__(endog, exog, check_rank, missing=missing,
                          offset=offset, exposure=exposure, **kwargs)
         if exposure is not None:
diff --git a/statsmodels/distributions/copula/archimedean.py b/statsmodels/distributions/copula/archimedean.py
--- a/statsmodels/distributions/copula/archimedean.py
+++ b/statsmodels/distributions/copula/archimedean.py
@@ -168,6 +168,7 @@ def pdf(self, u, args=()):
         return a * b ** c
 
     def logpdf(self, u, args=()):
+        # we skip Archimedean logpdf, that uses numdiff
         return super(ArchimedeanCopula, self).logpdf(u, args=args)
 
     def cdf(self, u, args=()):
@@ -268,6 +269,7 @@ def logpdf(self, u, args=()):
             return pdf
         else:
             # for now use generic from base Copula class, log(self.pdf(...))
+            # we skip Archimedean logpdf, that uses numdiff
             super(ArchimedeanCopula, self).logpdf(u, args)
 
     def cdfcond_2g1(self, u, args=()):
@@ -383,6 +385,7 @@ def cdf(self, u, args=()):
         return cdf
 
     def logpdf(self, u, args=()):
+        # we skip Archimedean logpdf, that uses numdiff
         return super(ArchimedeanCopula, self).logpdf(u, args=args)
 
     def tau(self, theta=None):
diff --git a/statsmodels/genmod/generalized_estimating_equations.py b/statsmodels/genmod/generalized_estimating_equations.py
--- a/statsmodels/genmod/generalized_estimating_equations.py
+++ b/statsmodels/genmod/generalized_estimating_equations.py
@@ -496,9 +496,10 @@ def __init__(self, endog, exog, groups, time=None, family=None,
                  exposure=None, dep_data=None, constraint=None,
                  update_dep=True, weights=None, **kwargs):
 
+        if type(self) is GEE:
+            self._check_kwargs(kwargs)
         if family is not None:
             if not isinstance(family.link, tuple(family.safe_links)):
-                import warnings
                 msg = ("The {0} link function does not respect the "
                        "domain of the {1} family.")
                 warnings.warn(msg.format(family.link.__class__.__name__,
@@ -543,6 +544,12 @@ def __init__(self, endog, exog, groups, time=None, family=None,
 
         self._init_keys.extend(["update_dep", "constraint", "family",
                                 "cov_struct"])
+        # remove keys added by super that are not supported
+        try:
+            self._init_keys.remove("freq_weights")
+            self._init_keys.remove("var_weights")
+        except ValueError:
+            pass
 
         # Handle the family argument
         if family is None:
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
@@ -289,14 +289,8 @@ def __init__(self, endog, exog, family=None, offset=None,
                  exposure=None, freq_weights=None, var_weights=None,
                  missing='none', **kwargs):
 
-        # TODO: this can be applied to many models
-        # 'n_trials' is specific to GLM Binomial
-        kwargs_allowed = ["missing", 'missing_idx', 'formula', 'design_info',
-                          "hasconst", 'n_trials']
-        kwargs_invalid = [i for i in kwargs if i not in kwargs_allowed]
-        if kwargs_invalid and type(self) is GLM:
-            warnings.warn("unknown kwargs" + repr(kwargs_invalid),
-                          UserWarning)
+        if type(self) is GLM:
+            self._check_kwargs(kwargs, ['n_trials'])
 
         if (family is not None) and not isinstance(family.link,
                                                    tuple(family.safe_links)):
@@ -1665,8 +1659,10 @@ def null(self):
         model = self.model
         exog = np.ones((len(endog), 1))
 
-        kwargs = model._get_init_kwds()
+        kwargs = model._get_init_kwds().copy()
         kwargs.pop('family')
+        for key in getattr(model, '_null_drop_keys', []):
+            del kwargs[key]
         start_params = np.atleast_1d(self.family.link(endog.mean()))
         oe = self.model._offset_exposure
         if not (np.size(oe) == 1 and oe == 0):
diff --git a/statsmodels/regression/linear_model.py b/statsmodels/regression/linear_model.py
--- a/statsmodels/regression/linear_model.py
+++ b/statsmodels/regression/linear_model.py
@@ -47,7 +47,10 @@
 # need import in module instead of lazily to copy `__doc__`
 from statsmodels.regression._prediction import PredictionResults
 from statsmodels.tools.decorators import cache_readonly, cache_writable
-from statsmodels.tools.sm_exceptions import InvalidTestWarning
+from statsmodels.tools.sm_exceptions import (
+    InvalidTestWarning,
+    ValueWarning,
+    )
 from statsmodels.tools.tools import pinv_extended
 from statsmodels.tools.validation import string_like
 
@@ -498,6 +501,8 @@ class GLS(RegressionModel):
 
     def __init__(self, endog, exog, sigma=None, missing='none', hasconst=None,
                  **kwargs):
+        if type(self) is GLS:
+            self._check_kwargs(kwargs)
         # TODO: add options igls, for iterative fgls if sigma is None
         # TODO: default if sigma is none should be two-step GLS
         sigma, cholsigmainv = _get_sigma(sigma, len(endog))
@@ -620,7 +625,14 @@ def fit_regularized(self, method="elastic_net", alpha=0.,
         # Need to adjust since RSS/n term in elastic net uses nominal
         # n in denominator
         if self.sigma is not None:
-            alpha = alpha * np.sum(1 / np.diag(self.sigma)) / len(self.endog)
+            if self.sigma.ndim == 2:
+                var_obs = np.diag(self.sigma)
+            elif self.sigma.ndim == 1:
+                var_obs = self.sigma
+            else:
+                raise ValueError("sigma should be 1-dim or 2-dim")
+
+            alpha = alpha * np.sum(1 / var_obs) / len(self.endog)
 
         rslt = OLS(self.wendog, self.wexog).fit_regularized(
             method=method, alpha=alpha,
@@ -689,6 +701,8 @@ class WLS(RegressionModel):
 
     def __init__(self, endog, exog, weights=1., missing='none', hasconst=None,
                  **kwargs):
+        if type(self) is WLS:
+            self._check_kwargs(kwargs)
         weights = np.array(weights)
         if weights.shape == ():
             if (missing == 'drop' and 'missing_idx' in kwargs and
@@ -871,11 +885,18 @@ class OLS(WLS):
 
     def __init__(self, endog, exog=None, missing='none', hasconst=None,
                  **kwargs):
+        if "weights" in kwargs:
+            msg = ("Weights are not supported in OLS and will be ignored"
+                   "An exception will be raised in the next version.")
+            warnings.warn(msg, ValueWarning)
         super(OLS, self).__init__(endog, exog, missing=missing,
                                   hasconst=hasconst, **kwargs)
         if "weights" in self._init_keys:
             self._init_keys.remove("weights")
 
+        if type(self) is OLS:
+            self._check_kwargs(kwargs, ["offset"])
+
     def loglike(self, params, scale=None):
         """
         The likelihood function for the OLS model.
diff --git a/statsmodels/regression/quantile_regression.py b/statsmodels/regression/quantile_regression.py
--- a/statsmodels/regression/quantile_regression.py
+++ b/statsmodels/regression/quantile_regression.py
@@ -75,6 +75,7 @@ class QuantReg(RegressionModel):
     '''
 
     def __init__(self, endog, exog, **kwargs):
+        self._check_kwargs(kwargs)
         super(QuantReg, self).__init__(endog, exog, **kwargs)
 
     def whiten(self, data):
diff --git a/statsmodels/robust/robust_linear_model.py b/statsmodels/robust/robust_linear_model.py
--- a/statsmodels/robust/robust_linear_model.py
+++ b/statsmodels/robust/robust_linear_model.py
@@ -107,6 +107,7 @@ class RLM(base.LikelihoodModel):
 
     def __init__(self, endog, exog, M=None, missing='none',
                  **kwargs):
+        self._check_kwargs(kwargs)
         self.M = M if M is not None else norms.HuberT()
         super(base.LikelihoodModel, self).__init__(endog, exog,
                                                    missing=missing, **kwargs)


Test patch:
diff --git a/statsmodels/discrete/tests/test_discrete.py b/statsmodels/discrete/tests/test_discrete.py
--- a/statsmodels/discrete/tests/test_discrete.py
+++ b/statsmodels/discrete/tests/test_discrete.py
@@ -45,6 +45,7 @@
     ConvergenceWarning,
     PerfectSeparationError,
     SpecificationWarning,
+    ValueWarning,
 )
 
 from .results.results_discrete import Anes, DiscreteL1, RandHIE, Spector
@@ -392,12 +393,13 @@ def setup_class(cls):
         res2 = Spector.probit
         cls.res2 = res2
 
-    @pytest.mark.xfail(reason="res2 has no predict attribute",
-                       raises=AttributeError, strict=True)
-    def test_predict(self):
-        assert_almost_equal(self.res1.model.predict(self.res1.params),
-                            self.res2.predict,
-                            DECIMAL_4)
+    def test_init_kwargs(self):
+        endog = self.res1.model.endog
+        exog = self.res1.model.exog
+        z = np.ones(len(endog))
+        with pytest.warns(ValueWarning, match="unknown kwargs"):
+            # unsupported keyword
+            Probit(endog, exog, weights=z)
 
 
 class TestProbitBFGS(CheckBinaryResults):
diff --git a/statsmodels/gam/tests/test_penalized.py b/statsmodels/gam/tests/test_penalized.py
--- a/statsmodels/gam/tests/test_penalized.py
+++ b/statsmodels/gam/tests/test_penalized.py
@@ -12,6 +12,8 @@
 from numpy.testing import assert_allclose, assert_equal, assert_
 import pandas as pd
 
+import pytest
+
 import patsy
 
 from statsmodels.discrete.discrete_model import Poisson, Logit, Probit
@@ -91,6 +93,10 @@ def test_fitted(self):
         assert_allclose(res1.fittedvalues, res2.fitted_values,
                         rtol=self.rtol_fitted)
 
+    @pytest.mark.smoke
+    def test_null_smoke(self):
+        self.res1.llnull
+
 
 class TestTheilPLS5(CheckGAMMixin):
 
@@ -120,6 +126,9 @@ def test_cov_robust(self):
         assert_allclose(np.asarray(res1.cov_params()),
                         res2.Ve * self.covp_corrfact, rtol=1e-4)
 
+    def test_null_smoke(self):
+        pytest.skip("llnull not available")
+
 
 class TestGLMPenalizedPLS5(CheckGAMMixin):
 
diff --git a/statsmodels/genmod/tests/test_gee.py b/statsmodels/genmod/tests/test_gee.py
--- a/statsmodels/genmod/tests/test_gee.py
+++ b/statsmodels/genmod/tests/test_gee.py
@@ -1803,7 +1803,6 @@ def test_plots(close_figures):
 
     model = gee.GEE(exog, endog, groups)
     result = model.fit()
-
     fig = result.plot_added_variable(1)
     assert_equal(isinstance(fig, plt.Figure), True)
     fig = result.plot_partial_residuals(1)
@@ -1988,16 +1987,16 @@ def test_quasipoisson(reg):
     grp = np.kron(np.arange(100), np.ones(n // 100))
 
     model1 = gee.GEE(y, x, family=families.Poisson(), groups=grp,
-                     cov_type="naive")
+                     )
     model2 = gee.GEE(y, x, family=families.Poisson(), groups=grp,
-                     cov_type="naive")
+                     )
 
     if reg:
         result1 = model1.fit_regularized(pen_wt=0.1)
         result2 = model2.fit_regularized(pen_wt=0.1, scale="X2")
     else:
-        result1 = model1.fit()
-        result2 = model2.fit(scale="X2")
+        result1 = model1.fit(cov_type="naive")
+        result2 = model2.fit(scale="X2", cov_type="naive")
 
     # The parameter estimates are the same regardless of how
     # the scale parameter is handled
diff --git a/statsmodels/genmod/tests/test_glm.py b/statsmodels/genmod/tests/test_glm.py
--- a/statsmodels/genmod/tests/test_glm.py
+++ b/statsmodels/genmod/tests/test_glm.py
@@ -12,7 +12,6 @@
     assert_array_less,
     assert_equal,
     assert_raises,
-    assert_warns,
 )
 import pandas as pd
 from pandas.testing import assert_series_equal
@@ -32,6 +31,7 @@
 from statsmodels.tools.sm_exceptions import (
     DomainWarning,
     PerfectSeparationError,
+    ValueWarning,
 )
 from statsmodels.tools.tools import add_constant
 
@@ -1736,7 +1736,8 @@ def setup_class(cls):
 
     def test_init_kwargs(self):
         family_link = sm.families.Gamma(sm.families.links.log())
-        with assert_warns(UserWarning):
+
+        with pytest.warns(ValueWarning, match="unknown kwargs"):
             GLM(self.endog, self.exog, family=family_link,
                 weights=self.weight,  # incorrect keyword
                 )
diff --git a/statsmodels/regression/tests/test_predict.py b/statsmodels/regression/tests/test_predict.py
--- a/statsmodels/regression/tests/test_predict.py
+++ b/statsmodels/regression/tests/test_predict.py
@@ -264,7 +264,7 @@ def test_predict_remove_data():
     # GH6887
     endog = [i + np.random.normal(scale=0.1) for i in range(100)]
     exog = [i for i in range(100)]
-    model = OLS(endog, exog, weights=[1 for _ in range(100)]).fit()
+    model = WLS(endog, exog, weights=[1 for _ in range(100)]).fit()
     # we need to compute scale before we remove wendog, wexog
     model.scale
     model.remove_data()
diff --git a/statsmodels/regression/tests/test_regression.py b/statsmodels/regression/tests/test_regression.py
--- a/statsmodels/regression/tests/test_regression.py
+++ b/statsmodels/regression/tests/test_regression.py
@@ -753,6 +753,7 @@ def setup_class(cls):
         cls.results = []
         cls.results.append(OLS(y, x).fit())
         cls.results.append(WLS(y, x, w).fit())
+        # scaling weights does not change main results (except scale)
         cls.results.append(GLS(y, x, 100 * w).fit())
         cls.results.append(GLS(y, x, np.diag(0.1 * w)).fit())
 
@@ -800,6 +801,7 @@ def setup_class(cls):
         w_inv = 1.0 / w
         cls.results = []
         cls.results.append(WLS(y, x, w).fit())
+        # scaling weights does not change main results (except scale)
         cls.results.append(WLS(y, x, 0.01 * w).fit())
         cls.results.append(GLS(y, x, 100 * w_inv).fit())
         cls.results.append(GLS(y, x, np.diag(0.1 * w_inv)).fit())
@@ -1433,17 +1435,21 @@ def test_regularized_refit():
 
 
 def test_regularized_predict():
+    # this also compares WLS with GLS
     n = 100
     p = 5
     np.random.seed(3132)
     xmat = np.random.normal(size=(n, p))
     yvec = xmat.sum(1) + np.random.normal(size=n)
     wgt = np.random.uniform(1, 2, n)
-
-    for klass in WLS, GLS:
-        model1 = klass(yvec, xmat, weights=wgt)
-        result1 = model1.fit_regularized(alpha=2.0, L1_wt=0.5, refit=True)
-
+    model_wls = WLS(yvec, xmat, weights=wgt)
+    # TODO: params is not the same in GLS if sigma=1 / wgt, i.e 1-dim, #7755
+    model_gls1 = GLS(yvec, xmat, sigma=np.diag(1 / wgt))
+    model_gls2 = GLS(yvec, xmat, sigma=1 / wgt)
+    res = []
+    for model1 in [model_wls, model_gls1, model_gls2]:
+        result1 = model1.fit_regularized(alpha=20.0, L1_wt=0.5, refit=True)
+        res.append(result1)
         params = result1.params
         fittedvalues = np.dot(xmat, params)
         pr = model1.predict(result1.params)
@@ -1453,6 +1459,16 @@ def test_regularized_predict():
         pr = result1.predict()
         assert_allclose(fittedvalues, pr)
 
+    assert_allclose(res[0].model.wendog, res[1].model.wendog, rtol=1e-10)
+    assert_allclose(res[0].model.wexog, res[1].model.wexog, rtol=1e-10)
+    assert_allclose(res[0].fittedvalues, res[1].fittedvalues, rtol=1e-10)
+    assert_allclose(res[0].params, res[1].params, rtol=1e-10)
+
+    assert_allclose(res[0].model.wendog, res[2].model.wendog, rtol=1e-10)
+    assert_allclose(res[0].model.wexog, res[2].model.wexog, rtol=1e-10)
+    assert_allclose(res[0].fittedvalues, res[2].fittedvalues, rtol=1e-10)
+    assert_allclose(res[0].params, res[2].params, rtol=1e-10)
+
 
 def test_regularized_options():
     n = 100