  Several methods have been proposed in the spatial statistics literature for
the analysis of big data sets in continuous domains. However, new methods for
analyzing high-dimensional areal data are still scarce. Here, we propose a
scalable Bayesian modeling approach for smoothing mortality (or incidence)
risks in high-dimensional data, that is, when the number of small areas is very
large. The method is implemented in the R add-on package bigDM. Model fitting
and inference is based on the idea of "divide and conquer" and use integrated
nested Laplace approximations and numerical integration. We analyze the
proposal's empirical performance in a comprehensive simulation study that
consider two model-free settings. Finally, the methodology is applied to
analyze male colorectal cancer mortality in Spanish municipalities showing its
benefits with regard to the standard approach in terms of goodness of fit and
computational time.
