Several different hierarchical Bayesian models can be used for the estimation of spatial risk patterns based on spatially aggregated count data. Typically, the resulting posterior distributions of the model parameters cannot be expressed in closed forms, and MCMC approaches are required for inference. However, implementations of hierarchical Bayesian models for such areal data are error-prone. Also, different implementation methods exist, and a surprisingly large variability may develop between the methods as well as between the different MCMC runs of one method. This paper has four main goals: (1) to present a point by point annotated code of two commonly used models for areal count data, namely the BYM and the Leroux models (2) to discuss technical variations in the implementation of a formula-driven sampler and to assess the variability in the posterior results from various alternative implementations (3) to give graphical tools to compare sample(r)s which complement existing convergence diagnostics and (4) to give various practical tips for implementing samplers.