  Despite the amount of research on disease mapping in recent years, the use of
multivariate models for areal spatial data remains limited due to difficulties
in implementation and computational burden. These problems are exacerbated when
the number of small areas is very large. In this paper, we introduce an
order-free multivariate scalable Bayesian modelling approach to smooth
mortality (or incidence) risks of several diseases simultaneously. The proposal
partitions the spatial domain into smaller subregions, fits multivariate models
in each subdivision and obtains the posterior distribution of the relative
risks across the entire spatial domain. The approach also provides posterior
correlations among the spatial patterns of the diseases in each partition that
are combined through a consensus Monte Carlo algorithm to obtain correlations
for the whole study region. We implement the proposal using integrated nested
Laplace approximations (INLA) in the R package bigDM and use it to jointly
analyse colorectal, lung, and stomach cancer mortality data in Spanish
municipalities. The new proposal permits the analysis of big data sets and
provides better results than fitting a single multivariate model.
