Repository: statsmodels/statsmodels (Python)

Issue:
missing endog values in OrderedModel
#### Describe the bug


OrderedModel produces one threshold per occurrence of nan in endog.

If not a pandas Categorical, Labels in OrderedModel are computed using `np.unique` [here](https://github.com/statsmodels/statsmodels/blob/a91eea01799b2fa1bcafe0dea0f00076091d9f5e/statsmodels/miscmodels/ordinal_model.py#L130). However `np.unique` does not squash nans:

```python
In [2]: np.unique([np.nan, np.nan, np.nan, 1, 2, 1, 2])
Out[2]: array([1, 2, nan, nan, nan])
```

I'm not sure if this is a bug in `numpy.unique()` or the intended behavior there, but the result is, instead of either 1 or 0 output values for nan (I'm not sure which would make the most sense, probably none?), a label and threshold is produced for *each* nan:

<img width="387" alt="Screen Shot 2021-02-27 at 14 37 20" src="https://user-images.githubusercontent.com/151929/109388759-56088e80-7909-11eb-953c-5b28a6c70aa7.png">

Since the same input cast to a pandas categorical produces a ValueError:

> ValueError: missing values in categorical endog are not supported

Presumably the same error should be raised in this case as well.

#### Code Sample, a copy-pastable example if possible


```python
import numpy as np
import pandas as pd
import scipy.stats as stats

from statsmodels.miscmodels.ordinal_model import OrderedModel

N = 100
df = pd.DataFrame(
    {
        "target": np.random.randint(1, 3, N).astype(float),
        "factor1": np.random.random(N),
        "factor2": np.random.random(N),
    }
)
df.loc[1:4, "target"] = float("nan")

mod_prob = OrderedModel(df["target"], df[["factor1", "factor2"]], distr="probit")
res_prob = mod_prob.fit(method="bfgs")
res_prob.summary()
```

which produces the output above, with one output per missing value.

Casting the same column to a categorical results in a sensible error:

```python
df["target"] = df["target"].astype(pd.CategoricalDtype(categories=[1, 2, 3], ordered=True))
...
ValueError: missing values in categorical endog are not supported
```


Since this was a pandas dataframe, calling `dropna(subset=["target"])` prior to passing to the model worked just fine, and produced the desired result.


<details>

**Note**: As you can see, there are many issues on our GitHub tracker, so it is very possible that your issue has been posted before. Please check first before submitting so that we do not have to handle and close duplicates.

**Note**: Please be sure you are using the latest released version of `statsmodels`, or a recent build of `master`. If your problem has been fixed in an unreleased version, you might be able to use `master` until a new release occurs. 

**Note**: If you are using a released version, have you verified that the bug exists in the master branch of this repository? It helps the limited resources if we know problems exist in the current master so that they do not need to check whether the code sample produces a bug in the next release.

</details>


If the issue has not been resolved, please file it in the issue tracker.

#### Expected Output

Either of these would make sense to me:

- early error that nans won't work (simplest and seems to match an existing error check in categorical endog handling)
- handle and drop missing values using `np.isnan` (unclear if this would break things to not have a label for an output)


#### Output of ``import statsmodels.api as sm; sm.show_versions()``

<details>

[paste the output of ``import statsmodels.api as sm; sm.show_versions()`` here below this line]

```

INSTALLED VERSIONS
------------------
Python: 3.8.6.final.0
OS: Darwin 20.3.0 Darwin Kernel Version 20.3.0: Thu Jan 21 00:07:06 PST 2021; root:xnu-7195.81.3~1/RELEASE_X86_64 x86_64
byteorder: little
LC_ALL: None
LANG: en_US.UTF-8

statsmodels
===========

Installed: v0.13.0.dev0+220.ga91eea017 (src/statsmodels/statsmodels)

Required Dependencies
=====================

cython: 3.0a6 (/Users/minrk/conda/lib/python3.8/site-packages/Cython)
numpy: 1.19.1 (/Users/minrk/conda/lib/python3.8/site-packages/numpy)
scipy: 1.5.2 (/Users/minrk/conda/lib/python3.8/site-packages/scipy)
pandas: 1.1.2 (/Users/minrk/conda/lib/python3.8/site-packages/pandas)
    dateutil: 2.8.1 (/Users/minrk/conda/lib/python3.8/site-packages/dateutil)
patsy: 0.5.1 (/Users/minrk/conda/lib/python3.8/site-packages/patsy)

Optional Dependencies
=====================

matplotlib: 3.3.2 (/Users/minrk/conda/lib/python3.8/site-packages/matplotlib)
    backend: MacOSX
cvxopt: Not installed
joblib: Not installed

Developer Tools
================

IPython: 8.0.0.dev (/Users/minrk/dev/ip/ipython/IPython)
    jinja2: 2.11.2 (/Users/minrk/conda/lib/python3.8/site-packages/jinja2)
sphinx: 3.2.1 (/Users/minrk/conda/lib/python3.8/site-packages/sphinx)
    pygments: 2.7.1 (/Users/minrk/conda/lib/python3.8/site-packages/pygments)
pytest: 6.1.1 (/Users/minrk/conda/lib/python3.8/site-packages/pytest)
virtualenv: 20.0.33 (/Users/minrk/conda/lib/python3.8/site-packages/virtualenv)
```

</details>



Hints:

Patch:
diff --git a/statsmodels/miscmodels/ordinal_model.py b/statsmodels/miscmodels/ordinal_model.py
--- a/statsmodels/miscmodels/ordinal_model.py
+++ b/statsmodels/miscmodels/ordinal_model.py
@@ -6,6 +6,8 @@
 License: BSD-3
 """
 
+import warnings
+
 from statsmodels.compat.pandas import Appender
 
 import numpy as np
@@ -135,6 +137,10 @@ def __init__(self, endog, exog, offset=None, distr='probit', **kwds):
                 unique, index = np.unique(self.endog, return_inverse=True)
                 self.endog = index
                 labels = unique
+                if np.isnan(labels).any():
+                    msg = ("NaN in dependent variable detected. "
+                           "Missing values need to be removed.")
+                    raise ValueError(msg)
             elif self.endog.ndim == 2:
                 if not hasattr(self, "design_info"):
                     raise ValueError("2-dim endog not supported")
@@ -146,7 +152,7 @@ def __init__(self, endog, exog, offset=None, distr='probit', **kwds):
                 # self.endog = self.endog.argmax(1)
 
         if self.k_constant > 0:
-            raise ValueError("there should not be a constant in the model")
+            raise ValueError("There should not be a constant in the model")
 
         self._initialize_labels(labels, k_levels=k_levels)
 
@@ -187,7 +193,6 @@ def _check_inputs(self, endog, exog):
         """
 
         if not isinstance(self.distr, stats.rv_continuous):
-            import warnings
             msg = (
                 f"{self.distr.name} is not a scipy.stats distribution."
             )


Test patch:
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
@@ -1176,16 +1176,16 @@ def check_score_hessian(results):
     # avoid checking score at MLE, score close to zero
     sc = results.model.score(params * 0.98, scale=1)
     # cs currently (0.9) does not work for all families
-    llfunc = lambda x: results.model.loglike(x, scale=1)
+    llfunc = lambda x: results.model.loglike(x, scale=1)  # noqa
     sc2 = approx_fprime(params * 0.98, llfunc)
-    assert_allclose(sc, sc2, rtol=0.05)
+    assert_allclose(sc, sc2, rtol=1e-4, atol=1e-4)
 
     hess = results.model.hessian(params, scale=1)
     hess2 = approx_hess(params, llfunc)
-    assert_allclose(hess, hess2, rtol=0.05)
-    scfunc = lambda x: results.model.score(x, scale=1)
+    assert_allclose(hess, hess2, rtol=1e-4)
+    scfunc = lambda x: results.model.score(x, scale=1)  # noqa
     hess3 = approx_fprime(params, scfunc)
-    assert_allclose(hess, hess3, rtol=0.05)
+    assert_allclose(hess, hess3, rtol=1e-4)
 
 
 def test_gradient_irls():
diff --git a/statsmodels/miscmodels/tests/test_ordinal_model.py b/statsmodels/miscmodels/tests/test_ordinal_model.py
--- a/statsmodels/miscmodels/tests/test_ordinal_model.py
+++ b/statsmodels/miscmodels/tests/test_ordinal_model.py
@@ -283,7 +283,8 @@ def test_formula_categorical(self):
         assert hasattr(modf2.data, "frame")
         assert not hasattr(modf2, "frame")
 
-        with pytest.raises(ValueError):
+        msg = "Only ordered pandas Categorical"
+        with pytest.raises(ValueError, match=msg):
             # only ordered categorical or numerical endog are allowed
             # string endog raises ValueError
             OrderedModel.from_formula(
@@ -375,7 +376,7 @@ def test_setup(self):
         # test over parameterized model with implicit constant
         formula = "apply ~ 0 + pared + public + gpa + C(dummy)"
 
-        with pytest.raises(ValueError):
+        with pytest.raises(ValueError, match="not be a constant"):
             OrderedModel.from_formula(formula, data, distr='logit')
 
         # ignore constant, so we get results without exception
@@ -482,3 +483,27 @@ def test_attributes(self):
         for attr in attributes:
             assert_allclose(getattr(resp, attr), getattr(res_logit, attr),
                             rtol=1e-4)
+
+
+def test_nan_endog_exceptions():
+    nobs = 15
+    y = np.repeat(np.arange(3), nobs // 3)
+    x = np.column_stack((np.ones(nobs), np.arange(nobs)))
+    with pytest.raises(ValueError, match="not be a constant"):
+        OrderedModel(y, x, distr='logit')
+
+    y_nan = y.astype(float)
+    y_nan[0] = np.nan
+    with pytest.raises(ValueError, match="NaN in dependent variable"):
+        OrderedModel(y_nan, x[:, 1:], distr='logit')
+
+    if hasattr(pd, "CategoricalDtype"):
+        df = pd.DataFrame({
+            "endog": pd.Series(
+                y, dtype=pd.CategoricalDtype([1, 2, 3], ordered=True)),
+            "exog": x[:, 1]
+            })
+
+        msg = "missing values in categorical endog"
+        with pytest.raises(ValueError, match=msg):
+            OrderedModel(df["endog"], df[["exog"]])