Geographic Disparity in Chronic Obstructive Pulmonary Disease (COPD) Mortality Rates among the Taiwan Population

Chronic obstructive pulmonary disease (COPD) causes a high disease burden among the elderly worldwide. In Taiwan, the long-term temporal trend of COPD mortality is declining, but the geographical disparity of the disease is not yet known. Nationwide COPD age-adjusted mortality at the township level during 1999–2007 is used for elucidating the geographical distribution of the disease. With an ordinary least squares (OLS) model and geographically weighted regression (GWR), the ecologic risk factors such as smoking rate, area deprivation index, tuberculosis exposure, percentage of aborigines, density of health care facilities, air pollution and altitude are all considered in both models to evaluate their effects on mortality. Global and local Moran’s I are used for examining their spatial autocorrelation and identifying clusters. During the study period, the COPD age-adjusted mortality rates in males declined from 26.83 to 19.67 per 100,000 population, and those in females declined from 8.98 to 5.70 per 100,000 population. Overall, males’ COPD mortality rate was around three times higher than females’. In the results of GWR, the median coefficients of smoking rate, the percentage of aborigines, PM10 and the altitude are positively correlated with COPD mortality in males and females. The median value of density of health care facilities is negatively correlated with COPD mortality. The overall adjusted R-squares are about 20% higher in the GWR model than in the OLS model. The local Moran’s I of the GWR’s residuals reflected the consistent high-high cluster in southern Taiwan. The findings indicate that geographical disparities in COPD mortality exist. Future epidemiological investigation is required to understand the specific risk factors within the clustering areas.


Introduction
According to the global burden of disease reported by the World Health Organization [1], chronic obstructive pulmonary disease (COPD) was the fourth leading cause of death in 2004, and caused 3 million deaths, which corresponded to 5.1% of all deaths globally. The prevalence of COPD worldwide was around 63.6 million cases, and it was also an important leading cause of moderate to severe disability among the elderly. One projection showed that COPD will become the third leading cause of death in 2020 [2]. In Taiwan, the temporal trend of the age-adjusted combined bronchitis, emphysema and asthma mortality rate has been declining sharply, from 33.2 per 100,000 in 1981 to 3.9 per 100,000 in 2007 [3]. However, the specific COPD International Classification of Diseases, Ninth Revision (ICD-9) code 496 was not included in computing COPD-related mortality. Thus, overall COPD-related mortality was greatly underestimated in Taiwan.
The estimated prevalence of moderate to severe COPD among adults aged §30 years was 5.4% in Taiwan, while the average prevalence in the 12 Asia-Pacific countries including Taiwan was 6.3%, and in the United States it was 4.8% [4]. With the increasing aging trend for the next few decades, the COPD prevalence among the elderly is expected to surge, and will cause a critical impact on health care and long-term care systems. One economic estimate of the direct medical cost and indirect cost such as lost work and productivity caused by COPD was that it cost a total of $38.8 billion in 2005 in the United States [5]. The higher severity of COPD hospitalized patients resulted in greater medical costs than those for patients with moderate conditions [6]. Thus, combating the threat of COPD will be the top priority of public health. The major risk factors of COPD are smoking [7], exposure to biomass smoke [8], air pollution [9,10], age [11], gender [12], socio-economic status [13], tuberculosis exposure [14][15][16] and altitude [17]. In addition to these risk factors, spatial clustering and geographic disparities are important information for public health intervention. In the United States, as of 2011 the COPD hospitalization rates were clustered in the eastern and south central areas [18]. In China, the COPD mortality has been shown to be significantly higher among the people living in northern and rural China [19]. Thus, it is important to elucidate the geographical distribution of COPD and figure out possible high risk areas for further epidemiologic investigation and prevention.
In the past, COPD-related studies in Taiwan rarely mentioned the geographical distribution of COPD mortality and possible health disparities. In addition, we will use the complete disease definition of COPD to explore COPD-related mortality. The purpose of this study was to analyze the spatio-temporal clusters of COPD-related mortality and to explore their correlations with the local-specific risk factors.

Ethics
This study was approved by the institutional review board (IRB) of Academia Sinica (IRB#: AS-IRB-BM 13057). The databases we used were all stripped of identifying information and thus informed consent was not needed.

Age-adjusted Mortality Rates
The cause of death data at the township level from 1999 to 2007 were collected from the Department of Health in Taiwan, and the mid-year population data were acquired from the Census database of the Ministry of the Interior. The definitions of COPD deaths used the International Classification of Diseases, Ninth Revision (ICD-9) codes 490 (bronchitis, not specified as acute or chronic), 491 (chronic bronchitis), 492 (emphysema), and 496 (chronic airway obstruction, not elsewhere classified). The direct age adjustment method was applied for age-adjusted mortality rates. The reference population was the year-2000 Taiwan population. In addition, direct township-level data on pulmonary tuberculosis incidence and prevalence in the 1990's were not available. Thus, we used 1994-1999 age-adjusted pulmonary tuberculosis (ICD-9: 011) mortality as a proxy indicator for tuberculosis exposure. The reason why we selected the period 1994-1999 was because the administrative units could be matched with the current spatial units; otherwise the spatial unit might not be consistent for further analysis.

Smoking Rates
The smoking rates were obtained from the earliest National Health Interview Survey in 2001, which mainly collected information on residents' health care utilization, expenditure and their health behaviors by a nationwide multistage stratified systematic sampling scheme. Males' and females' smoking rates were calculated separately. The spatial unit was the city or county level. Therefore, the heterogeneities within a city or county would be underestimated due to the limitation of the data.

Air Pollution
The daily air pollution level from 1994 to 1999 was obtained from 73 Environmental Protection Administration (EPA) monitoring stations in Taiwan. The annual average concentration of ambient PM with an aerodynamic diameter !10 mm (PM 10, mg/ m 3 ) and the gaseous pollutants SO 2 (ppb), NO 2 (ppb) , CO(ppm) were computed for each station. Then, we used the empirical Bayesian Kriging [20] in ArcGIS (ArcMap, version10.2; ESRI Inc., Redlands, CA, USA) to interpolate the concentrations throughout Taiwan, and used the zonal statistics function to calculate the average concentration for each township.

Densities of Health Care Facilities
Because of the privacy issue, we could not get the patients' exact home address information or which hospitals the patients went to. Thus, we calculated the mean number of the contracted hospitals or clinics in each township from the National Health Insurance Database, which covered nearly 99% of the population in Taiwan [21], as the numerator of the density of health care facilities. To take the size of the different townships into account, the denominator was the size of each township. The density of health care facilities was our proxy for representing the indirect accessibility index.

Area Deprivation Index and Percentage of Aborigines
The indicators of social economic status (SES) for each township here were area deprivation index and percentage of aborigines. We used population and housing census data in 2000 from the Survey Research Data Archive, Academia Sinica, Taiwan. The area deprivation index was composed of two parts, including the standardized proportion of elementary occupations and the standardized drop-out rate of students aged 15-17 from schooling [22]. The percentages of aborigines were directly computed from the census data.

Average Altitude in Each Township
The digital elevation data were downloaded from the ASTER GDEM website (http://gdem.ersdac.jspacesystems.or.jp/search. jsp), which is a publicly available dataset for global digital elevation with 30 m * 30 m resolution [23]. Then, we used the zonal statistics function of ArcGIS to calculate the average altitude in each township. The spatial distribution of original digital elevation and the average altitude are shown in Figure S1.

Statistical Analysis
To understand the relationship between COPD mortality and possible risk factors including smoking, air pollution, density of health care facilities, area deprivation index, pulmonary tuberculosis mortality, the percentage of aborigines and the altitude, we used Geographically Weighted Regression (GWR) to estimate and consider the spatial variability of their relationships [24]. Unlike conventional ordinary least squares (OLS), GWR model is a type of local statistic in which the parameter estimations vary over space. The assumption of spatial non-stationarity is made in GWR, which means the correlations between the independent variable and dependent variables are not the same for every area [25]. In the first step, we apply OLS analysis with SPSS 20.0 (IBM Corp., Armonk, NY, USA) to observe the correlations between COPD mortality and risk factors. Then, we then use free software GWR 4.0 to run the GWR analysis (http://www.st-andrews.ac. uk/geoinformatics/gwr/gwr-downloads/) [26]. The coefficients of each explanatory variable and the R-square in both OLS and GWR will be summarized in the results. The residual maps and the spatial distribution of the significant explanatory variables after GWR will be displayed in the results with ArcGIS (ArcMap, version10.2; ESRI Inc., Redlands, CA, USA). Finally, we use global Moran's I, which is a global test statistic for spatial autocorrelation. The interval of global Moran's I is between-1 and 1. A higher positive Moran's I indicates that values in the neighboring areas tend to cluster, while a lower negative Moran's I implies that higher and lower values are interspersed. When Moran's I is close to 0, there is no spatial clustering, meaning that the data are randomly distributed [27]. We tested the global spatial autocorrelation of age-adjusted COPD mortality and their residuals after GWR analysis. A Local Moran's I (LISA) cluster map of the residuals was created for identifying the clusters which cannot be explained by the current risk factors.

Results
During 1999-2007 (Figure 1), the COPD age-adjusted mortality rates for males declined from 26.83 to 19.67 per 100,000 population, and those for females declined from 8.98 to 5.70 per 100,000 population. Overall, males' COPD mortality rate was around three times higher than females'. In Figure 2, the ageadjusted mortality rates were found to be clustering in southern Taiwan (Males: Moran's I = 0.22, p,0.05; Females: Moran's Although the temporal trend of COPD mortality in Taiwan was going down, geographical disparities still occurred in some townships. Thus, we applied the GWR to identify the effects of risk factors and possible clusters. In Table 1, the ordinary least square model was applied to model the relationship between COPD age-adjusted mortality and the risk factors. For the males, the male smoking rate, area deprivation index, percentage of male aborigines, CO and PM 10 all have significant positive correlations with COPD morality (p, 0.05). The male tuberculosis mortality, densities of health care facilities, NO 2 and altitude have significant inverse correlation with COPD mortality (p,0.05). The overall adjusted R-square was 49.4%. For the females, the percentage of female aborigines, CO and PM 10 all have significant positive correlations with COPD morality (p,0.05). Female smoking rate, female pulmonary tuberculosis mortality, NO 2 and the altitude have significant inverse correlation with COPD mortality (p,0.05). The overall adjusted R-square was 60.9%. The collinearity was checked by the variance inflation factor (VIF), and the values were all below 7.5, which indicates that there is no collinearity problem.
In Table 2, the summary statistics of the coefficients for every township are listed. For the males, the medians of all factors' coefficients are positively correlated with male COPD mortality except for densities of health care facilities and SO 2 . The overall adjusted R-square was 72.2%. For the females, the medians of female smoking rate, percentage of female aborigines, SO 2 , PM 10 and the altitude were all positively correlated with female COPD mortality. The overall adjusted R-square was 77.9%.
The residual maps after GWR are shown in Figure 3. For both males and females, the spatial autocorrelation of the residuals was not statistically significant, which indicated there was no global spatial clustering after GWR. Global Moran's I indexes were2 0.008 (p = 0.60) for males and 20.007 (p = 0.65) for females, respectively. Furthermore, we used local Moran's I to find out the local high-high clusters of residuals. For both the males and the females (Figure 4), the high-high clusters were all in southern Taiwan.
In order to observe the significant risk factors' effects on the COPD mortality rates, we displayed the selected geographical distribution of the four major risk factors' coefficients in Figure 5 and Figure 6. For the males (Figure 5), the effect of the male smoking rate was higher in central Taiwan, but the percentage of male aborigines, CO and PM 10 had higher effects in southern Taiwan. For the females (Figure 6), the female smoking rate had a higher effect in central Taiwan, but the percentage of female aborigines and PM 10 had higher effects in southern Taiwan.

Discussion
This paper presents a 9-year spatial-temporal distribution of age-adjusted chronic obstructive pulmonary disease mortality, and identifies possible cluster areas for future epidemiologic investigation. To our knowledge, this is the first paper to address the spatial and health disparity issue of COPD mortality in Taiwan. This study can therefore be a pioneer project laying out the foundation and scheme for broader and longer COPD analysis in Taiwan. The official report of cause of the deaths only reported bronchitis, emphysema and asthma before 2007. The specific diagnosis of chronic airway obstruction (ICD-9: 496) was not included. However, the total number of deaths from chronic airway obstruction during 1999-2007 was 6.8 times more than those from bronchitis (ICD-9: 490), chronic bronchitis (ICD-9: 491) and emphysema (ICD-9: 492). Therefore, we included all these four diagnoses into our disease definition of COPD [28]. The overall temporal trends of COPD mortality have been declining in Taiwan. In the United States, the COPD mortality in males also declined during 2000-2005, but increased in females [29]. The COPD mortality in the Southeast Asia region was estimated at 40.0 per 100,000 population and 79.8 per 100,000 for the Western Pacific region [30]. The average age-adjusted COPD mortality during 1999-2007 in Taiwan was 15.88 per 100,000, which was lower than the average of those two regions. Two possible reasons explaining the declining trend in Taiwan were that the 99% health insurance coverage reduced the economic barrier to health care, and there was a declining trend of smoking rates in the past decade [31]. The conventional OLS model explained total variance of 49.4% for the males and 60.9% for the females. After applying the GWR model, the adjusted R-square increased 22.8% for the males and 17% for females (Males: 72.2%; Females: 77.9%). The global Moran's I indexes of the residuals were not significant, and the Moran's I indexes were close to 0, which meant nearly no clustering pattern. With the local Moran's I, the consistent clusters for both genders were distributed in a few townships of southern Taiwan.
In Figure 3 and Figure 4, the residuals for both genders were high and clustering in southern Taiwan, especially in Pingtung county. There are some possible reasons for these clusters. The most important might be exposure to tuberculosis. Pingtung county had the highest TB incidence (97.5 per 100,000 population) and the second highest TB mortality rate (6.2 per 100,000) in Taiwan in 2011 [32]. It also had the highest TB incidence among the population aged 55-64 (128.0 per 100,000 population) and §65 (456.4 per 100,000 population) [32]. However, a longer delay on TB diagnosis and treatment was observed in southern Taiwan [33], which might enhance the risk of getting COPD [15]. Information on the exposure rates to TB in the townships during the 1990's or before the 1990's was lacking.     Therefore, we used 1994-1999 TB age-adjusted mortality as a proxy indicator. However, an inverse correlation was found between TB mortality and COPD mortality. That might have resulted from the competing risks of death among the elderly [34]. The distribution of age-adjusted pulmonary tuberculosis mortality rates during 1994-1999 is shown in Figure S2. The effect of the smoking rate was lower in southern Taiwan, which might be underestimated, because the survey results were at the city or county level rather than the township level. The heterogeneity among the townships was unknown, and this was our major limitation. In the southern clusters, nearly 90% of the residents were aborigines, who had a high prevalence of health risk factors such as smoking, drinking and betel nut chewing, causing higher mortality and reducing life expectancy [35].
The socioeconomic status was found to affect the risk of developing COPD and COPD mortality [13]. In our GWR results, the lower densities of health care facilities (both genders), the higher area deprivation index (females), and the higher percentage of aborigines (both genders, Figure S3) were all correlated with an increase in COPD mortality. Although the health insurance coverage reached 99% in Taiwan, access to health care facilities in rural villages was still low. The average socioeconomic status was relatively low among the mountainous aborigines [36]. These factors might hinder proper management of the disease course and result in a worse prognosis. The same pattern was also observed in China, where COPD mortality rates were higher in rural areas than in urban ones (relative risk: 2.14 [1.86-2.46]) [19].
Although ambient air pollution was not conclusively shown to cause COPD [10], the epidemiological data showed that outdoor air pollution had short-term and chronically adverse effects on COPD-related hospital admission and mortality, especially for the elderly [37,38]. In Taiwan, the higher levels of ambient pollutants increase the risk of hospital admissions [39,40]. In this study, we found the air pollutants such as CO, PM 10 and SO 2 had positive correlations with COPD mortality, and NO 2 had negative correlations with COPD mortality. However, the highest average concentrations of CO were in northern Taiwan ( Figure S4), but these might be correlated with other high concentrations of air pollutants in southern Taiwan, where the clusters of COPD mortality were found.
A previous study in the United States [17] found that COPD mortality had a dose-response correlation with the increasing altitude. The possible biological mechanism was that for persons with respiratory problems, a lower ambient oxygen level might exacerbate hypoxia, pulmonary hypertension and even cause death. In our findings, the OLS results showed a mild inverse correlation between COPD mortality and altitude because most townships (316/358, 88.3%) are located at altitudes less than 500 m. The median values and third quartiles of the GWR estimation in the altitude were positive correlation. This reflects the fact that some townships had positive correlation but some did not. Thus, we further divided the altitudes into four groups (0-500 m, 501-1,000 m, 1001-1500 m, 1501+ m) and compared the COPD mortalities among them (Table S1). It was found that the lowest altitude group had lower COPD mortalities than the other groups. However, there was no clear dose-response relationship due to fewer townships being located at higher altitudes. Therefore, altitude's effect on COPD mortality was not conclusive in this study.

Limitations
This study demonstrated some correlations between health behavior, socio-economic status, air pollution and COPD mortality with ecological data. There were some limitations in this study. The first one is lack of a precise smoking rate and tuberculosis prevalence data at the township level and in the 1980's or 1990's. These two indicators were declining overall as a temporal trend. Thus, the long-term cumulative effects on COPD mortality rates might be underestimated. In addition, the ambient outdoor air pollution was averaged as one concentration of each air pollutant in each township. The exacerbation of COPD mortality resulting from short-term exposure to elevated concentrations could not be observed from this model. Furthermore, we did not have any information on indoor biomass fuel use, especially from the clustering areas. Therefore, it will be worth conducting further epidemiological studies on those clustering areas to know the specific risk factors locally.

Conclusion
This study displays the temporal and spatial distribution of COPD mortality in Taiwan. After applying a GWR model and adjusting the known and available risk factors such as smoking, socio-economic status, and air pollution, the spatial variability of COPD mortality in most of the townships can be explained. There are still some abnormal clusters in southern Taiwan. This will need further research to disentangle the local risk factors in order to reduce the disease burdens there. Figure S1 The spatial distribution of altitude in Tai