I'm faced with the following statistics problem that I thought a Baysian hierarchical model would give useful results, but I'm not sure how to apply it. A summarised description:

I have daily sales numbers of a range of products (let's call them product A, product B, etc) in a range of shops at different locations (shop 1, shop 2, etc). For each product-shop combination I want to estimate the distribution of sales per day, assumed to be normal (so not trying to take temporal behaviour into account yet, only looking for the average sales rate and uncertainty on it).

How do I model this, making sure that the estimate is informed by both the general behaviour of the specific shop (How well do other products sell in this shop?), but also the general selling rates of the specific product over the range of shops? (If product A is very popular across shops, it seems reasonable to assume this should influence our inference about product A in shop 2).

A simplified variety of this problem, with for example only sales data of product A over a range of shops, seems 'easily' solvable with a 2-stage hierarchical model. But I don't seem to find a guide on how to implement both dimensions of information.

Any suggestions?

  • $\begingroup$ To clarify when you say implement do you mean in terms of the model formulation or how to actually code it up and implement it in software (or both) ? $\endgroup$
    – gowerc
    Mar 18, 2018 at 17:52
  • $\begingroup$ Probably should've clarified that! Model formulation first of all. But I'm also using this problem to get more accustomed to software implementation of bayesian statistical modelling, so happy to get some input there too, preferably python suggestions, but R is fine too. :) $\endgroup$
    – Toekan
    Mar 18, 2018 at 17:56

1 Answer 1


I'm not a Bayesian or hierarchical expert (only done an introductory course) but I would probably model the problem like this

$$ X_{ij} \sim N( \mu_{ij} , \sigma ) \\ \mu_{ij} = \mu_0 + \alpha_i + \beta_j \\ \alpha_i \sim N(0 , \sigma_a ) \\ \beta_i \sim N(0 , \sigma_b ) $$


$X_{ij}$ is the number of item i sold from shop j
$\alpha_i$ is the change in number of mean sales associated to item i
$\beta_i$ is the change in number of mean sales associated to shop j
$\mu_0$ is some global average mean sales number

I would then start by placing flat priors on $\mu_0$ , $\sigma_a$ and $\sigma_b$. Obviously this is quite a basic model in that it assumes complete independence of sales between shop and item.

You could extend this by letting $\alpha_i$ and $\beta_j$ come from a bi-variate normal distribution where you attempt to model the correlation between the two, for example if a shop is particularly amazing at selling items they might have less of a per item difference. Or alternatively you might want to include an interaction term as certain shops might be bad at selling certain items and you might want to capture this.


The above model has 2 glaring flaws though, the first being that all items share a common variance parameter which is just highly unlikely for example imagine the variance of the number of tv sales vs the number of apple sales. The second is that currently we fit a flat increase/decrease per shop regardless of what the item is. In order to attempt to address these issues we could fit the following model:

$$ X_{ij} \sim N( \mu_{ij} , \sigma_i ) \\ \mu_{ij} = \mu_0 + \alpha_i + \beta_{ij} \\ \begin{bmatrix} \beta_{1j} \\ \beta_{2j} \\ \vdots \\ \beta_{Ij} \end{bmatrix} \sim N_I( 0 , \Sigma) \\ \sigma_i \sim \Gamma(a,b) \\ \alpha_i \sim N( 0 , \sigma_a) \\ $$

Under this model you get a global mean for each individual item as $ \mu_0 + \alpha_i$ and then some adjustment to this mean based upon which item it is and which shop its being sold from $\beta_{ij}$. The specification of $\Sigma$ defines the relationship and correlation between each item and shop. A high covariance would imply that shops that are good at selling 1 item are good at selling others.

  • $\begingroup$ Amazing, thanks! Unfortunately I won't be able to code this up immediately, as this problem is not priority in my project, but I've encountered it several times (in different guises) and everytime I wondered how to solve this. I hope I do manage to implement it soonish. :) One more question that isn't immediately clear to me, how does this handle missing data? e.g. some products aren't available in all shops. $\endgroup$
    – Toekan
    Mar 19, 2018 at 21:56
  • $\begingroup$ Hi @Toekan, that is a very good question and one I don't have a good answer for. In short it doesn't handle missing data. A simple solution is to just assign the items the shops don't sell to 0 though this may negatively bias your global mean for each item. Another option could be to model each item independently only including the stores that do sell it however you then lose out on the correlation information (though this isn't an issue if you arn't interested in the performance of the individual shops). Alternatives are to try and impute the missing data for example MICE. $\endgroup$
    – gowerc
    Mar 20, 2018 at 18:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.