Merge 5bdaad4 into 2322d35

statsmodels/regression/mixed_linear_model.py

@@ -156,6 +156,8 @@
 from statsmodels.tools import data as data_tools
 from statsmodels.tools.sm_exceptions import ConvergenceWarning
 
+_warn_cov_sing = "The random effects covariance matrix is singular."
+
 
 def _dot(x, y):
     """
@@ -413,7 +415,10 @@ def get_packed(self, use_sqrt, has_fe=False):
 
         if self.k_re > 0:
             if use_sqrt:
-                L = np.linalg.cholesky(self.cov_re)
+                try:
+                    L = np.linalg.cholesky(self.cov_re)
+                except np.linalg.LinAlgError:
+                    L = np.diag(np.sqrt(np.diag(self.cov_re)))
                 cpa = L[self._ix]
             else:
                 cpa = self.cov_re[self._ix]
@@ -1232,7 +1237,10 @@ def fit_regularized(self, start_params=None, method='l1', alpha=0,
 
         # Get the Hessian including only the nonzero fixed effects,
         # then blow back up to the full size after inverting.
-        hess = self.hessian(params_prof)
+        hess, sing = self.hessian(params_prof)
+        if sing:
+            warnings.warn(_warn_cov_sing)
+
         pcov = np.nan * np.ones_like(hess)
         ii = np.abs(params_prof) > ceps
         ii[self.k_fe:] = True
@@ -1258,27 +1266,49 @@ def fit_regularized(self, start_params=None, method='l1', alpha=0,
 
         return MixedLMResultsWrapper(results)
 
-    def get_fe_params(self, cov_re, vcomp):
+    def get_fe_params(self, cov_re, vcomp, tol=1e-10):
         """
         Use GLS to update the fixed effects parameter estimates.
 
         Parameters
         ----------
-        cov_re : array_like
+        cov_re : array_like (2d)
             The covariance matrix of the random effects.
+        vcomp : array_like (1d)
+            The variance components.
+        tol : float
+            A tolerance parameter to determine when covariances
+            are singular.
 
         Returns
         -------
-        The GLS estimates of the fixed effects parameters.
+        params : ndarray
+            The GLS estimates of the fixed effects parameters.
+        singular : bool
+            True if the covariance is singular
         """
 
         if self.k_fe == 0:
-            return np.array([])
+            return np.array([]), False
+
+        sing = False
 
         if self.k_re == 0:
             cov_re_inv = np.empty((0, 0))
         else:
-            cov_re_inv = np.linalg.inv(cov_re)
+            w, v = np.linalg.eigh(cov_re)
+            if w.min() < tol:
+                # Singular, use pseudo-inverse
+                sing = True
+                ii = np.flatnonzero(w >= tol)
+                if len(ii) == 0:
+                    cov_re_inv = np.zeros_like(cov_re)
+                else:
+                    vi = v[:, ii]
+                    wi = w[ii]
+                    cov_re_inv = np.dot(vi / wi, vi.T)
+            else:
+                cov_re_inv = np.linalg.inv(cov_re)
 
         # Cache these quantities that do not change.
         if not hasattr(self, "_endex_li"):
@@ -1292,15 +1322,29 @@ def get_fe_params(self, cov_re, vcomp):
         xtxy = 0.
         for group_ix, group in enumerate(self.group_labels):
             vc_var = self._expand_vcomp(vcomp, group_ix)
+            if vc_var.size > 0:
+                if vc_var.min() < tol:
+                    # Pseudo-inverse
+                    sing = True
+                    ii = np.flatnonzero(vc_var >= tol)
+                    vc_vari = np.zeros_like(vc_var)
+                    vc_vari[ii] = 1 / vc_var[ii]
+                else:
+                    vc_vari = 1 / vc_var
+            else:
+                vc_vari = np.empty(0)
             exog = self.exog_li[group_ix]
             ex_r, ex2_r = self._aex_r[group_ix], self._aex_r2[group_ix]
-            solver = _smw_solver(1., ex_r, ex2_r, cov_re_inv, 1 / vc_var)
+            solver = _smw_solver(1., ex_r, ex2_r, cov_re_inv, vc_vari)
             u = solver(self._endex_li[group_ix])
             xtxy += np.dot(exog.T, u)
 
-        fe_params = np.linalg.solve(xtxy[:, 0:-1], xtxy[:, -1])
+        if sing:
+            fe_params = np.dot(np.linalg.pinv(xtxy[:, 0:-1]), xtxy[:, -1])
+        else:
+            fe_params = np.linalg.solve(xtxy[:, 0:-1], xtxy[:, -1])
 
-        return fe_params
+        return fe_params, sing
 
     def _reparam(self):
         """
@@ -1451,15 +1495,18 @@ def loglike(self, params, profile_fe=True):
 
         # Move to the profile set
         if profile_fe:
-            fe_params = self.get_fe_params(cov_re, vcomp)
+            fe_params, sing = self.get_fe_params(cov_re, vcomp)
+            if sing:
+                self._cov_sing += 1
         else:
             fe_params = params.fe_params
 
         if self.k_re > 0:
             try:
                 cov_re_inv = np.linalg.inv(cov_re)
             except np.linalg.LinAlgError:
-                cov_re_inv = None
+                cov_re_inv = np.linalg.pinv(cov_re)
+                self._cov_sing += 1
             _, cov_re_logdet = np.linalg.slogdet(cov_re)
         else:
             cov_re_inv = np.zeros((0, 0))
@@ -1579,7 +1626,12 @@ def score(self, params, profile_fe=True):
                 has_fe=False)
 
         if profile_fe:
-            params.fe_params = self.get_fe_params(params.cov_re, params.vcomp)
+            params.fe_params, sing = \
+                self.get_fe_params(params.cov_re, params.vcomp)
+
+            if sing:
+                msg = "Random effects covariance is singular"
+                warnings.warn(msg)
 
         if self.use_sqrt:
             score_fe, score_re, score_vc = self.score_sqrt(
@@ -1644,7 +1696,8 @@ def score_full(self, params, calc_fe):
         try:
             cov_re_inv = np.linalg.inv(cov_re)
         except np.linalg.LinAlgError:
-            cov_re_inv = None
+            cov_re_inv = np.linalg.pinv(cov_re)
+            self._cov_sing += 1
 
         score_fe = np.zeros(self.k_fe)
         score_re = np.zeros(self.k_re2)
@@ -1818,6 +1871,9 @@ def hessian(self, params):
         -------
         hess : 2d ndarray
             The Hessian matrix, evaluated at `params`.
+        sing : boolean
+            If True, the covariance matrix is singular and a
+            pseudo-inverse is returned.
         """
 
         if type(params) is not MixedLMParams:
@@ -1828,8 +1884,14 @@ def hessian(self, params):
         fe_params = params.fe_params
         vcomp = params.vcomp
         cov_re = params.cov_re
+        sing = False
+
         if self.k_re > 0:
-            cov_re_inv = np.linalg.inv(cov_re)
+            try:
+                cov_re_inv = np.linalg.inv(cov_re)
+            except np.linalg.LinAlgError:
+                cov_re_inv = np.linalg.pinv(cov_re)
+                sing = True
         else:
             cov_re_inv = np.empty((0, 0))
 
@@ -1852,10 +1914,16 @@ def hessian(self, params):
         for group_ix, group in enumerate(self.group_labels):
 
             vc_var = self._expand_vcomp(vcomp, group_ix)
+            vc_vari = np.zeros_like(vc_var)
+            ii = np.flatnonzero(vc_var >= 1e-10)
+            if len(ii) > 0:
+                vc_vari[ii] = 1 / vc_var[ii]
+            if len(ii) < len(vc_var):
+                sing = True
 
             exog = self.exog_li[group_ix]
             ex_r, ex2_r = self._aex_r[group_ix], self._aex_r2[group_ix]
-            solver = _smw_solver(1., ex_r, ex2_r, cov_re_inv, 1 / vc_var)
+            solver = _smw_solver(1., ex_r, ex2_r, cov_re_inv, vc_vari)
 
             # The residuals
             resid = self.endog_li[group_ix]
@@ -1957,7 +2025,7 @@ def hessian(self, params):
         hess[self.k_fe:, 0:self.k_fe] = hess_fere
         hess[self.k_fe:, self.k_fe:] = hess_re
 
-        return hess
+        return hess, sing
 
     def get_scale(self, fe_params, cov_re, vcomp):
         """
@@ -1982,7 +2050,8 @@ def get_scale(self, fe_params, cov_re, vcomp):
         try:
             cov_re_inv = np.linalg.inv(cov_re)
         except np.linalg.LinAlgError:
-            cov_re_inv = None
+            cov_re_inv = np.linalg.pinv(cov_re)
+            warnings.warn(_warn_cov_sing)
 
         qf = 0.
         for group_ix, group in enumerate(self.group_labels):
@@ -2079,7 +2148,7 @@ def fit(self, start_params=None, reml=True, niter_sa=0,
         self.reml = reml
         self.cov_pen = cov_pen
         self.fe_pen = fe_pen
-
+        self._cov_sing = 0
         self._freepat = free
 
         if full_output:
@@ -2155,7 +2224,7 @@ def fit(self, start_params=None, reml=True, niter_sa=0,
             params, self.k_fe, self.k_re, use_sqrt=self.use_sqrt, has_fe=False)
         cov_re_unscaled = params.cov_re
         vcomp_unscaled = params.vcomp
-        fe_params = self.get_fe_params(cov_re_unscaled, vcomp_unscaled)
+        fe_params, sing = self.get_fe_params(cov_re_unscaled, vcomp_unscaled)
         params.fe_params = fe_params
         scale = self.get_scale(fe_params, cov_re_unscaled, vcomp_unscaled)
         cov_re = scale * cov_re_unscaled
@@ -2171,7 +2240,10 @@ def fit(self, start_params=None, reml=True, niter_sa=0,
         # Hessian with respect to the random effects covariance matrix
         # (not its square root).  It is used for obtaining standard
         # errors, not for optimization.
-        hess = self.hessian(params)
+        hess, sing = self.hessian(params)
+        if sing:
+            warnings.warn(_warn_cov_sing)
+
         hess_diag = np.diag(hess)
         if free is not None:
             pcov = np.zeros_like(hess)

statsmodels/regression/tests/test_lme.py

@@ -157,7 +157,10 @@ def test_compare_numdiff(self, use_sqrt, reml, profile_fe):
             # about the Hessian for the square root
             # transformed parameter).
             if (profile_fe is False) and (use_sqrt is False):
-                hess = -model.hessian(rslt.params_object)
+                hess, sing = model.hessian(rslt.params_object)
+                if sing:
+                    pytest.fail("hessian should not be singular")
+                hess *= -1
                 params_vec = rslt.params_object.get_packed(
                     use_sqrt=False, has_fe=True)
                 loglike_h = loglike_function(
@@ -469,7 +472,7 @@ def test_dietox_slopes(self):
         data = pd.read_csv(fname)
         model = MixedLM.from_formula(
             "Weight ~ Time", groups="Pig", re_formula="1 + Time", data=data)
-        result = model.fit(method='powell')
+        result = model.fit(method="cg")
 
         # fixef(r)
         assert_allclose(
@@ -496,7 +499,7 @@ def test_dietox_slopes(self):
         data = pd.read_csv(fname)
         model = MixedLM.from_formula(
             "Weight ~ Time", groups="Pig", re_formula="1 + Time", data=data)
-        result = model.fit(method='powell', reml=False)
+        result = model.fit(method='cg', reml=False)
 
         # fixef(r)
         assert_allclose(result.fe_params, np.r_[15.73863, 6.93902], rtol=1e-5)
@@ -1235,6 +1238,24 @@ def test_smw_logdet(self, p, q, r, s):
         check_smw_logdet(p, q, r, s)
 
 
+def test_singular():
+    #7051
+
+    np.random.seed(3423)
+    n = 100
+
+    data = np.random.randn(n, 2)
+    df = pd.DataFrame(data, columns=['Y', 'X'])
+    df['class'] = pd.Series([i % 3 for i in df.index], index=df.index)
+
+    with pytest.warns(Warning) as wrn:
+        md = MixedLM.from_formula("Y ~ X", df, groups=df['class'])
+        mdf = md.fit()
+        mdf.summary()
+        if not wrn:
+            pytest.fail("warning expected")
+
+
 def test_get_distribution():
 
     np.random.seed(234)