# Biased (?) results from Poisson regression model

**URL:** <https://discourse.pymc.io/t/biased-results-from-poisson-regression-model/4349>\
**Category:** Questions\
**Created:** [January 13, 2020, 8:10pm UTC](https://discourse.pymc.io/t/biased-results-from-poisson-regression-model/4349 "2020-01-13T20:10:47Z")\
**Posts on this page:** 1\
**Showing post:** 5

<div class="post-metadata">

**Author:** ![birnbera](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/birnbera/32/2262_2.png) [@birnbera](https://discourse.pymc.io/u/birnbera)\
**Post date:** [January 17, 2020, 8:19pm UTC](https://discourse.pymc.io/t/biased-results-from-poisson-regression-model/4349/5 "2020-01-17T20:19:12Z")

</div>

Hi @AlexAndorra , Sorry for the slow reply on this. Thank you for the link to your notebook. I did some prior predictive sampling and, as you mentioned, it revealed that insane values were being sampled with sigmas even at a modest size (e.g. 10). After using more regularizing priors (e.g. `sigma=1`) I found that switching the model error term to use an `InverseGamma` prior with `alpha=5` and `beta=5`, worked much better from a computational perspective. It sped the sampling way up and eliminated even the occasional divergence during sampling.

Even though the inverse gamma distribution is the conjugate prior for the scale parameter of a normal distribution, I wouldn’t have expected that to be particularly important for MCMC sampling purposes.

Finally, instead of using a single error term for the all of the `log_rate`'s, I replicated the prior for each observation. Originally I had thought that the scalar error term would be broadcasted up to the shape of `mu`, but it looks like it fits a single error to all of the `log_rate`'s. This makes sense now that I realized it, but I’m still not sure I did the right thing here… I’m worried that I’m confusing the error associated with the `log_rate` parameter and the residual of the fit, but when the posterior predictive check looks a lot better when fitting the `log_rate` for each observation using its own `sigma`. Any insight on this would be helpful.

The updated model:

```auto
with pm.Model() as model:
    model_intercept = pm.Normal('model_intercept', mu=0, sigma=1)
    factor_coefs = pm.Normal('factor_coefs', mu=0, sigma=1, shape=len(dummies.columns))
    covariate_coefs = pm.Normal('covariate_coefs', mu=0, sigma=1, shape=len(covariates.columns))

    mu = pm.Deterministic(
        'mu',
        model_intercept + 
        pm.math.dot(dummies, factor_coefs) + 
        pm.math.dot(covariates, covariate_coefs) +
        np.log(exposure)
    )

    err = pm.InverseGamma('err', alpha=5, beta=5, shape=n_obs)
    log_rate = pm.Normal('log_rate', mu=mu, sigma=err, shape=n_obs))
    numex = pm.Poisson('numex', mu=pm.math.exp(log_rate), observed=numex_obs)

```

---

_[View the full topic](https://discourse.pymc.io/t/biased-results-from-poisson-regression-model/4349)._