Delivering spatially comparable inference on the risks of multiple severities of respiratory disease from spatially misaligned disease count data

Population‐level disease risk varies between communities, and public health professionals are interested in mapping this spatial variation to monitor the locations of high‐risk areas and the magnitudes of health inequalities. Almost all of these risk maps relate to a single severity of disease outcome, such as hospitalization, which thus ignores any cases of disease of a different severity, such as a mild case treated in a primary care setting. These spatially‐varying risk maps are estimated from spatially aggregated disease count data, but the set of areal units to which these disease counts relate often varies by severity. Thus, the statistical challenge is to provide spatially comparable inference from multiple sets of spatially misaligned disease count data, and an additional complexity is that the spatial extents of the areal units for some severities are partially unknown. This paper thus proposes a novel spatial realignment approach for multivariate misaligned count data, and applies it to the first study delivering spatially comparable inference for multiple severities of the same disease. Inference is via a novel spatially smoothed data augmented MCMC algorithm, and the methods are motivated by a new study of respiratory disease risk in Scotland in 2017.

disease risk varies by severity, which is the novel insight provided by this paper.
Our motivating study focuses on respiratory disease in Scotland, and we have data on mild cases treated in primary care, moderately severe cases treated in hospital, and very severe cases resulting in death. Count data for the latter two are available for a set of known nonoverlapping areal units, whereas the primary care data have been spatially aggregated to doctors' surgeries whose catchment areas are partially unknown. Thus, the statistical challenge we address is the change of support problem (Gelfand et al., 2001), because we need to provide comparable inference for all disease severities on a common spatial scale.
Numerous methodologies have been proposed for overcoming spatial misalignment, including downscaling (Song et al., 2014), change of support (Nethery et al., 2019) and regression (Zhu et al., 2003) settings, and designing areal units that minimize spatial aggregation error (Bradley et al., 2017). For count data, Flowerdew and Green (1989) and Mugglin and Carlin (1998) focused on areal interpolation and Bradley et al. (2016) introduced change of support methodology, whereas  and Taylor et al. (2018) produce pseudospatially continuous inference from spatially aggregated data.
The methodological challenge addressed here extends this work to a multivariate setting, because we propose the first spatial realignment model for multivariate misaligned count data, some of which have partially unknown spatial supports. Our Bayesian hierarchical model delivers spatially comparable inference for all disease severities, and captures both spatial and between disease severity correlations via a latent Gaussian process modelled with a multivariate conditional autoregressive (MCAR, Gelfand and Vounatsou, 2003) prior distribution. Inference is based on a novel spatially smoothed data augmented (Tanner and Wong, 1987) Markov chain Monte Carlo (MCMC) algorithm, which extends the algorithm proposed by Taylor et al. (2018) that in our data context, suffers from computational convergence issues. The motivating study is outlined in Section 2, whereas our proposed methodology is presented in Section 3 and tested by simulation in Section 4. The study results are presented in Section 5, whereas Section 6 concludes the paper.

SCOTLAND RESPIRATORY DISEASE STUDY
Our methodology is motivated by a new study of respiratory disease in mainland Scotland in 2017, which includes cases that are: (i) relatively minor and treated in primary care; (ii) moderately severe and require hospitalization; and (iii) severe and result in death. Our modeling aims are to: (a) identify the extent to which the highest risk areas differ by severity; and (b) quantify how the magnitude of health inequalities vary by severity?

Disease data
Mild cases of disease are treated in primary care by doctors grouped within 869 general practice (GP) surgeries, and we obtained the total number of short-acting 2 agonists prescribed by each GP surgery. These medications are used to relieve the symptoms of respiratory disease such as asthma and chronic obstructive pulmonary disease (COPD), and the medications included are listed and justified in Web Appendix A. Following Blangiardo et al. (2016), these data are used as a proxy measure of respiratory disease not requiring hospitalization, because complete incidence data are unavailable. Moderately severe cases of respiratory disease are represented by the numbers of admissions to hospital (ICD-10 codes J00-J99), and these data are spatially aggregated counts for the populations living in each of the 1252 intermediate zones (IZs) that comprise mainland Scotland. IZs are a Scottish Government developed small-area geography, and have an average population of around 4000. Finally, we have data relating to respiratory deaths (again ICD-10 codes J00-J99) for the same set of IZs, and a summary of these count data is provided in Web Appendix A. We account for the differing population sizes and demographics across the spatial units using indirect standardization, which computes the expected numbers of disease events in each unit based on national age-sex-specific disease incidence rates. These expected counts are computed exactly for the hospitalization and death outcomes, because their national age-sex-specific rates are available from Public Health Scotland. In contrast, national rates of respiratory prescribing for different age-sex strata are not publicly available, so we estimate these unknown rates using data on the national age-sex-specific rates of asthma and COPD.
An exploratory measure of disease risk is the standardized incidence ratio (SIR), which is the observed numbers of disease cases/prescriptions divided by their expected numbers. Maps of the spatial distribution of the SIRs for each disease severity are displayed in Figure 1, where an SIR of 1.2 corresponds to a 20% increased risk compared to the national average. In the GP prescription map, each surgery is shown as a dot at the surgery location, because the residential locations and geographical extent of its patient population are partially unknown. The figures show that for all severities, the highest SIR values are in and around the city of Glasgow, which is well known to have some of the poorest health in western Europe.

2.2
A common spatial inferential scale Taylor et al. (2018) make inference on a regular grid when modeling data relating to known areal units, and base their spatial realignment on the area of intersection between each areal unit and grid square. However, the patient catchment area for each GP surgery is partially unknown, which makes it difficult to spatially allocate prescription counts based on areas of intersection. The only data available quantifying the spatial locations of the GP surgery practice populations are how they are distributed across IZs, which is the spatial scale at which we have hospitalization and death data. Specifically, we have data on how many people registered at each GP surgery live in each IZ, and a summary of these population intersections is provided in Web Appendix A. Therefore, as we have disease data for two of the three severities at the IZ level and indirect data relating the GP surgery patient populations to the IZ scale, we use IZs as our common inferential scale. We note that while we know how many people registered with each GP surgery live in each IZ, we do not know which of these individuals obtained prescriptions for respiratory disease. Hence, we cannot directly compute the number of prescriptions at the IZ scale, which necessitates the spatial realignment methodology that we propose. Spatial misalignment such as this can occur routinely when jointly modeling spatially aggregated data from different sources, such as when health care is provided by different providers.

Socioeconomic deprivation data
We quantify the social inequality in disease risk by estimating the effect of socioeconomic deprivation on each severity of disease. Socioeconomic deprivation is measured in Scotland by the Scottish Index of Multiple Deprivation (SIMD), which comprises indicators in the domains of geographical access to services, crime, education, employment, health, housing, and income. However, as our outcome variables are health related, we ignore indicators from the health domain. The remaining 16 indicators have pairwise correlations ranging between −0.81 and 0.98. Thus, we create one measure of socioeconomic deprivation for each domain of the SIMD using principal component analysis, and details are given in Web Appendix A.

Risk model
The first step is to estimate the expected numbers of GP prescriptions  1 () = { 1 ( 1 ), … , 1 ( )} at the IZ level using population-weighted interpolation similar to Flowerdew and Green (1993), which is described in Web Appendix B. In contrast, the observed numbers of GP prescriptions at the IZ scale, { 1 ( )}, are estimated within the inferential algorithm described below. At the IZ level, we model the multivariate disease counts using a Poisson loglinear structure, whose mean depends on covariates and a latent Gaussian process. The latter is commonly used in the spatial modeling literature (Gelfand and Schliep, 2016), but alternatives include the multivariate log-gamma distribution (Bradley et al., 2018) and a Pólya-Gamma augmentation (Bansal et al., 2021). The data likelihood model is given by where ( ) denotes a vector of known covariates. The regression parameters differ by disease severity, and each parameter is assigned a weakly informativeindependent Gaussian prior with mean zero and variance 100,000 to allow the data to speak for themselves.
The random effects = { ( 1 ), … , ( )} 3 ×1 , where ( ) = { 1 ( ), … , 3 ( )}, are modeled with an MCAR prior (Gelfand and Vounatsou, 2003), which captures both spatial and between outcome correlations. MCAR priors are based on a binary × spatial neighbor-hood matrix = ( ), and here = 1 if IZs ( ,  ) share a common border and = 0 otherwise (with = 0 ∀ ). The MCAR prior we use has joint distribution ∼ N[ , { ( , ) ⊗ −1 } −1 ], where 3×3 is the between severity conditional covariance matrix and ( , ) × = {diag( ) − } + (1 − ) is the spatial precision matrix corresponding to the conditional autoregressive prior proposed by Leroux et al. (2000). This prior is used rather than the BYM prior (Besag et al., 1991) because it has an explicit spatial dependence parameter that one can interpret, as well as only having a single random effect for each areal unit (the BYM model has two random effects for each areal unit). The spatial dependence implied by the MCAR model can be seen from its full conditional form for ( )| (− ) (where (− ) = ⧵ ( )) given by (2) Here, is a spatial dependence parameter, with = 1 corresponding to strong spatial dependence (the conditional expectation of ( ) is the mean of the random effects in neighboring IZs), while if = 0, then ( ( ), ( )) are independent. We specify a noninformative uniform prior for on the unit interval, that is, ∼ Uniform(0, 1), and a weakly informative Inverse-Wishart(4, ) prior for the cross-severity covariance matrix .

Inference
We fit the model using a spatially smoothed dataaugmented MCMC algorithm, which jointly updates the IZ-level model parameters = { 1 , … , 3 , , , } and the IZ-level GP prescription counts  1 () = { 1 ( 1 ), … , 1 ( )} conditional on the disease and covariate data  and the population data . Our algorithm extends the proposal of Taylor et al. (2018) that is summarized in Web Appendix B, and iterates the following steps.
(1) Sample from { |  1 (), , }, the conditional distribution of the model parameters given the current values of the IZ-level prescription counts  1 () and the data {, }.
Step (2) is the data augmentation step, and we undertake it at every th iteration of the MCMC algorithm because this lessens the dependence between successive values of used in this step. As described in Web Appendix B, we sample from {  1 ()| , , } using multinomial sampling steps, because each GP surgery's total number of prescriptions 1 ( ) has to be allocated to the IZs in an integer fashion. Letting 1 ( ∩  ) denote the number of prescriptions allocated to IZ  from the patients registered to GP surgery  , step (2) involves sampling from The IZ-level prescription counts are then computed as which allocate GP prescriptions to IZs proportionally to the product of the population overlap via ( ∩  ) and the IZ-level risk via 1 ( ). The latter are the current values obtained from step (1) of the algorithm, and most * = 0 as ( ∩  ) = 0 for most combinations of GP surgery and IZ. The choice of weights in (4) is not unique, because  and Taylor et al. (2018) proposed slightly different specifications. Our weights are most similar to those of Taylor et al. (2018) outlined in Web Appendix B, although ours are based on population intersection rather than land area intersection.
However, initial analyses showed that these weights can result in either a nonconverging MCMC algorithm or highly inaccurate parameter estimation. For example, in simulated data sets, the relative rates of GP prescriptions at the IZ level {ˆ1( )} were routinely estimated as either very close to 0 or very large at around 20, whereas the true values lay in the interval (0.25, 4). The cause of this problem is having to simultaneously estimate the model parameters and the IZ-level disease counts, and could be phrased as a rich get richer problem. This is because if a large value of 1 ( ) is simulated by chance, then the next sample of (5) This replaces 1 ( ) in (4) with the spatially smoothed alternative˜1( ), which leads to less extreme values of {˜1( )} in the multinomial weights and thus prevents the rich get richer phenomenon. This spatial smoothing is applied to the random effects 1 ( ), because it is these that are assumed to be spatially smooth by (2). Here, ∈ [0, 1] is a spatial smoothing parameter, and using this as a simple spatial smoother was suggested by Banerjee et al. (2004). Clearly setting = 1 in (5) is equivalent to (4). This spatial smoothing step is only applied in (3) and not at any other part of the MCMC algorithm, and we assess its performance in the simulation study in the next section for a range of values of .

SIMULATION STUDY
We present a simulation study to assess: (i) how accurate is the inference from our spatial realignment methodology? and (ii) which value of gives the best estimates? This study provides novel insight into spatial realignment for count data, because Bradley et al. (2016) only considered the ability of the models to estimate fitted values and not covariate effects, whereas Taylor et al. (2018) did not undertake any assessment of inferential accuracy.

Data generation
Data are generated to match the Scotland study with prescription data for = 869 GP surgeries { } and hospitalizaton and death data for = 1252 IZs { }, and our aim is to make inference at the IZ level. Initially, true IZ-level disease risks are generated from ( ) = exp{Educ( ) 1 + Acc( ) 2 + ( )}, where (Acc( ), Educ( )), respectively, denote the access and education domains of the SIMD from the motivating study. The true values are 2 = 0 for all three severities , whereas 1 1 = 0.075, 2 1 = 0.05, 3 1 = 0.025 to match the estimates from the real data. The random effects { ( )} are generated from a zero mean multivariate normal distribution, whose covariance structure is the Kronecker product of a spatial covariance matrix and a between severity covariance matrix. The former is defined by a spatial exponential covariance structure, where the spatial range parameter is fixed so that the correlation between two IZs whose centroids is 1 km apart is 0.9. The between severity correlation matrix is chosen to match the correlations observed in the real data, namely, 0.5 between GP prescriptions and hospitalizations, 0.45 between GP prescriptions and deaths, and 0.9 between hospitalizations and deaths. Finally, the variance of the random effects 2 is allowed to vary in our simulation design.
For the hospitalization and death outcomes, the expected counts ( ) come from the real data, whereas ( ) are generated from (1). In contrast, as the GP prescription data need to be aggregated to the GP surgery level { }, we first generate observed and expected disease counts for the intersection populations { ∩  }. In the majority of these ( ∩  ) = 0, and hence controls disease prevalence and is chosen to match the real data. Finally, we generate which has a similar mean model to (1). Then the observed and expected disease counts at the GP surgery level are computed via Thus, the simulated data being modeled are GP prescription data at the GP surgery level, and hospitalization, death, and covariate data at the IZ level.

Models and scenarios
Models (1)-(2) are fitted to the data with = 0.5, 0.7, 0.9, 0.95, 0.99, 1, where = 1 corresponds to a nonsmoothed data augmentation algorithm. Additionally, we compare two other models, the first being applied directly to the IZ-level data (denoted as IZ), which allows us to quantify the loss in estimation accuracy from having spatial misalignment. The second is a simple spatial realignment approach based on population-weighted interpolation (denoted as Naive), which allows us to evidence the utility of our methods against a simpler alternative. Full details of both these models are given in Web Appendix C. All models are applied to 100 simulated data sets generated under each of six scenarios, which include all pairwise combinations of the following two factors.
• Spatial risk variation. It is controlled by comparing 2 = 0.05 → ( ) ∈ [0.45, 2] and 2 = 0.2 → ( ) ∈ [0.2, 4]. • Accuracy of 1 ( ). In the motivating study, the GP prescription counts { 1 ( )} are estimated rather than computed exactly, due to a lack of data on national agesex-specific GP prescribing rates. We examine the effect that such estimation error has on model performance, by fitting the models using { 1 ( )} calculated with three levels of error. The first of these is no error, whereas the remaining two add zero-mean Gaussian random noise to the true values with standard deviations are = 3 and = 6. Both these correspond to small amounts of error because the mean value of { 1 ( )} is 1500. Inference from each model is based on a single MCMC chain run for 250,000 iterations. The first 50,000 iterations are removed as the burn-in period, and the remaining 200,000 are thinned by 100 to reduce their autocorrelation, resulting in 2000 posterior samples for inference. The data augmentation step is run every 30 iterations, and pilot analyses suggested that these specifications were sufficient for the MCMC algorithm to converge.

Results
The results for the GP prescription outcome are presented here because it is spatially misaligned with the set of IZs used for inference. The results for the hospitalization and death outcomes where there is no spatial misalignment are presented in Web Appendix C for completeness. The accuracy of the IZ-level disease risks { 1 ( )} and the two covariate effects ( 1 1 , 1 2 ) are summarized below by bias and relative mean absolute error (RMAE), the latter being presented as a percentage of the true value. The models produce generally unbiased estimates of disease risk across the scenarios, with biases being less than 0.6 in absolute value in almost all cases. The main exception to this is when = 1, which has four scenarios with larger positive biases (all those where > 0). The percentage RMAEs for { 1 ( )} are displayed in the top panel of Table 1, which shows that as expected having no spatial misalignment in the data (model IZ) leads to the best results (lowest RMAE) in all cases. In contrast, having spatially misaligned data causes around a threefold increase in RMAE when there is relatively little spatial variation in disease risk ( 2 = 0.05), and around a fivefold increase when there is large spatial variation in disease risk ( 2 = 0.2). However, in both cases, the percentage RMAE values are still relatively small, being mostly between 6% and 16% of the size of the true risks.
The Naive interpolation model consistently performs worse than the model proposed here for a range of values of , suggesting that the complexity of our data augmentation approach is warranted. The data augmentation approach used by Taylor et al. (2018) ( = 1) performs comparably with our spatially smoothed variant when there is no estimation error in { 1 ( )} ( = 0), but has much larger RMAEs when the expected counts are unknown exactly ( > 0) as in the motivating study. Typically, large values of less than 1 produce the smallest RMAEs, with = 0.95, 0.99 being best when = 0 and = 0.9, 0.95 generally being best when > 0. The slight exception to this is in the most extreme case where the risk surface is highly spatially variable ( 2 = 0, 2) and the error in { 1 ( )} is largest ( = 6), where smaller values of provide better results. However, in most cases, there are relatively little differences in RMAEs when ∈ [0.7, 0.95], suggesting that the results are relatively robust.
The results for the education variable are included to illustrate whether the models can accurately estimate covariate effects when they are present. Overall, the results show a very similar pattern to those for the risk estimates, and their RMAE values are displayed in the middle panel of Table 1. Unsurprisingly, the best results occur if there is no spatial misalignment (model IZ), whereas the simple interpolation model (Naive) exhibits hugely attenuated estimates close to zero in all cases, which causes the large RMAEs of around 70% of the true value. The standard data augmentation algorithm ( = 1) produces comparable results to our smoothed version if = 0, but exhibits large positive biases and hence large RMAE values when there are errors in the expected counts. The results for our smoothed data augmentation algorithm ( < 1) are not hugely different when ∈ [0.7, 0.95], although = 0.95 typically performs the best across the range of scenarios considered.
The results for the access variable in the bottom panel of Table 1 allow us to investigate whether the models can accurately estimate when a covariate has no relationship with disease risk. As the true value of 1 2 = 0, the results presented are raw mean absolute errors. The Naive interpolation model produces the smallest MAEs, which is artificial because its covariate effect estimates are greatly attenuated toward zero (see above). Our spatially smoothed model performs similarly for ∈ [0.7, 0.95], and its MAE values are not that much larger than those from the model with no spatial misalignment (IZ), suggesting that spatial misalignment has a relatively small effect in this setting. Additionally, in common with the previous results if = 1 and > 0, then the MAE values are greatly inflated compared to our spatially smoothed model. Finally, Web Appendix D presents smaller scale additional simulation studies that investigate the impact of changing additional aspects of this study.

Model fitting
We fit the spatial realignment model to the data with three different covariate combinations: (i) no covariates, (ii) TA B L E 1 Percentage relative (to the true value) mean absolute errors (RMAE) of the regression parameters and the risk estimates relating to the GP prescriptions outcome for each scenario and model. Note that as the true value of the regression parameter 1 2 relating to the access covariate is zero, these results are raw mean absolute errors covariates chosen by a model building strategy (denoted as full), and (iii) each SIMD covariate included in a separate model. The first two of these are included as a sensitivity analysis to see what impact it has on the resulting spatial risk maps (motivating question (a)), whereas the single covariate models allow us to quantify the consistency of the social inequalities in disease risk (motivating question (b)) across the SIMD domains. In building the full model, the education, employment, and income domains of the SIMD are collinear, with Pearson's correlation coefficients above 0.9. Therefore, we only included the education domain from these three, because it minimizes the Akaike information criterion (AIC) when incorporated in a simple Poisson log-linear model in conjunction with the access, crime, and housing variables. We fitted the model with = 0.8, 0.9, 0.95, 0.99, 1 for each of these covariates combinations, and inference in each case is based on 15,000 MCMC samples obtained from three parallel Markov chains. Each chain was run for 1,100,000 iterations with a burn-in period of 100,000, and the remaining samples were thinned by 200 to reduce their autocorrelation. Convergence of the Markov chains for selected parameters was visually assessed using traceplots and numerically assessed using the Gelman-Rubin diagnostic (Gelman et al., 2013), and when = 0.8, 0.9, 0.95, the results show no evidence against convergence. However, when = 0.99, 1, the Markov chains show evidence of nonconvergence, which is visually evident from the traceplots and numerically supported by Gelman-Rubin diagnostics well above 1.1. This lack of convergence is due to the rich get richer phenomenon described earlier, and we illustrate this by presenting the estimated risk map for the GP prescriptions outcome when = 1 in Web Appendix E. Thus, for the rest of this section, we only present the results for models that showed no evidence against convergence, namely, when = 0.8, 0.9, 0.95.

Overall model fit and estimated dependence structures
The overall fit of each model to the data is summarized in panel 1 of Table 2 by the widely applicable information criterion (WAIC), with the estimated number of independent parameters . in brackets. The table shows that the full model provides the best fit to the data for every value TA B L E 2 Summary of the models fitted to the data, including their overall fit as measured by the WAIC and the effective number of independent parameters . (in brackets), as well as posterior means and 95% credible intervals for the correlation parameters

Covariate
Spatial smoothing parameter Quantity model = .
None 29,470 (1580)    of , although the difference in WAIC between it and the model with only the education covariate is very small. Additionally, for each covariate model, = 0.95 provides the best overall fit to the data of the different values considered. Panel 2 in Table 2 displays posterior means and 95% credible intervals for the spatial dependence parameter , which shows that the data contain substantial spatial dependence ( close to 1) even after the covariate effects have been accounted for. The posterior mean does reduce slightly as increases, which is due to the reduced amount of spatial smoothing induced by the data augmentation step (5). The between severity correlations are summarized in panels 3-5 of the table, and display the between severity correlations for the random effects derived from . These correlations are highest for the model with no covariates, which is because in the covariate adjusted models, the covariates account for some of the correlations in the data. For the no covariate model, the GP prescription data are moderately correlated with both hospitalizations and deaths, with correlations between 0.4 and 0.5. In contrast, the correlation between hospitalizations and deaths is higher being around 0.93. This stronger correlation is not surprising, because a death from respiratory disease is often preceded by a stay in hospital for that individual. In contrast, some of the mild cases of respiratory disease (e.g., asthma) will be successfully managed without the patient requiring hospitalization, which explains the lower correlations in this case.

Spatial risk surfaces
The estimated risk surfaces from the full covariate model show little change when varying ∈ {0.8, 0.9, 0.95}, and the biggest differences naturally occur for the GP prescription outcome because these disease counts are estimated rather than known at the IZ scale. The mean absolute differences (over the IZs) between the estimated risk surfaces obtained from ∈ {0.8, 0.9, 0.95} range between 0.038 and 0.085, which compares to disease risk estimates that range between 0.1 and 4.3 across the IZs. The differences between the estimated risks from the full model and the model with no covariates are also relatively small, with mean average differences of 0.126 (GP prescriptions), 0.038 (hospitalizations), and 0.091 (deaths) when = 0.95. Thus, in what follows we present the results from the full model with = 0.95, because it best fits the data as measured by the WAIC and the simulation study showed it performed consistently well across the scenarios considered. The correlation between the posterior mean risk surfaces for hospitalizations and deaths is 0.947, whereas the corresponding correlations between GP prescription and the other two outcomes are lower at 0.635 with hospitalizations and 0.610 with deaths. Figure 2 displays maps of the estimated risk surfaces for the two largest cities Edinburgh and Glasgow for all three severities, and we focus on these cities because displaying Scotland-wide maps as in Figure 1 makes it hard to see the small-scale risk variation. The most striking feature of these maps is that Glasgow exhibits much higher risks than Edinburgh on average for all severities, which is due to the well-known Glasgow effect that Glasgow exhibits some of the poorest health in the western world.
The risk patterns are largely consistent across the three severities, and in Edinburgh, the high-risk areas are mainly in Craigmillar in the east and Broomhouse in the west. In Glasgow, the main high-risk areas include Drumchapel in the west, Castlemilk in the south, and a large part of the east end of the city. The other striking feature of the Glasgow maps is that visually, there appears to be more high-risk areas as the severity of disease increases, with the high-risk areas for GP prescriptions expanding when increasing in severity to hospitalizations and then deaths. This is most apparent in the north-east of the city in locations such as Riddrie, which has a risk 22% below the Scottish average (ˆ1( ) = 0.78) for GP prescriptions, compared to risks that are 59% above (ˆ2( ) = 1.59) and 61% above (ˆ3( ) = 1.61) the Scottish average for hospitalizations and deaths, respectively.

Health inequalities
The inequalities in disease risk are summarized in Table 3 for all severities, where Panel (A) presents total inequality, while Panel (B) displays social inequality. Total inequality is summarized by the standard deviations and interquartile ranges in the posterior mean risk surfaces, which come from the full model with = 0.95. The table shows that the level of total inequality reduces as the severity of disease increases, which is a consistent finding across both standard deviation and interquartile range metrics. For example, the interquartile range in risk is 0.602 for GP prescriptions, 0.579 for hospitalizations, and 0.469 for deaths, a reduction of 22% in the inequality of deaths compared to the inequality in GP prescriptions. The level of social inequality is summarized in Panel (B) of Table 3 by the relative risk between each SIMD covariate and disease risk. All results relate to = 0.95 and come from single covariate models, and to ensure comparability across the different SIMD domains, all relative risks relate to a one-standard-deviation increase in the covariate in question. For completeness, the estimated covariate effects from the full model are presented in Web Appendix E. With the exception of the access domain that specifically F I G U R E 2 Maps displaying the posterior mean risk estimates for Edinburgh (left) and Glasgow (right) for GP prescription rates (top), hospitalizations (middle) and deaths (right). This figure appears in color in the electronic version of this article, and any mention of color refers to that version TA B L E 3 Summary of the health inequalities in respiratory disease risk for all three severities. Panel (A) presents total inequality, whereas Panel (B) displays social inequality measures geographical access to services rather than general socioeconomic deprivation, all estimated relative risks and corresponding 95% credible intervals are greater than one. This shows consistent convincing evidence of social inequality in disease risk for all severities, as increases in the levels of socioeconomic deprivation are associated with significant increases in disease risk. The sizes of these effects are large and range between a 20.5% increased risk of hospitalization for the crime indicator, to a 50.2% increase in GP prescription rates when the income deprivation indicator increases by one standard deviation. The levels of social inequality seem to be highest for the GP prescription outcome, because it has the largest inequality for the crime, education, employment, and income domains of the SIMD.

DISCUSSION
This paper develops a new approach for modeling multivariate spatially misaligned disease count data with known and partially unknown spatial supports, which delivers inference on a common spatial scale using a spatially smoothed data-augmented MCMC algorithm. The simulation study shows that the spatial smoothing in the data augmentation step provides reliable inference across a range of scenarios, and consistently outperforms the data-augmented MCMC algorithm proposed by Taylor et al. (2018) in the data context considered here. The study suggests that ∈ [0.9, 0.95] provides the most accurate inference, but that the results do not vary greatly for values down to = 0.7. This preference for a small amount of spatial smoothing is because as gets closer to zero, the current estimate of risk used in the data augmentation step is overly smooth, leading to oversmoothing and poorer estimation in the IZ-level risk estimates. In contrast, as approaches 1, this oversmoothing reduces, leading to better estimation. However, if gets very close to 1, the MCMC algorithm is affected by the rich get richer phenomenon, leading to inaccurate parameter estimation. In the simulation study, the presence of spatial misalignment reduces inferential accuracy as expected, with the RMAE for disease risk increasing by between three-and fivefold depending on the scenario, whereas the corresponding result for covariate effects is only a 1.5fold increase. However, for disease risk, the absolute size of the RMAE in the presence of spatial misalignment is only between 6% and 16% of the true risk size, suggesting that reliable inference can still be obtained from spatially misaligned data. The novel insight provided by this paper is how disease risk varies by the severity, which extends the existing literature that ubiquitously focuses on a single severity of outcome. The first key finding from our Scotland respiratory disease study is that the spatial risk trends are largely consistent across the different severities, with mainly the same subregions exhibiting elevated risks for each severity. This is likely to be because areas that exhibit elevated risks will have increased exposures to factors such as smoking and air pollution, which thus affect all severities of outcome. Additionally, elevated risk areas are likely to self-perpetuate across severities, because having a larger proportion of people suffering from mild disease is likely to lead to this cohort having higher rates of more severe disease.
The second key finding from our Scotland study is that both total and social health inequalities reduce as the severity of disease outcome increases, with mild disease treated in primary care being much more unevenly distributed across different sections of society than moderate cases requiring hospitalization or severe cases resulting in death. This reduced inequality for deaths is likely to be because everybody dies eventually, and as respiratory disease is one of the most common causes of death, it affects people from all communities. In contrast, respiratory diseases such as asthma that can be treated in primary care are more likely to affect younger people (see table 8.2 in Scottish Government, 2018), and the incidence rates of such cases tend to be driven by risk factors such as those highlighted above whose magnitude varies spatially. However, this examination of social health inequalities was based on the estimated effects of covariates, which were obtained from a data set containing residual spatial autocorrelation. This autocorrelation was modeled by spatially correlated random effects as is standard practice in the literature, but Hodges and Reich (2010) have shown that confounding can happen in this context because the random effects can themselves be correlated with the covariates. A number of approaches such as Hughes and Haran (2013) have been proposed to deal with this potential issue when there is no spatial misalignment in the data, but the extension of this work to our spatial misalignment setting has not been explored and is an important avenue for future work.
Finally, this paper represents the start of a long-term focus on spatially misaligned count data modeling, which has received much less attention in the literature to date compared to spatial misalignment in a continuous data setting. Our use of IZs here as the common inferential scale was a pragmatic choice based on data availability, but may suffer from the well-known modifiable areal unit problem (MAUP, Openshaw and Taylor, 1979).  and Taylor et al. (2018) overcome this limitation by providing pseudospatially continuous inference on a regular grid, but they mainly focus on the simpler setting where the areal units have known spatial extents. In our study, the GP surgery catchment areas have partially unknown spatial supports, precluding the use of area of intersection in the spatial misalignment mechanism as used by the above authors. Therefore, future work will begin by investigating how to deliver pseudospatially continuous inference on a regular grid in the context of areal units with partially unknown spatial supports.

A C K N O W L E D G M E N T S
The authors gratefully acknowledge the contributions of two reviewers and an associate editor, whose suggestions have improved the motivation for and content of this paper.

D ATA AVA I L A B I L I T Y S TAT E M E N T
The General Practice prescription and Scottish Index of Multiple Deprivation data are available at https://www.opendata.nhs.scot/dataset/prescriptions-inthe-community and https://www.gov.scot/publications/ scottish-index-multiple-deprivation-2016/ , while the hospitalization and death data can be requested from Public Health Scotland (https://publichealthscotland.scot/).