Merge pull request statsmodels#5956 from bashtage/var-cov-params

BUG: Fix mutidimensional model cov_params when using pandas

Author: bashtage

statsmodels/base/data.py

@@ -55,6 +55,7 @@ class ModelData(object):
     appropriate form
     """
     _param_names = None
+    _cov_names = None
 
     def __init__(self, endog, exog=None, missing='none', hasconst=None,
                  **kwargs):
@@ -352,6 +353,26 @@ def param_names(self):
     def param_names(self, values):
         self._param_names = values
 
+    @property
+    def cov_names(self):
+        """
+        Labels for covariance matrices
+
+        In multidimensional models, each dimension of a covariance matrix
+        differs from the number of param_names.
+
+        If not set, returns param_names
+        """
+        # for handling names of covariance names in multidimensional models
+        if self._cov_names is not None:
+            return self._cov_names
+        return self.param_names
+
+    @cov_names.setter
+    def cov_names(self, value):
+        # for handling names of covariance names in multidimensional models
+        self._cov_names = value
+
     @cache_readonly
     def row_labels(self):
         exog = self.orig_exog
@@ -427,6 +448,8 @@ def wrap_output(self, obj, how='columns', names=None):
             return self.attach_generic_columns_2d(obj, names)
         elif how == 'ynames':
             return self.attach_ynames(obj)
+        elif how == 'multivariate_confint':
+            return self.attach_mv_confint(obj)
         else:
             return obj
 
@@ -448,6 +471,9 @@ def attach_rows(self, result):
     def attach_dates(self, result):
         return result
 
+    def attach_mv_confint(self, result):
+        return result
+
     def attach_generic_columns(self, result, *args, **kwargs):
         return result
 
@@ -533,8 +559,7 @@ def attach_columns_eq(self, result):
         return DataFrame(result, index=self.xnames, columns=self.ynames)
 
     def attach_cov(self, result):
-        return DataFrame(result, index=self.param_names,
-                         columns=self.param_names)
+        return DataFrame(result, index=self.cov_names, columns=self.cov_names)
 
     def attach_cov_eq(self, result):
         return DataFrame(result, index=self.ynames, columns=self.ynames)
@@ -565,6 +590,11 @@ def attach_dates(self, result):
             return DataFrame(result, index=self.predict_dates,
                              columns=self.ynames)
 
+    def attach_mv_confint(self, result):
+        return DataFrame(result.reshape((-1, 2)),
+                         index=self.cov_names,
+                         columns=['lower', 'upper'])
+
     def attach_ynames(self, result):
         squeezed = result.squeeze()
         # May be zero-dim, for example in the case of forecast one step in tsa

statsmodels/discrete/discrete_model.py

@@ -16,12 +16,12 @@
 W. Greene. `Econometric Analysis`. Prentice Hall, 5th. edition. 2003.
 """
 __all__ = ["Poisson", "Logit", "Probit", "MNLogit", "NegativeBinomial",
-           "GeneralizedPoisson", "NegativeBinomialP"]
+           "GeneralizedPoisson", "NegativeBinomialP", "CountModel"]
 
 from scipy.special import loggamma
 
 import numpy as np
-from pandas import get_dummies
+from pandas import get_dummies, MultiIndex
 
 from scipy.special import gammaln, digamma, polygamma
 from scipy import stats, special
@@ -2162,7 +2162,20 @@ class MNLogit(MultinomialModel):
     Notes
     -----
     See developer notes for further information on `MNLogit` internals.
-    """ % {'extra_params' : base._missing_param_doc}
+    """ % {'extra_params': base._missing_param_doc}
+
+    def __init__(self, endog, exog, **kwargs):
+        super(MNLogit, self).__init__(endog, exog, **kwargs)
+
+        # Override cov_names since multivariate model
+        yname = self.endog_names
+        ynames = self._ynames_map
+        ynames = MultinomialResults._maybe_convert_ynames_int(ynames)
+        # use range below to ensure sortedness
+        ynames = [ynames[key] for key in range(int(self.J))]
+        idx = MultiIndex.from_product((ynames[1:], self.data.xnames),
+                                      names=(yname, None))
+        self.data.cov_names = idx
 
     def pdf(self, eXB):
         """
@@ -4050,7 +4063,8 @@ def __init__(self, model, mlefit):
         self.J = model.J
         self.K = model.K
 
-    def _maybe_convert_ynames_int(self, ynames):
+    @staticmethod
+    def _maybe_convert_ynames_int(ynames):
         # see if they're integers
         issue_warning = False
         msg = ('endog contains values are that not int-like. Uses string '
@@ -4213,75 +4227,108 @@ def __init__(self, model, mlefit):
 
 class OrderedResultsWrapper(lm.RegressionResultsWrapper):
     pass
+
+
 wrap.populate_wrapper(OrderedResultsWrapper, OrderedResults)
 
+
 class CountResultsWrapper(lm.RegressionResultsWrapper):
     pass
+
+
 wrap.populate_wrapper(CountResultsWrapper, CountResults)
 
+
 class NegativeBinomialResultsWrapper(lm.RegressionResultsWrapper):
     pass
+
+
 wrap.populate_wrapper(NegativeBinomialResultsWrapper,
                       NegativeBinomialResults)
 
+
 class GeneralizedPoissonResultsWrapper(lm.RegressionResultsWrapper):
     pass
+
+
 wrap.populate_wrapper(GeneralizedPoissonResultsWrapper,
                       GeneralizedPoissonResults)
 
+
 class PoissonResultsWrapper(lm.RegressionResultsWrapper):
     pass
-    #_methods = {
-    #        "predict_prob" : "rows",
-    #        }
-    #_wrap_methods = lm.wrap.union_dicts(
-    #                            lm.RegressionResultsWrapper._wrap_methods,
-    #                            _methods)
+
+
 wrap.populate_wrapper(PoissonResultsWrapper, PoissonResults)
 
+
 class L1CountResultsWrapper(lm.RegressionResultsWrapper):
     pass
 
+
 class L1PoissonResultsWrapper(lm.RegressionResultsWrapper):
     pass
-    #_methods = {
+    # _methods = {
     #        "predict_prob" : "rows",
     #        }
-    #_wrap_methods = lm.wrap.union_dicts(
+    # _wrap_methods = lm.wrap.union_dicts(
     #                            lm.RegressionResultsWrapper._wrap_methods,
     #                            _methods)
+
+
 wrap.populate_wrapper(L1PoissonResultsWrapper, L1PoissonResults)
 
+
 class L1NegativeBinomialResultsWrapper(lm.RegressionResultsWrapper):
     pass
+
+
 wrap.populate_wrapper(L1NegativeBinomialResultsWrapper,
                       L1NegativeBinomialResults)
 
+
 class L1GeneralizedPoissonResultsWrapper(lm.RegressionResultsWrapper):
     pass
+
+
 wrap.populate_wrapper(L1GeneralizedPoissonResultsWrapper,
                       L1GeneralizedPoissonResults)
 
+
 class BinaryResultsWrapper(lm.RegressionResultsWrapper):
-    _attrs = {"resid_dev" : "rows",
-              "resid_generalized" : "rows",
-              "resid_pearson" : "rows",
-              "resid_response" : "rows"
+    _attrs = {"resid_dev": "rows",
+              "resid_generalized": "rows",
+              "resid_pearson": "rows",
+              "resid_response": "rows"
               }
     _wrap_attrs = wrap.union_dicts(lm.RegressionResultsWrapper._wrap_attrs,
                                    _attrs)
+
+
 wrap.populate_wrapper(BinaryResultsWrapper, BinaryResults)
 
+
 class L1BinaryResultsWrapper(lm.RegressionResultsWrapper):
     pass
+
+
 wrap.populate_wrapper(L1BinaryResultsWrapper, L1BinaryResults)
 
+
 class MultinomialResultsWrapper(lm.RegressionResultsWrapper):
-    _attrs = {"resid_misclassified" : "rows"}
+    _attrs = {"resid_misclassified": "rows"}
     _wrap_attrs = wrap.union_dicts(lm.RegressionResultsWrapper._wrap_attrs,
-            _attrs)
+                                   _attrs)
+    _methods = {'conf_int': 'multivariate_confint'}
+    _wrap_methods = wrap.union_dicts(lm.RegressionResultsWrapper._wrap_methods,
+                                     _methods)
+
+
 wrap.populate_wrapper(MultinomialResultsWrapper, MultinomialResults)
 
+
 class L1MultinomialResultsWrapper(lm.RegressionResultsWrapper):
     pass
+
+
 wrap.populate_wrapper(L1MultinomialResultsWrapper, L1MultinomialResults)

statsmodels/discrete/tests/test_discrete.py

@@ -13,16 +13,19 @@
 import warnings
 
 import numpy as np
-import pandas as pd
 from numpy.testing import (assert_, assert_raises, assert_almost_equal,
                            assert_equal, assert_array_equal, assert_allclose,
                            assert_array_less)
+import pandas as pd
+from pandas.testing import assert_index_equal
 import pytest
+from scipy import stats
 
 from statsmodels.discrete.discrete_model import (Logit, Probit, MNLogit,
-                                                Poisson, NegativeBinomial,
-                                                CountModel, GeneralizedPoisson,
-                                                NegativeBinomialP)
+                                                 Poisson, NegativeBinomial,
+                                                 CountModel,
+                                                 GeneralizedPoisson,
+                                                 NegativeBinomialP)
 from statsmodels.discrete.discrete_margins import _iscount, _isdummy
 import statsmodels.api as sm
 import statsmodels.formula.api as smf
@@ -2356,8 +2359,22 @@ def test_unchanging_degrees_of_freedom():
 
 def test_mnlogit_float_name():
     df = pd.DataFrame({"A": [0., 1.1, 0, 0, 1.1], "B": [0, 1, 0, 1, 1]})
-    result = smf.mnlogit(formula="A ~ B", data=df).fit()
     with pytest.warns(SpecificationWarning,
                       match='endog contains values are that not int-like'):
-        summ = result.summary().as_text()
+        result = smf.mnlogit(formula="A ~ B", data=df).fit()
+    summ = result.summary().as_text()
     assert 'A=1.1' in summ
+
+
+def test_cov_confint_pandas():
+    data = sm.datasets.anes96.load(as_pandas=True)
+    exog = sm.add_constant(data.exog, prepend=False)
+    res1 = sm.MNLogit(data.endog, exog).fit(method="newton", disp=0)
+    cov = res1.cov_params()
+    ci = res1.conf_int()
+    se = np.sqrt(np.diag(cov))
+    se2 = (ci.iloc[:, 1] - ci.iloc[:, 0]) / (2 * stats.norm.ppf(0.975))
+    assert_allclose(se, se2)
+    assert_index_equal(ci.index, cov.index)
+    assert_index_equal(cov.index, cov.columns)
+    assert isinstance(ci.index, pd.MultiIndex)

statsmodels/tsa/vector_ar/tests/test_var.py

@@ -11,6 +11,8 @@
 import sys
 
 import numpy as np
+import pandas as pd
+from pandas.testing import assert_index_equal
 import pytest
 
 
@@ -803,3 +805,15 @@ def test_exog(self):
 def test_deprecated_attributes_varresults(bivariate_var_result, attr):
     with pytest.warns(FutureWarning):
         getattr(bivariate_var_result, attr)
+
+
+def test_var_cov_params(bivariate_var_data):
+    df = pd.DataFrame(bivariate_var_data, columns=['x', 'y'])
+    mod = VAR(df)
+    res = mod.fit(2)
+    cov = res.cov_params()
+    assert isinstance(cov, pd.DataFrame)
+    exog_names = ('const', 'L1.x', 'L1.y', 'L2.x', 'L2.y')
+    index = pd.MultiIndex.from_product((exog_names, ('x', 'y')))
+    assert_index_equal(cov.index, cov.columns)
+    assert_index_equal(cov.index, index)