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

Previous studies have shown an association between mortality and ambient air pollution in South Korea. However, these studies may have been subject to bias, as they lacked adjustment for spatio-temporal structures. This paper addresses this research gap by examining the association between air pollution and cause-specific mortality in South Korea between 2012 and 2015 using a two-stage Bayesian spatio-temporal model. We used 2012–2014 mortality and air pollution data for parameter estimation (i.e., model fitting) and 2015 data for model validation. Our results suggest that the relative risks of total, cardiovascular, and respiratory mortality were 1.028, 1.047, and 1.045, respectively, with every 10-µg/m³ increase in monthly PM_2.5 (fine particulate matter) exposure. These findings warrant protection of populations who experience elevated ambient air pollution exposure to mitigate mortality burden in South Korea.

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. Summary of variables (2012–2014).

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.

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:

• Model 1:$\log \left({\theta }_{it}\right)={\beta }_0+{\beta }_1{X}_{it}+S_1\left(Z_{it}\right)+S_2\left(W_{it}\right)+S_3\left(Q_{it}\right)+\delta D_i$,
• Model 2: $\log \left({\theta }_{it}\right)={\beta }_0+{\beta }_1{X}_{it}+S_1\left(Z_{it}\right)+S_2\left(W_{it}\right)+S_3\left(Q_{it}\right)+\delta D_i+u_i+v_i+k_t+l_t+{\varphi }_{it}$.

The observed mortality for area $i$ and month $t$, $y_{it}$, followed a Poisson distribution with mean ${\theta }_{it}{N}_{it}$, where ${\theta }_{it}$ is the relative risk, and $N_{it}$ is the elderly population. The constant term ${\beta }_0$ indicates the intercept of the log of the relative risk that was common to all areas and months. Covariates $X_{it}, Z_{it}, W_{it}$, and $Q_{it}$ are the monthly average air pollutant concentrations of ${\mathrm{PM}}_{10}$ or ${\mathrm{PM}}_{2.5}$, temperature, humidity, and wind speed, respectively. The socio-economic covariate $D_i$ indicates the deprivation index. The smoothing function $S_i\left(·\right), 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 ${\beta }_1 \mathrm{and} \delta$ denote the regression coefficients $X_{it}$ and $D_i$, respectively. The random effects $u_i$ and $v_i$ are spatially uncorrelated and correlated terms, respectively, and $k_t$ and $l_t$ are temporally uncorrelated and correlated terms. Lastly, the random component ${\varphi }_{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 $\widehat{{\theta }_{it}}$, and the continuous-type residuals $\widehat{r_{it}}$ were obtained from

$$\widehat{r_{it}}=\log \left(\frac{y_{it}}{N_{it}}\right)-\log (\widehat{{\theta }_{it}})$$

(1)

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

$$\widehat{r_{it}}|S_{it}~\mathrm{N}\left(S_{it},{\sigma }_r^2\right)$$

(2)

$$S_{it}=u_i+v_i+k_t+l_t +{\varphi }_{it}$$

(3)

In the second stage, our model was expressed as follows, using $\widehat{S_{it}},$ the estimated $S_{it}$:

• Model 3:$\log \left({\theta }_{it}\right)={\beta }_0+{\beta }_1{X}_{it}+S_1\left(Z_{it}\right)+S_2\left(W_{it}\right)+S_3\left(Q_{it}\right)+\delta D_i+\widehat{S_{it}}.$

Regression coefficient estimation was performed at this stage.

In the Bayesian framework, we used non-informative prior distributions for the parameters. For the intercept ${\beta }_0$ and air pollutant coefficient ${\beta }_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 ${\sigma }_u$ and ${\sigma }_k$. The spatially correlated random term $v_i$ had a conditional autoregressive (CAR) [23] prior, and the temporally correlated random term $l_t$ had an autoregressive (AR)(1) prior. For the interaction term ${\varphi }_{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. Model performance.

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.

. 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=\frac{1}{N}{\sum }_{i=1}^I{\sum }_{t=1}^T{\left(y_{it}-\widehat{y_{it}}\right)}^2,$$

(4)

where $y_{it}$ is the observed value, and $\widehat{y_{it}}$ is the predicted value. DIC is defined as follows:

$$DIC=pD+\overline{D\left(\alpha \right)},$$

(5)

where $\overline{D\left(\alpha \right)}$ is the posterior mean of deviance $D\left(\alpha \right)=-2\log f\left(y|\alpha \right)$, and $pD$ is $\overline{D\left(\alpha \right)}-D\left(\widehat{\alpha }\right)$, where $\widehat{\alpha }$ is the posterior mean of the parameter $\alpha$.

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 ${\mathrm{PM}}_{10}$ and ${\mathrm{PM}}_{2.5}$ with p-values less than 0.05.

Model performance for total mortality and ${\mathrm{PM}}_{10}$ 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. Relative risks with 95% credible intervals for mortality in Model 3.

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)

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 ${\mathrm{PM}}_{10}$ 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 ${\mathrm{PM}}_{10}$ 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 $\mathrm{y}=\mathrm{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 ${\mathrm{PM}}_{10}$ and ${\mathrm{PM}}_{2.5}$ were associated with increases in various types of mortality. Other than the non-significant association of ${\mathrm{PM}}_{10}$ with respiratory death, these findings are consistent with those from several previous studies [29,30,31]. One previous study [29] determined that the daily ${\mathrm{PM}}_{10}$ concentration was positively associated with daily mortality. One of their findings was that ${\mathrm{PM}}_{10}$ 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 ${\mathrm{PM}}_{2.5}$ and death in South Korea, the effects of ${\mathrm{PM}}_{2.5}$ on mortality can be inferred from previous studies on ${\mathrm{PM}}_{10}$ 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.