Merge pull request #93 from ejolly/dev

Fix divide by 0 error for non-full rank square matrices

Author: ejolly

pymer4/models/Lm.py

@@ -289,9 +289,10 @@ def fit(
             # Cluster corrected dof (num clusters - num coef)
             # Differs from stats and statsmodels which do num cluster - 1
             # Ref: http://cameron.econ.ucdavis.edu/research/Cameron_Miller_JHR_2015_February.pdf
-            df = cluster.nunique() - x.shape[1]
+            df = cluster.nunique() - np.linalg.matrix_rank(x)
         else:
-            df = x.shape[0] - x.shape[1]
+            # Use dof calculation that accounts for square matrices to avoid 0 division error or pval lookup error
+            df = x.shape[0] - np.linalg.matrix_rank(x)
             if isinstance(weights, str) and wls_dof_correction:
                 if weight_groups.ngroups != 2:
                     w = "Welch-Satterthwait DOF correction only supported for 2 groups in the data"
@@ -408,6 +409,7 @@ def fit(
         self.rsquared_adj = rsquared_adj(
             self.rsquared, len(res), len(res) - x.shape[1], center_tss
         )
+        # self.rsquared_adj = np.nan
         half_obs = len(res) / 2.0
         ssr = np.dot(res, res.T)
         self.logLike = (-np.log(ssr) * half_obs) - (

pymer4/stats.py

@@ -638,6 +638,9 @@ def rsquared_adj(r, nobs, df_res, has_constant=True):
         float: adjusted coefficient of determination
     """
 
+    # Skip calculation if Dof is not a positive number to avoid 0 division error or calculating negative R^2 values
+    if df_res <= 0:
+        return np.nan
     if has_constant:
         return 1.0 - (nobs - 1) / df_res * (1.0 - r)
     else:

pymer4/utils.py

@@ -235,7 +235,8 @@ def _ols(x, y, robust, n_lags, cluster, all_stats=True, resid_only=False, weight
                 res, X, robust_estimator=robust, n_lags=n_lags, cluster=cluster
             )
         else:
-            sigma = np.sqrt(res.T.dot(res) / (X.shape[0] - X.shape[1]))
+            # Use dof calculation in the denominator that accounts for square matrices to avoid 0 division error
+            sigma = np.sqrt(res.T.dot(res) / (X.shape[0] - np.linalg.matrix_rank(X)))
             se = np.sqrt(np.diag(np.linalg.pinv(np.dot(X.T, X)))) * sigma
 
         t = b / se

pymer4/version.py

@@ -1,4 +1,4 @@
 """Specifies current version of pymer4 to be used by setup.py and __init__.py
 """
 
-__version__ = "0.7.7"
+__version__ = "0.7.8.dev0"