Public health spatial data are often recorded at different spatial scales (or geographic regions/divisions) and over different correlated variables. Motivated by data from the Dartmouth Atlas Project, we consider jointly analyzing average annual percentages of diabetic Medicare enrollees who have taken the hemoglobin A1c and blood lipid tests, observed at the hospital service area (HSA) and county levels, respectively. Capitalizing on bivariate relationships between these two scales is not immediate as counties are not nested within HSAs. It is well known that one can improve predictions by leveraging correlations across both variables and scales. There are very few methods available that simultaneously model multivariate and multiscale correlations. We propose three new hierarchical Bayesian models for bivariate multiscale spatial data, extending spatial random effects, multivariate conditional autoregressive (MCAR), and ordered hierarchical models through a multiscale spatial approach. We simulated data from each of the three models and compared the corresponding predictions, and found the computationally intensive multiscale MCAR model is more robust to model misspecification. In an analysis of 2015 Texas Dartmouth Atlas Project data, we produced finer resolution predictions (partitioning of HSAs and counties) than univariate analyses, determined that the novel multiscale MCAR and OH models were preferable via out-of-sample metrics, and determined the HSA with the highest within-HSA variability of hemoglobin A1c blood testing. Additionally, we compare the univariate multiscale models to the bivariate multiscale models and see clear improvements in prediction over univariate analyses.
