Biased (?) results from Poisson regression model

Questions
5

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:

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)
1 Like