Population based cross-sectional ecological study, using the “Basic Health Zones” as unit of analysis.
The study was conducted in the Valencia Community, an Autonomous Region on the Mediterranean coast of Spain, with approximately five million inhabitants. Like all of the healthcare system in Spain during the period studied, [22, 23] healthcare coverage is practically universal, being 97% of the population covered by the public Health Service of the Valencia region (the Valencia Health Agency, VHA). The VHA operates an extended network of hospitals (84% of the hospital beds in the region) and primary healthcare centres (PHC). Some key features of the healthcare system during the study period are the following: health care in this network is free of charge (except for drugs in non-retired people who have a 40% co-payment), hospital and primary care is supported by Regional Government budgets, and doctors and other healthcare workers, who enjoy a civil servant-like status, are paid basically by salary. In 2009, the VHA was organized in 23 Healthcare Areas (three of them managed by private companies through Public-Private Partnership agreements) and 240 Basic Health Zones. Healthcare Areas are geographical territories between 75,000 and 500,000 people (most of them between 150,000-250,000 people) served by one public hospital that provides inpatient and outpatient specialized care to the BHZs of its demarcation. BHZs are small geographical areas between 1,000 and 65,000 people (most of them between 10,000 and 25,000 people) commonly served by one Primary Healthcare Centre with a stable team of doctors, nurses and other healthcare workers. Due to these organizational characteristics (geographical planning, minimal accessibility barriers, and practical absence of economic incentives to providers), patients receive most of their primary care from the Primary Healthcare Center of the BHZ where they belong, and most of their specialized care, including outpatient consultations, from the hospital of the corresponding Healthcare Area.
Sources of data
The Population Information System, called SIP, a record of the population covered by the VHA that assigns an identification number to each individual, supplied the population denominator. Among other information, this dataset provides the Healthcare Area, the BHZ and Primary Healthcare Centre to which each individual belongs. The numerator (PPH admissions), was built using anonymized data from the Minimum Basic Hospital Discharge Dataset of the VHA from 2007 to 2009. This database, similar to the US Uniform Hospital Discharge Dataset, provides clinical and sociodemographic information on all hospital discharges in the VHA, and, basically, it is a synopsis of the patient episode of care, including diagnoses and procedures coded according to the International Classification of Diseases 9th revision Clinical Modification (ICD9CM). The VHA Hospital Discharge Dataset includes the patient’s BHZ of residence and, once anonymized, was transferred to the research team.
The study population consisted of all residents 15 years and over that were registered in the SIP register in the period 2007–2009. The data was aggregated by five-year age-sex groups and BHZ.
Main outcome measure
Age-sex standardized rates by 100,000 person-years in each BHZ, of six chronic PPH: diabetes short-term complications, chronic obstructive pulmonary disease (COPD), congestive heart failure (CHF), dehydration, angina admission and adult asthma. All 2007 to 2009 hospital admissions -readmissions included- of patients aged 15 years and over with a main diagnosis one of these PPH were selected and aggregated by age-sex groups and BHZ. For PPH operative definitions, we used the criteria of the Spanish validation [24, 25] of the US Agency for Healthcare Research and Quality (AHRQ) Prevention Quality Indicators . This Spanish version is similar to the US version, but some ICD9CM codes were adapted to the codification patterns most common in Spain and are fully described in a previous work [25, 26].
This study, observational in design, uses retrospective anonymized non-identifiable and non-traceable data provided upon request by the Health Department of the Valencia Regional Government (not the similar and freely available MBDS from the Spanish Ministry of Health). The authors declare that the transfer of the data to the research team met the requirements of the provider and, according to the CIOMS-WHO International Ethical Guidelines for Epidemiological Studies  and the Spanish personal data protection  and patients rights’ laws,  did not require Ethics Committee approval.
First, age-standardized rates were obtained for each BHZ and each single condition, together with global PPH frequencies and global age-standardized PPH rates. Variability among rates was quantified using the Extremal Quotient, excluding areas outside the percentiles 5 and 95 (EQ5-95), and the Empirical Bayes statistic (EB). Standardized Hospitalization Ratios (SHR) were estimated for each condition using the ratio of observed-to-expected cases (oij/eij, being oij and eij the observed and expected number of cases for BHZ i and PPH condition j), and correlations among these were also assessed. The expected number of cases per BHZ unit for each PPH, estimated by applying the rate for the whole region to the population at risk of each BHZ, represents the number of admissions for the condition under study that would have been observed in each BHZ under the hypothesis of constant rate across the whole Valencia region.
To assess the geographical variation in standardized hospitalization ratios, we used the Shared Component Model, which allows to analyze jointly several conditions by decomposing the spatial pattern of each one into two components: one shared by all conditions, and the other that is specific to each one. The SCM has as first-level assumption for the observed counts: Oij ~ Poisson (μij = eijρij), being ρij the unknown relative risk for the BHZ i in the condition j. The second level stage assumes a common structure of risks using a random effect that is shared by the six PPH conditions plus random effects specific to each of them. This is achieved assuming log(μij) = log(eij) + αj + δjφi + ϵij, where αj values are the intercepts for each j-th PPH condition, ϵij are the corresponding specific effects, and φi the random effect representing the shared component of the risk. The δj parameter is a scaling parameter that can be interpreted as a measure of the strength of the association between the shared term and the j-th PPH condition, which is comparable among PPH.
For the specific random effects ϵij, we assumed an exchangeable distribution, whereas for the common random effect φi, we assessed two different specifications: an exchangeable distribution, and a conditional autoregressive distribution (CAR). The hyper-prior specifications used to carried out the Bayesian estimation procedure were: a dflat distribution for αj and a normal distribution N(0,5.9) for log(δj) [21, 30]. For comparison purposes, we also fitted the so-called BYM model  for each condition, which has the same first-level assumption as the SCM models, but use independent random effects for each condition that take into account the spatial correlation through a CAR structure . Model comparisons between the two competing shared component models, based on DIC statistics,  suggested the superiority of the exchangeable distribution (DIC(pD) = 8156.3(887.1)) over the CAR prior assumption (DIC(pD) = 8205.5(937.5)). The SCM exchangeable distribution model has also better DIC than the sum of the DIC values of the individual BYM models, with a difference of 61.5 points, showing a clear advantage of the joint over the individual modeling. Therefore, we selected the SCM exchangeable specification as a definitive SCM model.
All models were implemented in R, version 2.13.1, via the library R2WinBUGS (R Development Core team, 2007), which connects with the software WinBUGS . The estimation procedure was carried out using Monte Carlo methods based on three Markov chains. A total of 49,500 iterations per chain were used, and after a burning period of 12,000 iterations, we kept every 75th for posterior inference. Convergence was determined using the Brooks and Gelman statistic and sequential and autocorrelation graphs. Scripts for SCM and BYM models are given in ‘Additional file 1’.