# Normal-Poisson Mixture: ARM 15.1

**URL:** <https://discourse.mc-stan.org/t/normal-poisson-mixture-arm-15-1/32069>\
**Category:** Modeling\
**Tags:** poisson, specification\
**Created:** [July 13, 2023, 9:57am UTC](https://discourse.mc-stan.org/t/normal-poisson-mixture-arm-15-1/32069 "2023-07-13T09:57:06Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![fditraglia](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/fditraglia/32/16656_2.png) [@fditraglia](https://discourse.mc-stan.org/u/fditraglia)\
**Post date:** [July 13, 2023, 9:57am UTC](https://discourse.mc-stan.org/t/normal-poisson-mixture-arm-15-1/32069/1 "2023-07-13T09:57:06Z")

</div>

I’d like to use STAN to fit a model similar to that in section 15.1 of Gelman and Hill’s “Data Analysis Using Regression and Multilevel / Hierarchical Models” (ARM):

\begin{align\*} y\_i &\sim \text{Poisson}(u\_i e^{X\_i\beta + \epsilon\_i})\\ \epsilon\_i &\sim \text{Normal}(0, \sigma\_\epsilon^2) \end{align\*}

where u\_i is an observed offset, y\_i is an observed outcome, X\_i is an observed vector of covariates, and \epsilon\_i is an unobserved error. I know how to fit this model without the \epsilon\_i term, so my question is about how to implement the normal mixture.

I found [this repo](https://github.com/stan-dev/example-models/tree/master/ARM) with STAN code for most of the examples from ARM, but it doesn’t contain code for any of the examples based on the police stops dataset. (Incidentally, a version of that dataset with noise added is available [here](http://www.stat.columbia.edu/~gelman/arm/examples/police/frisk_with_noise.dat) and some discussion is available [here](https://www.clayford.net/statistics/poisson-regression-ch-6-of-gelman-and-hill/)).

Any suggestions you can offer on how to fit the model described above using STAN would be most appreciated. If doing so is hard / impossible / inadvisable I’d also be interested to know why. I’m aware that a negative binomial model provides another way to allow over-dispersion in count data, and I know how to fit such a model in STAN. The reason I’m interested in this normal mixture of Poissons is that I ultimately want to build a more complicated model that uses this one as a building block, including correlated normal error terms across multiple “equations.” (Basically, an instrumental-variables type model but for count data.)

---

<div class="post-metadata">

**Author:** ![jsocolar](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jsocolar/32/2486_2.png) [@jsocolar](https://discourse.mc-stan.org/u/jsocolar)\
**Post date:** [July 13, 2023, 4:36pm UTC](https://discourse.mc-stan.org/t/normal-poisson-mixture-arm-15-1/32069/2 "2023-07-13T16:36:46Z")

</div>

Here we are taking the model without the \epsilon\_i and adding an observation level random effect. This will likely be fit most effectively in the non-centered parameterization, which would look like this:

```stan
parameters{
  real<lower = 0> sigma; // sd of the normal
  vector<multiplier = sigma>[N] errors; // the errors, non-centered parameterization
  ...
}

model{
  errors ~ normal(0, sigma);
  ...<compute the linear predictor on the log scale>...
  <linear predictor> += errors;
  ...<the Poisson likelihood>...
}

```

---

<div class="post-metadata">

**Author:** ![fditraglia](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/fditraglia/32/16656_2.png) [@fditraglia](https://discourse.mc-stan.org/u/fditraglia)\
**Post date:** [August 16, 2023, 2:16pm UTC](https://discourse.mc-stan.org/t/normal-poisson-mixture-arm-15-1/32069/3 "2023-08-16T14:16:58Z")

</div>

This is really helpful: thanks! I think I was getting hung up on the issue of the mixture likelihood lacking a closed-form, whereas it’s more natural to just think of this as a hierarchical model as your example code suggests. I’m still a little unclear on how to correctly use the affinely-transformed parameter feature `<multiplier = sigma>` so I ended up implementing the non-centered parameterization “the hard way.” I coded up a little example with a single regressor and it seems to work well in simulated data:

```plaintext
data {
  int<lower=0> N;
  vector[N] x;
  array[N] int<lower = 0> y;
}

parameters {
  real<lower = 0> sigma;
  real alpha; 
  real beta;
  vector[N] errors;  
}

model {
  sigma ~ exponential(1);
  beta ~ normal(0, 0.2);
  alpha ~ normal(0, 1.5);
  errors ~ normal(0, 1);
  y ~ poisson_log(alpha + beta * x + sigma * errors);
}

```