Regional aggregates of health outcomes over delineated administrative units such as counties or zip codes are widely used by epidemiologists to map mortality or incidence rates and better understand geographic variation. Disease mapping is an important statistical tool to assess geographic variation in disease rates and identify lurking environmental risk factors from spatial patterns. Such maps rely upon spatial models for regionally aggregated data, where neighboring regions tend to exhibit more similar outcomes than those farther apart. We contribute to the literature on multivariate disease mapping, which deals with measurements on multiple (two or more) diseases in each region. We aim to disentangle associations among the multiple diseases from spatial autocorrelation in each disease.

We propose two Multivariate Directed Acyclic Graph Autoregression (MDAGAR) models using conditional and joint probability laws respectively to accommodate spatial and inter-disease dependence. The hierarchical construction of conditional MDAGAR imparts flexibility and richness, interpretability of spatial autocorrelation and inter-disease relationships, and computational ease, but depends upon the order in which the diseases are modeled. To obviate this, we demonstrate how Bayesian model selection and averaging across orders are easily achieved using bridge sampling. We compare our method with a competitor using simulation studies and present an application to multiple cancer mapping using data from the Surveillance, Epidemiology, and End Results (SEER) Program. We also develop a joint MDAGAR model using latent factors, which avoids the disease ordering issue in conditional modelling.

Based on multivariate disease mapping, one often seeks to identify ``difference boundaries'' that separate adjacent regions with significantly different spatial effects. We adopt a Bayesian multiple-comparison approach for this problem, where we compare all pairs of random effects between neighboring regions. We develop a class of multivariate areally-referenced Dirichlet process (MARDP) models that endow the spatial random effects with a discrete probability law. Within the MARDP framework, the joint MDAGAR model is applied to accommodate spatial and inter-disease dependence for spatial components. We evaluate our method through simulation studies and subsequently present an application to detect difference boundaries for multiple cancers using data from the SEER Program.