Large Poisson model with Individual effects is too slow

[bgoodri]

@pauper42

An important assumption I’m making is that I need within week variation to identify the “true” price response (let’s take this as given for the purpose of the question). I could do this by including product-week dummy variables but doing this is problematic. To avoid the “incidental parameters problem” associated with these individual effects, and also to reduce the number of parameters to be estimated I use a transformation of the Poisson that leads to an equivalent likelihood function that looks like the product of multinomials (Lancaster 2004). This leads to a model in which there are no dummies but the other parameters of the Poisson are identified through within week variation for each product.

Lancaster reparameterizations are critical when using Stan, although no one does them. However, his 2004 code was written for BUGS and that makes things a bit more confusing when adapting them to Stan. Basically, you want to put the parameters you integrated out back in. In the Poisson case, Lancaster writes something equivalent to

log(alphastar_g) = log(alpha_g)  + log_sum_exp(x[g] * beta);

so in the parameters block you declare

parameters {
  real log_alphastar[groups];
  ...
}

and in the transformed parameters block you map to

transformed parameters {
  real log_alpha[groups];
  for (g in 1:groups) 
    log_alpha[g] = log_alphastar[g] - log_sum_exp(x[g] * beta);
}

Finally, in the model block you use the log_alpha in conjunction with a poisson_log_lpmf likelihood with any prior you want on log_alphastar. This is actually faster than doing the multinomial thing when you call poisson_log_lpmf just once. It also makes things easier to do zero-inflation or overdispersion.