Repository: statsmodels/statsmodels (Python)

Issue:
Zero-Inflated Poisson/Negative Binomial: not predicting for arrays of different size?
**edit** same bug as #4757

Hi,

I am modelling a discrete skewed distribution with several of the available methods in `statsmodels`, namely:  `Poisson`, `Negative Binomial`, `Zero Inflated Poisson`, and `Zero Inflated Negative Binomial`. What I am doing is training with 60% of my samples, and testing the model with the remaining 40%. For this 40%, I obtain the RMSE to have a performance metric for each method.

There are no problems with the Poisson and the Negative Binomial: I use the `.fit()` method on my training samples, and the `.predict()` method on the test samples. However, I observed that for the Zero-Inflated models, **I cannot predict for the testing set**, as you can see in this snippet:

 ```
    print("Model: Zero Inflated Poisson")
    zip_mod = sm.ZeroInflatedPoisson(ytrain, xtrain).fit(method="newton", maxiter=50)
    **zip_mean_pred = zip_mod.predict(xtest)**
    zip_ppf_obs = stats.poisson.ppf(q=0.95, mu=zip_mean_pred)
    zip_rmse = np.sqrt(mean_squared_error(ytest, zip_ppf_obs))
    print("RMSE Zero-Inflated Poisson", zip_rmse)

    print("Model: Zero Inflated Neg. Binomial")
    zinb_mod = sm.ZeroInflatedNegativeBinomialP(ytrain, xtrain).fit(method="newton", maxiter=50)
    zinb_pred = zinb_mod.predict(xtest)
    zinb_rmse = np.sqrt(mean_squared_error(ytest, zinb_pred))
    print("RMSE Zero-Inflated Negative Binomial: ", zinb_rmse)
```

The error that I get says: 
```
Traceback (most recent call last):
  File "/home/irene/.../four_models.py", line 569, in <module>
    test_poi_nb_zip_zinb_tiny_subset(meta_m, m)
  File "/home/irene/.../four_models.py", line 325, in test_poi_nb_zip_zinb_tiny_subset
    **zinb_pred = zinb_mod.predict(xtest)**
  File "/home/irene/.conda/envs/defaultpy34/lib/python2.7/site-packages/statsmodels/base/model.py", line 854, in predict
    **kwargs)
  File "/home/irene/.conda/envs/defaultpy34/lib/python2.7/site-packages/statsmodels/discrete/count_model.py", line 444, in predict
    prob_zero = (1 - prob_main) + prob_main * np.exp(llf)
**ValueError: operands could not be broadcast together with shapes (480,) (320,)** 
```
Perhaps I am doing something wrong, but I do not understand why these two methods require a test sample for the predict method as big as the one used in the fitting. Is there any way to pass to these methods a test array with a different size than the training? I tried training and testing with the same dataset and it works.

Thanks a lot for all the work! :-) 


Hints:
I forgot to mention that the same procedure that you see here for the zero inflated models is applied above it to the normal poisson and negative binomial with not problem. Thanks.
Hi, there is a good chance that this is a bug. Given that the models are new, there might be cases that are not included in the test suite and that could be buggy.

Can you make a reproducible example? 
It can have random number, just so I can reproduce the shape mismatch.

(temporarily prio-high, at least until I know what's going on.)
Based on a quick look at the source.
If exog_infl is not specified in predict, then the one of the model is used which will conflict with the predict exog.

Try to specify `exog_infl`. 
In the case of constant inflation probability as in your case, you need `exog_infl=np.ones(len(exog_predict))`

the same is true for `exposure`, which also uses the model exposure if it is available but no predict exposure is specified.
Thanks for the quick response! I made the following dummy code causing the error I was mentioning, I hope it is useful to trace the bug in the code.
```
import numpy as np
import statsmodels.api as sm
from scipy import stats
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split

X = np.random.randint(99, size=(800, 21))
Y = np.random.randint(2, size=(800, 1))

xtrain, xtest, ytrain, ytest = train_test_split(X, Y, train_size=0.60, random_state=42)

print(xtrain.shape, ytrain.shape, xtest.shape, ytest.shape)

print("Model: Zero Inflated Poisson")
zip_mod = sm.ZeroInflatedPoisson(ytrain, xtrain).fit(method="newton", maxiter=50)
zip_mean_pred = zip_mod.predict(xtest)
zip_ppf_obs = stats.poisson.ppf(q=0.95, mu=zip_mean_pred)
zip_rmse = np.sqrt(mean_squared_error(ytest, zip_ppf_obs))

print("Model: Zero Inflated Neg. Binomial")
zinb_mod = sm.ZeroInflatedNegativeBinomialP(ytrain, xtrain).fit(method="newton", maxiter=50)
zinb_pred = zinb_mod.predict(xtest)
zinb_rmse = np.sqrt(mean_squared_error(ytest, zinb_pred))

print("RMSE ZIP", zip_rmse)
print("RMSE ZINB: ", zinb_rmse)
```


If I add to the `.predict()` methods the `exog_infl` as follows:
```
zip_mean_pred = zip_mod.predict(xtest, exog_infl=np.ones(len(xtest)))
zinb_pred = zinb_mod.predict(xtest, exog_infl=np.ones(len(xtest)))
```

I get the following error:
```
Traceback (most recent call last):
  File "/home/irene/.../test_zip_zinb.py", line 29, in <module>
    zip_mean_pred = zip_mod.predict(xtest, exog_infl=np.ones(len(xtest)))
  File "/home/irene/.conda/envs/defaultpy34/lib/python2.7/site-packages/statsmodels/base/model.py", line 854, in predict
    **kwargs)
  File "/home/irene/.conda/envs/defaultpy34/lib/python2.7/site-packages/statsmodels/discrete/count_model.py", line 415, in predict
    prob_main = 1 - self.model_infl.predict(params_infl, exog_infl)
  File "/home/irene/.conda/envs/defaultpy34/lib/python2.7/site-packages/statsmodels/discrete/discrete_model.py", line 447, in predict
    return self.cdf(np.dot(exog, params))
ValueError: shapes (320,) and (1,) not aligned: 320 (dim 0) != 1 (dim 0)
```
The `.shape` of the `np.ones(len(xtest))` is `(320,)` just in case it is important.

I just tried it. exog_infl in predict needs to be 2dim

`zip_mean_pred = zip_mod.predict(xtest, exog_infl=np.ones((len(xtest), 1)))`
this works, but the example has convergence problems or hessian inversion problem if I use method='nm' as optimizer. (The data doesn't generate anything close to ZIP or ZINB)

general aside: 
Estimation of a NegativBinomial(P) model and ZI variant **needs** overdispersion relative to Poisson in the data. If there is no overdispersion, then the loglike and derivative computation end with a zero division. Poisson is the limit of NegativeBinomial as overdispersion disappears, but computationally we need to stay away from that boundary case.
GeneralizedPoisson and it's ZI version allow for overdispersion, no excess dispersion and a small amount of underdispersion, i.e. the Poisson special case is in the interior of the GP parameter space.
(In the zi version excess dispersion refers to the non-zero inflated part.)

I tried your fix and it works! 
```
print("Model: Zero Inflated Poisson")
zip_mod = sm.ZeroInflatedPoisson(ytrain, xtrain).fit(method="nm", maxiter=50)
zip_mean_pred = zip_mod.predict(xtest, exog_infl=np.ones((len(xtest), 1)))
zip_ppf_obs = stats.poisson.ppf(q=0.95, mu=zip_mean_pred)
zip_rmse = np.sqrt(mean_squared_error(ytest, zip_ppf_obs))

print("Model: Zero Inflated Neg. Binomial")
zinb_mod = sm.ZeroInflatedNegativeBinomialP(ytrain, xtrain).fit(method="nm", maxiter=50)
zinb_pred = zinb_mod.predict(xtest, exog_infl=np.ones((len(xtest), 1)))
zinb_rmse = np.sqrt(mean_squared_error(ytest, zinb_pred))
```
It gives some warnings as you said, like this:
```
HessianInversionWarning: Inverting hessian failed, no bse or cov_params available
  'available', HessianInversionWarning)
```
But it processes my testing dataset. Tried the "nm" and "newton" optimizers.
So now i can compute the RMSE:
```
RMSE ZIP: 1.5642490210960658
RMSE ZINB:  0.50817991749367819
```
Thanks for the super fast reaction, much appreciated!

PS - Regarding this "overdispersion", is there any paper/manual you recommend to learn more about it? :-)

> PS - Regarding this "overdispersion", is there any paper/manual you recommend to learn more about it?

If you do a google search for "overdispersion" or  "overdispersion poisson", then you will find some introductory explanation,
e.g. I never read this
http://biometry.github.io/APES//LectureNotes/2016-JAGS/Overdispersion/OverdispersionJAGS.html

Zero-inflation can cause overdispersion if we just estimate a poisson model. If we estimate a zero-inflated model, then the zero inflation is handled separately but there can still be a (over-, under-) dispersion issue for the not inflated part.

I don't have directly a reference because I haven't looked at it recently, and most of my references are more technical for specific issues (e.g. what we don't have yet is testing zero-inflated Poisson versus zero-inflated NegativeBinomial)

Here are links to two of my notebooks which are a mix of something I started to write for "public" consumption, and then often ended up investigating some specific problems. Those were written when I reviewed and finished up the count model PR and I'm not sure those notebooks are up to date with current statsmodels.
https://gist.github.com/josef-pkt/c932904296270d75366a24ee92a4eb2f
https://gist.github.com/josef-pkt/b528c099592bdb0239e6e750d365f7fc

Patch:
diff --git a/statsmodels/discrete/count_model.py b/statsmodels/discrete/count_model.py
--- a/statsmodels/discrete/count_model.py
+++ b/statsmodels/discrete/count_model.py
@@ -60,11 +60,13 @@ def __init__(self, endog, exog, exog_infl=None, offset=None,
 
         if exog_infl is None:
             self.k_inflate = 1
+            self._no_exog_infl = True
             self.exog_infl = np.ones((endog.size, self.k_inflate),
                                      dtype=np.float64)
         else:
             self.exog_infl = exog_infl
             self.k_inflate = exog_infl.shape[1]
+            self._no_exog_infl = False
 
         if len(exog.shape) == 1:
             self.k_exog = 1
@@ -387,19 +389,35 @@ def predict(self, params, exog=None, exog_infl=None, exposure=None,
         Notes
         -----
         """
+        no_exog = False
         if exog is None:
+            no_exog = True
             exog = self.exog
 
         if exog_infl is None:
-            exog_infl = self.exog_infl
+            if no_exog:
+                exog_infl = self.exog_infl
+            else:
+                if self._no_exog_infl:
+                    exog_infl = np.ones((len(exog), 1))
+        else:
+            exog_infl = np.asarray(exog_infl)
+            if exog_infl.ndim == 1 and self.k_inflate == 1:
+                exog_infl = exog_infl[:, None]
 
         if exposure is None:
-            exposure = getattr(self, 'exposure', 0)
+            if no_exog:
+                exposure = getattr(self, 'exposure', 0)
+            else:
+                exposure = 0
         else:
             exposure = np.log(exposure)
 
         if offset is None:
-            offset = 0
+            if no_exog:
+                offset = getattr(self, 'offset', 0)
+            else:
+                offset = 0
 
         params_infl = params[:self.k_inflate]
         params_main = params[self.k_inflate:]


Test patch:
diff --git a/statsmodels/discrete/tests/test_count_model.py b/statsmodels/discrete/tests/test_count_model.py
--- a/statsmodels/discrete/tests/test_count_model.py
+++ b/statsmodels/discrete/tests/test_count_model.py
@@ -144,19 +144,18 @@ def test_exposure(self):
 
         assert_allclose(res3.params, self.res1.params, atol=1e-6, rtol=1e-6)
         fitted1 = self.res1.predict()
-        fitted3 = self.res1.predict()
+        fitted3 = res3.predict()
         assert_allclose(fitted3, fitted1, atol=1e-6, rtol=1e-6)
 
         ex = model1.exog
         ex_infl = model1.exog_infl
         offset = model1.offset
         fitted1_0 = self.res1.predict(exog=ex, exog_infl=ex_infl,
-                                      offset=offset)
+                                      offset=offset.tolist())
         fitted3_0 = res3.predict(exog=ex, exog_infl=ex_infl,
                                  exposure=np.exp(offset))
         assert_allclose(fitted3_0, fitted1_0, atol=1e-6, rtol=1e-6)
 
-
         ex = model1.exog[:10:2]
         ex_infl = model1.exog_infl[:10:2]
         offset = offset[:10:2]
@@ -175,6 +174,11 @@ def test_exposure(self):
         assert_allclose(fitted1_2, fitted1[:10:2], atol=1e-6, rtol=1e-6)
         assert_allclose(fitted3_2, fitted1[:10:2], atol=1e-6, rtol=1e-6)
 
+        # without specifying offset and exposure
+        fitted1_3 = self.res1.predict(exog=ex, exog_infl=ex_infl)
+        fitted3_3 = res3.predict(exog=ex, exog_infl=ex_infl)
+        assert_allclose(fitted3_3, fitted1_3, atol=1e-6, rtol=1e-6)
+
 
 class TestZeroInflatedModelPandas(CheckGeneric):
     @classmethod
@@ -276,6 +280,18 @@ def test_predict_prob(self):
                                              res.predict(), 0.05).T
         assert_allclose(pr, pr2, rtol=0.05, atol=0.05)
 
+    def test_predict_options(self):
+        # check default exog_infl, see #4757
+        res = self.res
+        n = 5
+        pr1 = res.predict(which='prob')
+        pr0 = res.predict(exog=res.model.exog[:n], which='prob')
+        assert_allclose(pr0, pr1[:n], rtol=1e-10)
+
+        fitted1 = res.predict()
+        fitted0 = res.predict(exog=res.model.exog[:n])
+        assert_allclose(fitted0, fitted1[:n], rtol=1e-10)
+
 
 @pytest.mark.slow
 class TestZeroInflatedGeneralizedPoisson(CheckGeneric):