Space–Time Relationship between Short-Term Exposure to Fine and Coarse Particles and Mortality in a Nationwide Analysis of Korea: A Bayesian Hierarchical Spatio-Temporal Model

Abstract

Since air pollution is particularly detrimental to vulnerable people, such as senior citizens and infants, understanding how ambient air pollution affects health is very important for making relevant policy [2,3]. [...]there is a spatial confounding bias because space-time varying air pollution covariates and space-time random effects are simultaneously included in space-time health modeling [17]. In this paper, to better estimate air pollution effects without the influence of spatial confounding bias, we adopted the two-stage model proposed in [15]. [...]most of the epidemiologic studies of air pollution conducted in South Korea have focused on how temporally varying air pollutant concentrations affected mortality in individual cities [18,19], but such studies should be conducted over entire areas of Korea to establish national-level air pollution control policy for health improvement. [...]the random component ϕit is the space-time interaction term.

1. Introduction

Air pollution has been viewed as a threat to human health since the onset of rapid industrialization. It is widely known that air pollutants, including particulate matter with an aerodynamic diameter less than 10 µm (PM₁₀) or 2.5 µm (PM_2.5; fine particulate matter), is linked to total, cardiovascular, and respiratory mortality [1]. Since air pollution is particularly detrimental to vulnerable people, such as senior citizens and infants, understanding how ambient air pollution affects health is very important for making relevant policy [2,3]. Several studies have been conducted regarding the association of PM₁₀ and health in South Korea. One study [3] showed that post-neonatal infants are most susceptible to PM₁₀ in terms of mortality, especially respiratory mortality. Post-neonatal mortality increased by 14.2% for each 42.9-µg/m³ rise in PM₁₀. Another study [4] showed that with a 43.12-µg/m³ increase in PM₁₀ concentration, daily non-accidental, respiratory, cardiovascular, and cerebrovascular mortality increased by 3.7%, 13.9%, 4.4%, and 6.3%, respectively.

Numerous epidemiological studies have investigated the associations of PM₁₀ and PM_2.5 with mortality worldwide [5,6,7,8,9]. One study [5] concluded that PM₁₀ concentration is associated with total, cardiovascular, and respiratory daily mortality in 16 cities in China using a Poisson regression model. Another study [6] showed how air pollutants affect human health in 25 European cities across 12 countries. They determined that a decrease in PM_2.5 level could have led to a gain in life expectancy, postponing many deaths in the city. A European study [7] determined that the daily total mortality rose as PM₁₀ and black smoke concentrations increased with a Poisson regression model. The study [8] investigated the short-term positive effects of PM₁₀ and PM_2.5 on all-cause, cardiorespiratory, and other-cause mortality in the United States using a Bayesian hierarchical model. Another study in the United States [9] showed that increases in PM₁₀ and PM_2.5 concentrations led to increased numbers of post-neonatal infant deaths in 96 cities with populations greater than 250,000. Most of the domestic and foreign studies mentioned above focused on the effects of air pollutants on health outcomes in specific locations. However, to establish national-level air pollution control policy, a study should consider the entire spatial domain of a country.

Spatio-temporal modeling of the association between air pollution and health has recently received attention in environmental epidemiological studies [10,11]; however, there are some challenges in this domain. First, there is frequently a misalignment problem, due to the different dataset sources. This is because air pollution monitoring data are collected at monitoring stations, while health data are collected at an aggregated areal level. There are several approaches to overcome this problem. For example, after spatial modeling of air pollution data, spatial kriging can be conducted, and areal-level air pollution estimates can be produced [12,13,14,15]. Alternatively, one study [16] proposed a Markov Cube Kriging method to improve computational efficiency. In this paper, we applied a commonly used numerical interpolation method to obtain areal-level air pollution data. Second, there is a spatial confounding bias because space-time varying air pollution covariates and space-time random effects are simultaneously included in space-time health modeling [17]. For example, air pollution may be related to space-time random effects, which would bias the estimate on the effects of air pollution. In this paper, to better estimate air pollution effects without the influence of spatial confounding bias, we adopted the two-stage model proposed in [15].

Until now, most of the epidemiologic studies of air pollution conducted in South Korea have focused on how temporally varying air pollutant concentrations affected mortality in individual cities [18,19], but such studies should be conducted over entire areas of Korea to establish national-level air pollution control policy for health improvement. No studies have examined the relationship between temporally varying air pollutants and mortality throughout South Korea. Additionally, there have been almost no studies on the link between mortality and PM_2.5 in South Korea [18]. Herein, we investigated the associations of PM₁₀ and PM_2.5 with mortality throughout South Korea from 2012 to 2015. A two-stage Bayesian spatio-temporal hierarchical model was employed to better estimate the effects of air pollution on mortality outcomes, as well as to better predict the mortality associated with ambient air pollution.

2. Data and Methods 2.1. Study Domain

This study was conducted over the entire area of South Korea and accounted for the period from 2012 to 2015. There are 250 administrative districts in South Korea (Figure 1) and, because some administrative districts changed during the study period, we conducted analyses based on the administrative areas as of 2012. We used the data from 2012–2014 for model fitting, and data from 2015 for forecasting and evaluation.

2.2. Data Description

We obtained monthly mortality data for 250 administrative districts using the Micro Data Integrated Service (MDIS) from Statistics Korea [20]. We used three types of mortality data: total mortality (all deaths except ICD codes V01–Y98), cardiovascular mortality (I00–I99), and respiratory mortality (J00–J99).

We acquired monthly average concentrations of PM₁₀ and PM_2.5 and meteorological data (temperature, humidity, and wind speed) for 250 districts for the years 2012–2015. The district-level air pollution data were produced based on Pun’s interpolation method [21], a data-assimilation method that combines data from monitoring stations and numerical model outputs.

The regional deprivation index (RDI) was also added to the covariates to control for socio-economic status effects. Because there was little temporal variation in the regional deprivation index, we only used data from 2010. Higher index values indicated that the area was more economically deprived. Because deaths occur more often in areas with relatively large elderly populations, we used the number of people over the age of 65 as an offset. Table 1 shows the summary statistics for the 36 months of data (2012–2014) and 250 administrative districts.

2.3. Statistical Analysis

We proposed a two-stage Bayesian hierarchical spatio-temporal model (Model 3) to capture the spatial and temporal dynamics. Two competing models (Model 1, Model 2) were used to compare performances of the proposed model, Model 3. Mathematical expressions for Model 1 and Model 2 are as follows:

1. Model 1: log(_θit)=_β0+_{β1 Xit}+_S1(_Zit)+_S2(_Wit)+_S3(_Qit)+δ_Di ,

2. Model 2: log(_θit)=_β0+_{β1 Xit}+_S1(_Zit)+_S2(_Wit)+_S3(_Qit)+δ_Di+_ui+_vi+_kt+_lt+_ϕit .

The observed mortality for area i and month t , _yit , followed a Poisson distribution with mean _{θit Nit} , where _θit is the relative risk, and _Nit is the elderly population. The constant term _β0 indicates the intercept of the log of the relative risk that was common to all areas and months. Covariates _Xit, _Zit, _Wit , and _Qit are the monthly average air pollutant concentrations of _PM10 or _PM2.5 , temperature, humidity, and wind speed, respectively. The socio-economic covariate _Di indicates the deprivation index. The smoothing function _Si(·), i=1,2,3 denotes a natural cubic spline function to explain the nonlinear effects of meteorological variables on mortality. The degrees of freedom of the natural cubic spline were 3, 2, and 3 for temperature, humidity, and wind speed, respectively. The parameters _β1 and δ denote the regression coefficients _Xit and _Di , respectively. The random effects _ui and _vi are spatially uncorrelated and correlated terms, respectively, and _kt and _lt are temporally uncorrelated and correlated terms. Lastly, the random component _ϕit is the space-time interaction term.

To deal with the spatial confounding bias problem [17] in Model 2, a two-stage model [22] was considered. In the first stage, the Poisson regression model with only covariates (Model 1) was used. Using this model, we acquired the estimated relative risk _θit^ , and the continuous-type residuals _rit^ were obtained from

_rit^=log(_yitNit)−log(_θit^)

To capture the extra spatial and temporal variations in the residuals, we considered the following model:

_rit^|_Sit~N(_Sit,_σr2)

_Sit=_ui+_vi+_kt+_lt +_ϕit

In the second stage, our model was expressed as follows, using _Sit^, the estimated _Sit :

1. Model 3: log(_θit)=_β0+_{β1 Xit}+_S1(_Zit)+_S2(_Wit)+_S3(_Qit)+δ_Di+_Sit^.

Regression coefficient estimation was performed at this stage.

In the Bayesian framework, we used non-informative prior distributions for the parameters. For the intercept _β0 and air pollutant coefficient _β1 , we assumed normal distribution with zero 0 and variance 1,000,000, which is a fairly flat prior. For random effects, spatially and temporally uncorrelated terms had independent and identical normal distributions with zero mean hyperparameters _σu and _σk . The spatially correlated random term _vi had a conditional autoregressive (CAR) [23] prior, and the temporally correlated random term _lt had an autoregressive (AR)(1) prior. For the interaction term _ϕit , we considered four of the types proposed in [24]. The interaction term with different spatial trends for each time unit showed the best performance. Uniform distributions with lower bound 0 and upper bound 10 were specified for all hyperparameters.

The three models described above were fitted for the 2012–2014 data. Based on the results from the fitted models, we forecasted death counts for 2015 for all administrative regions in South Korea. We followed the forecasting scheme used in [25].

Bayesian analyses were carried out using the WinBUGS statistical package [26] Two parallel Monte Carlo Markov Chains (MCMC) were used with different initial values. To assess sample convergence, we utilized trace plots, auto-correlation plots, and the Gelman–Rubin statistic [27]. After burn-in, we generated 2500 samples for each chain with thin 50, resulting in a total of 5000 posterior samples. Including the burn-in period, it took around 30 hours to obtain 5000 posterior samples with a CPU Intel Xeon gold 5118 2.3 GHz and RAM 32 GB computer. The open source software R [28] was used to produce the figures in this paper.

We compared the performance of Models 1–3 to identify the model with the best performance, which is shown in Table 2. For comparison, the deviance information criterion (DIC) and mean squared prediction error (MSPE) were used to evaluate model fitness in the Bayesian framework and prediction performance, respectively. The mathematical expression for MSPE is as follows:

MSPE=1N_{∑i=1I ∑t=1T} ^{(_yit−_yit^)2},

where _yit is the observed value, and _yit^ is the predicted value. DIC is defined as follows:

DIC=pD+D(α)¯,

where D(α)¯ is the posterior mean of deviance D(α)=−2logf(y|α) , and pD is D(α)¯−D(α^) , where α^ is the posterior mean of the parameter α .

To investigate the impact of the degrees of freedom in the spline functions of meteorological variables, we performed a sensitivity analysis by changing degrees of freedom from 2 to 12. The models differed very little in terms of regression coefficient estimates of air pollution and model performance.

3. Results

To show the spatial distributions of air pollutants and mortality in South Korea over the period of 2012–2014, we used the average values of air pollutant concentrations for 36 months and the total death counts for each region over the 36-month period, which are presented in Figure 2. Each map is partitioned into four colors based on quantile. As shown in Figure 2, the total number of deaths from 2012 to 2014 was highest in Nowon-gu, in Seoul, at 7030, and lowest in Gyeryong-si, in Chungcheongnam-do, at 394. Moran’s I statistic was applied to measure the spatial dependencies of mortality and air pollution data at the administrative region level. The null hypothesis of Moran’s I statistic is that there is no spatial correlation among the data. As the adjacency matrix from Moran’s I statistic, various neighborhood structures including binary, row-standardized, and globally standardized techniques were considered. Binary neighborhood structure used a matrix with a value of 1 if it was adjacent, or 0 otherwise. Row-standardized neighborhood structure was used to divide the binary matrix by the sum of the rows. Globally standardized neighborhood structure was used to divide the binary matrix by the mean of the sum of the rows. Moran’s I statistics using row-standardized neighborhood structure of the total, cardiovascular, and respiratory mortality data were 0.27, 0.25, and 0.14, respectively. The results of Moran’s I test for all neighborhood structures showed that all mortality data had statistically significant spatial dependency, with a p-value less than 0.05. Similarly, Moran’s I test indicated spatial variations in _PM10 and _PM2.5 with p-values less than 0.05.

Model performance for total mortality and _PM10 is shown in Table 2. Smaller MSPE and DIC values indicate a better model. While the difference in MSPE between Model 2 and Model 3 was not large, Model 3 yielded a smaller DIC than the other models. Therefore, Model 3 was chosen as the best model. Based on these results, we fitted only Model 3 for cardiovascular and respiratory deaths.

Table 3 gives the posterior summaries of the estimated coefficients of the explanatory variables. The relative risks of air pollutants are based on a 10-μg/m³ increase. The estimated relative risks (RR) of _PM10 and deprivation index for total mortality were 1.011 and 1.001, respectively. Since 95% credible intervals of all variables did not include 1, all estimated coefficients except _PM10 for respiratory mortality were statistically significant. When the relative risk was larger than 1, we inferred that the number of deaths increases as the corresponding variable increases. Figure 3 shows the calibration plots of the observed and estimated numbers of deaths for 2012–2014 for the best model, Model 3. The points are close to the line y=x (red line in Figure 3), indicating that the estimation matched the observations well. Figure 4 shows the calibration plots for the observed and forecasted death counts in 2015 using Model 3 with the 2012–2014 data. It shows that the forecasted numbers of deaths were very similar to the observed numbers of deaths. To see the forecasting performance in detail, we calculated the Pearson residuals of the observed and forecasted values, which are shown in Figure S1 in the supplementary material.

4. Discussion

This study focuses on verifying the association between air pollution and human health over the entire area of South Korea for the years 2012–2015. To adjust for other effects on mortality, we also used meteorological factors and regional deprivation index as covariates.

After controlling for the effects of weather variables, the deprivation index was positively associated with different types of mortality. A previous study [19] found that socioeconomic status, as represented by marital status, education level, and occupation, was linked to different causes of death. Because the deprivation index is a comprehensive index encompassing all types of social or economic statuses for small geographic areas, areas with higher values for this index represent areas with greater deprivation. Deprived areas generally have higher death rates; therefore, our findings are consistent with those of [19].

After controlling for other factors, increases in the concentrations of _PM10 and _PM2.5 were associated with increases in various types of mortality. Other than the non-significant association of _PM10 with respiratory death, these findings are consistent with those from several previous studies [29,30,31]. One previous study [29] determined that the daily _PM10 concentration was positively associated with daily mortality. One of their findings was that _PM10 was positively associated with all-cause, cardiovascular, and respiratory mortalities after adjusting for weather variables [29]. Although there are very few existing studies addressing the relationship between _PM2.5 and death in South Korea, the effects of _PM2.5 on mortality can be inferred from previous studies on _PM10 and human health.

One of our study’s strengths is that it covered the entire area of South Korea from 2012 to 2015. Although some studies have considered different geographic areas in South Korea, most did not evaluate spatial and temporal dynamics. For instance, [18] focused on the seven metropolitan cities in South Korea, but analyzed each geographic area independently. In contrast, our study area encompassed all 250 South Korean administrative areas, and we also considered spatial and temporal variations. We utilized spatial and temporal random effects to capture the spatial and temporal dependency structures not captured by the covariates. As far as we are concerned, this study is the first attempt to illustrate the link between air pollutants and human health using a Bayesian two-stage spatio-temporal model.

Despite the many strengths of this study, there are several limitations. First, the regression coefficient of particulate matter is the same for all administrative areas in our proposed model. However, since the effect of air pollution on death might vary across regions, region-specific regression coefficients may be more suitable. This will be one of the future directions of our research into air pollution and human health. Second, stratification of age and sex might lead to a more detailed understanding of air pollution and health. Several previous studies have conducted age- and sex-specific analyses and obtained meaningful findings [5]. Lastly, the results of this paper may be biased due to the possibility of exposure misclassification resulting from the aggregated nature of the data. For example, we assigned air pollution averages to all cases in each district for a given month, which may not truly capture the monthly average exposure for the 30 days prior to the day of mortality. Ideally, individual level data with precise location and date of death are needed to assess daily time-lagged effects of air pollution exposure [32,33,34]. Future research should be geared to address the above limitations.

5. Conclusions

In conclusion, this study contributes to an understanding of the relationship between air pollution and human health using a Bayesian two-stage spatio-temporal model. Our findings show that air pollution positively affects total, cardiovascular, and respiratory mortalities. Additionally, this study indicates the need for considering spatio-temporal dynamics in epidemiological studies.

Figure 1. Map of administrative districts in South Korea as of 2012.

Figure 2. Spatial distributions of air pollutants and mortality. (a) and (b) are monthly averages of air pollutants from 2012 to 2014; (c), (d), and (e) are monthly sums of deaths from 2012 to 2014.

Figure 3. Calibration plots for the observed and estimated values of 2012-2014 using Model 3.

Figure 4. Calibration plots for the observed and forecasted values of 2015 using Model 3.

Type	Variable	Unit	Mean	SD	Med	Q1	Q3	Min.	Max
Mortality	Total	person	78.8	43.2	73	44	106	2	236
	Cardiovascular	person	19.3	11.4	17	11	26	0	80
	Respiratory	person	7.7	4.8	7	4	10	0	35
Air pollutant	PM10	µg/m3	37.7	12.1	36.8	28.8	46	3.7	83.2
	PM2.5	µg/m3	29.4	10.5	28.3	21.5	36.7	1.6	69.6
Meteorological data	Temperature	°C	13.2	9.6	14.1	5.2	21.7	−7.9	29.7
	Humidity	%	71	9	71.1	64.6	78.1	44.2	92.6
	Wind speed	m/s	2.9	0.9	2.7	2.3	3.3	1.3	8.2
Extra data	RDI		−0.1	8.4	−1.3	−6.8	6.6	−22.5	16.5
	Population	person	202,196	154,313	170,220	60,151	308,111	18,036	669,068

RDI: Regional deprivation index, Med: median, Q1: first quantile, Q3: third quantile.

Model	MSPE	Deviance	pD	DIC
Model 1	184.22	72,690	11	72,701
Model 2	77.09	63,597	548	64,145
Model 3	79.25	63,098	10	63,108

MSPE: Mean squared prediction error, pD: the effective number of parameters, DIC: Deviance information criterion.

Air pollutant	Explanatory Variables	Total	Cardiovascular	Respiratory
		Mean (95% CI)	Mean (95% CI)	Mean (95% CI)
PM10	Air pollutant	1.011 (1.008, 1.013)	1.014 (1.009, 1.019)	0.998 (0.991, 1.006)
	Deprivation index	1.001 (1.001, 1.002)	1.003 (1.002, 1.004)	1.016 (1.015, 1.017)
PM2.5	Air pollutant	1.028 (1.026, 1.031)	1.047 (1.042, 1.053)	1.045 (1.036, 1.053)
	Deprivation index	1.001 (1.001, 1.001)	1.003 (1.002, 1.004)	1.016 (1.015, 1.017)

Supplementary Materials

The following are available online at https://www.mdpi.com/1660-4601/16/12/2111/s1, Figure S1: Plots of Pearson residuals versus the forecasted values for different types of mortality and air pollutants.

Author Contributions

Conceptualization, J.C.; Methodology, D.K. and J.C.; Software, D.K., Y.J., and H.C.; Validation, D.K.; Formal Analysis, D.K., Y.J., and H.C.; Investigation, J.C., S.-s.H., and Y.K.; Resources, S.-s.H. and Y.K.; Data Curation, J.C., S.-s.H., and Y.K.; Writing-Original Draft Preparation, J.C., D.K., Y.J., and H.C.; Writing-Review & Editing, J.C. and D.K.; Visualization, D.K., Y.J., and H.C.; Supervision, J.C.; Project Administration, J.C. and S.-s.H.; Funding Acquisition, J.C. and S.-s.H.

Funding

This research was supported by the National Strategic Project-Fine Particle of the National Research Foundation of Korea (NRF) funded by the Ministry of Science and ICT (MSIT), the Ministry of Environment (ME), and the Ministry of Health and Welfare (MOHW) (NRF-2017M3D8A1092013) and by the Bio & Medical Technology Development Program of the National Research Foundation (NRF) funded by the Ministry of Science & ICT (No. 2017M3A9G8084539).

Conflicts of Interest

The authors declare no conflict of interest.

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

References

1. Samet, J.M.; Dominici, F.; Curriero, F.C.; Coursac, I.; Zeger, S.L. Fine particulate air pollution and mortality in 20 US cities, 1987-1994. N. Engl. J. Med. 2002, 343, 1742-1749.

2. Kwon, H.; Cho, S.; Chun, Y.; Lagarde, F.; Pershagen, G. Effects of the Asian dust events on daily mortality in Seoul, Korea. Environ. Res. 2002, 90, 1-5.

3. Ha, E.; Lee, J.; Kim, H.; Hong, Y.; Lee, B.; Park, H.; Christiani, D. Infant susceptibility of mortality to air pollution in Seoul, South Korea. Pediatrics 2003, 111, 284-290.

4. Kim, H.; Kim, Y.; Hong, Y. The lag-effect pattern in the relationship of particulate air pollution to daily mortality in Seoul, Korea. Int. J. Biometeorol. 2003, 48, 25-30.

5. Chen, R.; Kan, H.; Chen, B.; Huang, W.; Bai, Z.; Song, G.; Pan, G. Association of particulate air pollution with daily mortality: The China Air Pollution and Health Effects Study. Am. J. Epidemiol. 2012, 175, 1173-1181.

6. Pascal, M.; Corso, M.; Chanel, O.; Declercq, C.; Badaloni, C.; Cesaroni, G.; Medina, S. Assessing the public health impacts of urban air pollution in 25 European cities: Results of the Aphekom project. Sci. Total Environ. 2013, 449, 390-400.

7. Katsouyanni, K.; Touloumi, G.; Samoli, E.; Gryparis, A.; Le Tertre, A.; Monopolis, Y.; Rossi, G.; Zmirou, D.; Ballester, F.; Boumghar, A.; et al. Confounding and Effect Modification in the Short-Term Effects of Ambient Particles on Total Mortality: Results from 29 European Cities within the APHEA2 Project. Epidemiology 2001, 12, 521-531.

8. Dominici, F.; Peng, R.D.; Zeger, S.L.; White, R.H.; Samet, J.M. Particulate air pollution and mortality in the United States: Did the risks change from 1987 to 2000? Am. J. Epidemiol. 2007, 166, 880-888.

9. Woodruff, T.; Darrow, L.; Parker, J. Air pollution and postneonatal infant mortality in the United States, 1999-2002. Environ. Health Perspect. 2007, 116, 110-115.

10. Carlin, B.P.; Gelfand, A.E.; Banerjee, S. Hierarchical Modeling and Analysis for Spatial Data; Chapman and Hall CRC Press: Boca Raton, FL, USA, 2014.

11. Cressie, N.; Wikle, C.K. Statistics for Spatio-Temporal Data; John Wiley & Sons: Hoboken, NJ, USA, 2015.

12. Sahu, S.K.; Bakar, K.S. A comparison of Bayesian models for daily ozone concentration levels. Stat. Methodol. 2012, 9, 144-157.

13. Bakar, K.S.; Sahu, S.K. spTimer: Spatio-temporal bayesian modelling using R. J. Stat. Softw. 2015, 63, 1-32.

14. Del, S.S.; Ranalli, M.G.; Bakar, K.S.; Cappelletti, D.; Moroni, B.; Crocchianti, S.; Salvatori, R. Bayesian Spatiotemporal Modeling of Urban Air Pollution Dynamics. Top. Methodol. Appl. Stat. Inference 2016, 95-103.

15. Lawson, A.B.; Choi, J.; Cai, B.; Hossain, M.; Kirby, R.S.; Liu, J. Bayesian 2-stage space-time mixture modeling with spatial misalignment of the exposure in small area health data. J. Agric. Biol. Environ. Stat. 2012, 17, 417-441.

16. Liang, D.; Kumar, N. Time-space Kriging to address the spatiotemporal misalignment in the large datasets. Atmos. Environ. 2013, 72, 60-69.

17. James, S.H.; Brian, J.R. Adding spatially-correlated errors can mess up the fixed effect you love. Am. Stat. 2010, 64, 325-334.

18. Lee, H.; Kim, H.; Honda, Y.; Lim, Y.H.; Yi, S. Effect of Asian dust storms on daily mortality in seven metropolitan cities of Korea. Atmos. Environ. 2013, 79, 510-517.

19. Son, J.Y.; Lee, J.T.; Kim, H.; Yi, O.; Bell, M.L. Susceptibility to air pollution effects on mortality in Seoul, Korea: A case-crossover analysis of individual-level effect modifiers. J. Expo. Sci. Environ. Epidemiol. 2012, 22, 227.

20. Microdata Integrated Service Home Page. Available online: https://mdis.kostat.go.kr (accessed on 13 June 2019).

21. Pun, K.; Seigneur, C. Using CMAQ to Interpolate among CASTNET Measurements; Atmospheric and Environmental Research Inc.: San Ramon, CA, USA, 2006.

22. Choi, J.; Lawson, A.B. A Bayesian two-stage spatially dependent variable selection model for space-time health data. Stat. Methods Med. Res. 2018, 1-13.

23. Besag, J. Spatial interaction and the statistical analysis of lattice systems. J. R. Stat. Soc. Ser. B Stat. Methodol. 1974, 36, 192-236.

24. Knorr-Held, L. Bayesian modelling of inseparable space-time variation in disease risk. Stat. Med. 2000, 19, 2555-2567.

25. Bowman, D.D.; Liu, Y.; McMahan, C.S.; Nordone, S.K.; Yabsley, M.J.; Lund, R.B. Forecasting United States heartworm Dirofilaria immitis prevalence in dogs. Parasites Vectors 2016, 9, 540.

26. MRC Biostatistics Unit Home Page. Available online: https://www.mrc-bsu.cam.ac.uk/software/bugs/the-bugs-project-winbugs/ (accessed on 13 June 2019).

27. Gelman, A.; Rubin, D.B. Inference from iterative simulation using multiple sequences. Stat. Sci. 1992, 7, 457-472.

28. The R Project for Statistical Computing Home Page. Available online: https://www.r-project.org (accessed on 13 June 2019).

29. Park, A.; Hong, Y.; Kim, H. Effect of changes in season and temperature on mortality associated with air pollution in Seoul, Korea. J. Epidemiol. Community Health 2011, 65, 368-375.

30. Pope, C.A., III; Burnett, R.T.; Thun, M.J.; Calle, E.E.; Krewski, D.; Ito, K.; Thurston, G.D. Lung cancer, cardiopulmonary mortality, and long-term exposure to fine particulate air pollution. Jama 2002, 287, 1132-1141.

31. Jerrett, M.; Burnett, R.T.; Ma, R.; Pope, C.A., III; Krewski, D.; Newbold, K.B.; Thurston, G.; Shi, Y.; Finkelstein, N.; Calle, E.E.; et al. Spatial analysis of air pollution and mortality in Los Angeles. Epidemiology 2005, 16, 727-736.

32. World Health Organization. Quantification of Health Effects of Exposure to Air Pollution: Report on A WHO Working Group; WHO Regional Office for Europe: Copenhagen, Denmark, 2000.

33. Kan, H.; London, S.J.; Chen, G.; Zhang, Y.; Song, G.; Zhao, N.; Jiang, L.; Chen, B. Season, sex, age, and education as modifiers of the effects of outdoor air pollution on daily mortality in Shanghai, China: The Public Health and Air Pollution in Asia (PAPA) Study. Environ. Health Perspect. 2008, 116, 1183-1188.

34. Burnett, R.T.; Brook, J.; Dann, T.; Delocla, C.; Philips, O.; Cakmak, S.; Vincent, R.; Goldberg, M.S.; Krewski, D. Association between particulate-and gas-phase components of urban air pollution and daily mortality in eight Canadian cities. Inhal. Toxicol. 2000, 12, 15-39.