Asian Dust Storm Elevates Children’s Respiratory Health Risks: A Spatiotemporal Analysis of Children’s Clinic Visits across Taipei (Taiwan)

Concerns have been raised about the adverse impact of Asian dust storms (ADS) on human health; however, few studies have examined the effect of these events on children’s health. Using databases from the Taiwan National Health Insurance and Taiwan Environmental Protection Agency, this study investigates the documented daily visits of children to respiratory clinics during and after ADS that occurred from 1997 to 2007 among 12 districts across Taipei City by applying a Bayesian structural additive regressive model controlled for spatial and temporal patterns. This study finds that the significantly impact of elevated children’s respiratory clinic visits happened after ADS. Five of the seven lagged days had increasing percentages of relative rate, which was consecutively elevated from a 2-day to a 5-day lag by 0.63%∼2.19% for preschool children (i.e., 0∼6 years of age) and 0.72%∼3.17% for school children (i.e., 7∼14 years of age). The spatial pattern of clinic visits indicated that geographical heterogeneity was possibly associated with the clinic’s location and accessibility. Moreover, day-of-week effects were elevated on Monday, Friday, and Saturday. We concluded that ADS may significantly increase the risks of respiratory diseases consecutively in the week after exposure, especially in school children.


Introduction
Every winter and spring season, Asian dust storms (ADS) are frequently a major health problem in many cities throughout the world. Annual dust emissions total approximately 43 million tons, with springtime dust emissions accounting for almost half of the annual totals [1]. ADS usually originate in the deserts of Mongolia and China, lift soil particles into the atmosphere, and carry them east and southeast towards Korea, Japan, Taiwan, and even the Philippines [2,3]. In addition to highly-elevated ambient particulate matter (PM) concentration, the PM composition during ADS differs from that during non-ADS periods. Reports from central Taiwan show that the concentration of the crustal elements and sea salt species in PM 2.5-10 (i.e., particulate matters with granulometric diameter between 2.5 and 10 mm) during ADS exceeds the mean concentration of these same entities during non-ADS periods by at least 2 units of factors [4]. Moreover, during ADS in China, higher PM concentrations were measured in Korea and Japan, and the ADS particles were also observed to be rich in some elements that include aluminum, iron, and calcium [5][6][7].
Several environmental epidemiologic studies have demonstrated evidence of an adverse relationship between airborne particles and human health. Airborne particles uniquely impact the mortality and hospitalization rates of several diseases, most notably respiratory diseases. Most research has concluded that the effects of PM, especially composed of fine particles, are more severe than that of other air pollutants because they can be inhaled deeply into the lungs [8]. A European study revealed that each 10 mg/m 3 increase of PM 10 (i.e., particulate matters with granulometric diameter less than 10 mm) raises the relative rate of clinic visit by 1.20% for asthma and 1.00% for chronic obstructive pulmonary disease [9]. In addition, the effect of fine particulate air pollution on mortality may contribute to higher instances of stroke and respiratory deaths [10]. Besides PM, the adverse health impact of CO, O 3 , NO 2 and SO 2 has also been well-documented historically [11][12][13].
The health impact of ADS has been discussed by comparing mortality rates or emergency room admissions between ADS and non-ADS events [14,15]. In addition, researchers have investigated health related lag effects associated with ADS, and they have discovered deferred influence on hospital admissions [16][17][18]. Yang et al. [19] found a statistically significant association between ADS events and primary intracerebral hemorrhagic stroke admissions 3 days after the dust storms. ADS also increased the risk of both asthma and cerebrovascular admissions from a 1-day lag to a 3-day lag [20,21]. Among the research population, children are particularly sensitive to airborne exposure [22]. However, few studies have focused on children's health from the perspective of ADS events [23,24]. Most of the previous studies regarding the health impact analysis of ADS utilize relatively few health observations, and therefore, inferences from these studies may be limited and conservative [25]. Chien et al. [24] showed the elevated rate of children's respiratory clinic usage during one week following ADS; however, no temporal lag effect structure of health impact was discussed.
Geographic heterogeneity has been a salient factor in ambient pollutant distributions [26,27] as well as its associations with health outcomes [12,28]. Nevertheless, few studies have assessed the spatial variation of an ADS's impact on human health. In order to address this issue, this study applies a unidirectional approach [29] under a spatiotemporal model framework to diagnose the spacetime disparity of children's respiratory clinic visits. The influence of ADS on children's health was examined by considering temporal lag effects starting from the end of each ADS event as well as the spatial variation over study areas. This study specifically investigates the daily clinic visits of children with respiratory diseases in 12 districts in Taipei City from 1997 to 2007.

Children's Clinic Data
Initiated in March 1995, Taiwan's National Health Insurance (NHI) program contacted more than 97% of hospitals and clinics nationwide within its first year of inception, enrolling more than 96% of Taiwanese residents. The Taiwan National Health Research Institute maintains the NHI program database, and has established a standard procedure that assures the quality and accuracy of claims data [30]. The NHI database includes ambulatory care expenditures by visit as well as the registries of contracted medical facilities nationwide. The procedure and diagnostic codes are used to retrieve cause-specific data according to diagnosis-related groups or International Classification of Diseases, Ninth Revision, Clinical Modification (ICD-9-CM) classification codes by the Bureau of National Health Insurance. Regarding personal privacy and confidentiality, all individually identifiable health information has been encrypted prior to release, e.g., personal identification or hospital identification numbers.

Dust Storm Data
ADS often occur in northern and northwestern China, impacting Taiwan only under certain atmospheric circumstances. Before 2000, the Department of Atmospheric Science at Chinese Culture University (CCU) was tasked with the responsibility of characterizing, defining, and monitoring ADS events in Taiwan. They came up with the following criteria for identifying an ADS event: 1) dust storm events with PM 10 concentrations .100 mg/ m 3 observed by any air quality monitoring stations located in Wanli, Guanyin, Danshui, and Yilan, and 2) dust storm events with visibility less than 1 km for 24 hours in any of three neighboring First Global GARP Experiment-type ground stations [31]. After 2000, the Taiwan Environmental Protection Agency (TWEPA) became the official organization to define, monitor, and predict ADS. They utilize three distinct steps in order to categorize storm events in Taiwan as ADS and to predict their arrival time. First, the Weather Integration and Nowcasting System is consulted to confirm the occurrence of ADS in Mongolia and China. Second, the Moderate Resolution Imaging Spectroradiometer remote data and several models of ADS are used to track the transport of ADS and to predict the probability of the arrival of ADS in Taiwan. Third, if the data confirm that ADS may blow to Taiwan, the TWEPA will issue an early warning with the estimated arrival date [32]. This study highlights and analyzes 76 storm events that were thus categorized as ADS during the period of 1997-2007 (see Table 1). These ADS events span a total period of 172 dust storm days. Proportionally, these dust storm days account for 4.28% of the entire study period that lasted for a total time span of 4017 days. In addition to the ADS data, ambient pollutants concentrations and temperature have been regularly monitored at the TWEPA stations across Taiwan since 1994. The temperature measurements used in this analysis are the daily observations at the Jhongshan air quality monitoring station located in the most populated area of Taipei City (see Figure 1).

Study Area
Taipei City, located in northern Taiwan, is the country's capital and largest metropolitan area with a population of more than six million inhabitants. Topographically, it is the second largest basin in Taiwan, bounded by the Yangming Mountains to the north, the Linkou mesa to the west, and Snow Mountains to the southeast. The basin region along the rivers is more populated than regions near the mountains, creating serious air pollution problems due to its heavy traffic. In the past decade, Taipei's metropolitan rapid transit system (MRT) has effectively reduced cross-district traveling times within Taipei, and therefore, it has increased the accessibility to major medical facilities. Figure 1 displays the Taipei City map and essential urban information relevant to this study, including geographic topography, districts, major medical facilities, and the MRT. The clinic visit data are aggregated with respect to the 12 districts across Taipei City.

Spatiotemporal Modeling
In this study, assuming Y st is the number of daily children's respiratory clinic visits at calendar time t[(1,2,:::,4017) in district s[(1,2,:::,12), this outcome variable would follow a Poisson distribution by Y st Dm ste POI(m st ) with the expected value E(Y st )~m st and variance V(Y st )~Q s m st , where Q s is an overdispersion parameter representing the variation of clinic visits unable to be calculated by statistical models. A vector of dust storm lag index (DSLI) contains eight dummy variables to represent 0lag to 7-lag days of ADS. Moreover, a vector of day-of-week (DOW) dummy variables from Monday to Saturday is used to control short-term temporal autoregressive correlations. The Bspline with a second order random walk prior is used in a time smoother f (Time) for calendar time to consider long-term autoregressive correlations and in a temperature smoother f (TP) to control nonlinear weather confounding effect. Therefore, a model framework can be constructed by applying a Bayesian structural additive regression (STAR) modeling approach: where b 0 is an intercept for interpreting the overall association for all districts. The parameters b 1 and b 2 are a 166 vector with six coefficients of day-of-week variables and a 168 vector with eight coefficients of DSLI variables. The offset is the logarithm of the district-level population based on the 2000 Census.
In particular, a spatial function f spat (s) was appended to consider potential spatial autocorrelations among the 12 districts. It is actually a Markov random field [33] achieved by a conditional autoregressive prior to following a normal distribution f spat (s)D f spat (s'),s=s e N( P s'[Vs f spat (s')=N s ,s 2 s =N s ), where N s is the number of adjacent districts s' connected to district s, and s' [ V s represents that the district s' belongs to the set of the neighboring districts V s of district s. Note that the geographic information system in this study is boundary data, and that the definition of two connected districts means that they share parts of boundaries.
In the spatial function, the unknown variance, s 2 s , and smoothing parameter are assumed to follow an inverse Gamma distribution IG(0.001, 0.001). The spatial effect can be explained by the relative rate (RR) in district s compared with the mean value for the whole population controlling for spatial autocorrelations [34]. Maps of the spatial function at the district level were presented to visualize the geographic distribution of RR. Spatial effects were also classified into three groups according to their posterior probabilities with respect to the number 1 in the following manner: (i) 80% of the posterior distribution below 1 representing a significantly lower RR than the mean level for the Taipei City area; (ii) 80% of the posterior distribution above 1 representing a significantly higher RR than the mean level for the Taipei City area; (iii) the other districts representing a non-significant difference of RR compared to the mean level for the Taipei City area.
All unknown parameters in this STAR model were estimated by an empirical Bayesian approach using a restricted maximum likelihood (REML) estimation technique [35]. The REML method is a substitution of the maximum likelihood method that avoids the loss of degrees of freedom while estimating fixed effects, and prevents estimators toward zero with biases. Demographic analyses were generated using SAS v9.2 software [36], and the model-based estimation procedure was accomplished in the software package BayesX v2.01 [37]. A p-value ,0.05 was considered statistically significant in estimated coefficients. Table 2  After a 1-day lag, clinic visits decreased in most areas, however this falling trend was inconsistent. Some areas had a higher average during the 5-day lag or 7-day lag. Moreover, according to the daily records at the Jhongshan air quality monitoring station, the average concentration of PM 10 was 90.64 mg/m 3 during dust storm days, which was significantly higher than the average concentration during non-dust storm days (p-value ,0.0001), as shown in Table 3. The average concentration of O 3 during dust storm days was 5.61 ppb higher than that during non-dust storm days significantly (p-value ,0.0001). This finding confirms that during ADS events as verified by the CCU and TWEPA, higher concentrations of two important air pollutants (i.e., PM 10 and O 3 ) were measured. The association between ADS and children's respiratory clinic visits was not positive until the second day after ADS, suggesting the percentage increase of RR for preschool children was 22.53% (p-value ,0.0001; 95% CI = 22.69, 22.36), while it was much lower in school children by 26.28% (p-value ,0.0001; 95% CI = 26.50, 26.06). The negative association lasted through the 1-day lag, and became positive from the 2-day lag to the 7-day lag, except for the 6-day lag. Among lag days with positive associations, preschool children had the highest 2.19% (p-value ,0.0001; 95% CI = 1.95, 2.43) at the 3-day lag; meanwhile for school children, RR at the 7-day lag reached its highest percentage increase by 3.20% (p-value ,0.0001, 95% CI = 2.81, 3.60). Regardless of age stratification, the strongest association happened at the 3-day lag with a 2.40% (p-value ,0.0001; 95% CI = 2.20, 2.59) increase in RR for all children. Figure 1 depicts the distribution of spatial effects attributed to children's respiratory clinic visits in Taipei City. In most districts with preschool children or school children, positive risk of increased clinic visits was prevalent. The range of spatial effect in preschool children was (20.37, 0.27), and it was wider than that for school children (20.22, 0.13). Out of all the 12 districts studied, the Jhongshan District displayed the strongest spatial effect contributed to respiratory clinic visits for both preschool and school children. The positive spatial effect in preschool children's clinic visits was uniformly distributed in the following districts: the Beitu, Shihlin, Neihu, Nangang, and Wanshan District; however, no specific pattern described the spatial heterogeneity in school children's clinic visits. The maps of 80% posterior probability show that 6 of 12 districts demonstrated a significantly positive spatial effect in both preschool children and school children. Combing two groups, 7 of 12 districts had a significantly positive spatial effect, which locations are identical to the finding in preschool children.

Discussion
Due to high PM concentrations and unusual PM compositions during ADS, ADS and their occurrences have been considered to pose a high risk to human health. Recent studies have demonstrated potential health risks associated with ADS in terms of higher mortality rates and hospital admissions [21,38]. Some studies have evaluated the biological plausibility of the ADS to induce adverse health effects. These studies have shown that the particles in ADS can exert toxicological effects on the respiratory system, such as causing pulmonary inflammation and inducing cytotoxicity in rat alveolar cells [39][40][41][42][43][44]. However, other epidemiological studies have noted that the relationship between ADS and adverse health effects, particularly respiratory diseases, is at best uncertain or statistically insignificant [3,21,[45][46][47][48]. One plausible explanation of these inconsistencies may be that the health assessment measures used to evaluate the health impact of the ADS did not adequately capture the health effects. For instance, some of the health measures utilized only accounted for severe cases that required inpatient care, and consequentially, negative environmental events such as ADS may not necessarily induce such severe medical conditions. Furthermore, these previous analyses were based upon observations from limited hospitals [21,25,45,48], and reflected severe health conditions that were exhibited by the most vulnerable individuals within a general population. In contrast, this study provides complete ambulatory and emergency service utilization information from the NHI, and hopefully, better captures the potential health impact of ADS events on the general population.
In recent decades, the importance of spatiotemporal analysis has been emphasized in environmental epidemiological research, especially in quantifying uncertainties in space-time health and exposure data as well as capturing the resulting impact on the estimates of these associations [28,49,50]. Although temporal health impacts caused by ADS have been extensively investigated,  previous studies have seldom considered the spatial heterogeneity of such health data in which the geographic disparity may be derived from geographical variations of medical resources and exposure levels.
In response, this study implemented the STAR modeling approach for the spatiotemporal analysis of children's clinic visits related to ADS. The model identifies temporal patterns of a time process by accounting for linear and nonlinear explanatory variables similar to many time series models, such as the generalized additive model [51]. Moreover, the STAR model reveals the spatial heterogeneity independent of temporal variations by using Markov random fields. In addition, the Bayesian framework of the STAR model allows feasibility to account for the parameter's uncertainty. This novel approach provides a more comprehensive perspective on the impact of ADS on children's clinic visits for respiratory illnesses.
Previous studies have demonstrated different results in the temporal lag effects of adverse human health related ADS. For instance, hospital admissions were prominent 2 days after ADS in one asthma study [45]. Also, in another study, a positive influence was noted between ADS and ischemic stroke hospital admissions on the third day following a dust storm event [19]. At the 1-day lag, the relative risk of the association between ADS and cardiovascular disease hospital admissions is also increased [14]. However, these associations were statistically insignificant. In contrast, hospital admissions records in Taipei City noted a significant increase in ischaemic heart disease admissions at the 2day lag and asthma admissions at the 3-day lag [21]. Meanwhile, total respiratory diseases at the 3-day lag and upper respiratory tract infection in males at the 4-day lag were significant in Minqin City, China [48]. As noted, these findings were mostly based upon the analysis of more severe health measures. This study conducted a population-based study and found that children's respiratory health can be affected by ADS. This impact significantly occurred during most days within a week after a dust storm event. The elevated rates for children's respiratory clinic visits after a dust storm began at the 2-day lag and attained its highest impact at the 3-day lag for preschool children and the 7-day lag for school children.
Several of the following considerations may provide plausible explanations for such an increased rate of children's respiratory clinic visits: First, the impact of the ADS may not necessarily incite immediate respiratory illness or illness severe enough to necessitate the patient to seek medical services. A latency period may exist between the adverse environmental influence and the onset of respiratory symptoms requiring the need for medical services. Second, the adverse weather conditions that exist during ADS, such as strong winds and low visibility [52], often prevent citizens from going out. Third, the increasing popularity of ADS forecasting by the media and governmental agency may increase the population awareness of ADS and their potential health effects. Fourth, since over-the-counter pharmaceuticals are easily accessible and inexpensive in Taiwan, many Taiwanese residents may prefer initiating treatment of their symptoms and their children's symptoms with these products before seeking medical treatment at a clinic, especially if they perceive that their condition is not serious. However, if unsuccessful, treatment with these over-thecounter medications could also account for the lag time noticed in the children's respiratory clinic visits following ADS. Table 4 notes that the consecutive elevated risks may only apply to the children because of their vulnerability to ambient pollutants. Further studies should assess the ADS health impact on other age groups. Table 4 also shows that school children were affected by ADS much easier than preschool children. This may be explained by the fact that Taiwanese schools were not suspended during ADS, and school children may have experienced higher exposures and elevated concentrations of heavy metals and ambient PM compositions than their preschool counterparts who more likely stayed at home [4,53,54]. This type of exposure has been closely associated with the reduction of pulmonary functions in children [23]. Further studies should investigate the temporal fluctuation of daily relative risks that may result from the cross-infection of respiratory diseases among children or the influence payment policy of the Taiwanese NHI system. Interestingly, the day-of-week has been considered as a meaningful confounding factor for clinic visits in Taiwan. It is important to know that local ambulatory service is commonly rendered on a ''first-come, first-serve'' basis, and that physician appointments are not necessary for the regular weekday schedule. Access to any level of healthcare facility or provider is therefore unconstrained during the week. However, weekend medical services must be justified by severe symptoms and result in higher co-payment requirements (out-of-pocket amount) from the NHI. Thus, the service schedule and payment system are important factors affecting the timing of medical-care-seeking behavior. Thus, the temporal pattern of clinic visits is closely associated with this weekend effect.
In Taiwan, the majority of the medical services in hospitals and clinics are closed from Saturday afternoon until Monday morning. Therefore, there is a strong incentive to visit clinics on Fridays and Saturday (before weekend effect), and on Monday (after weekend effect). The day-of-week clinic visit pattern observed in this study is quite consistent with medical care-seeking pattern that has resulted under the current national health care delivery system. Moreover, the highly elevated RR on Monday essentially comprised those patients seeking medical treatment from Saturday to Monday, especially in the case of children who are incapable of accessing clinic care independently and are dependent on a parent's working schedule which allows for only nighttime availability [55].
The spatial heterogeneity of clinical visits also reflects children's respiratory clinic visits (see Figure 2). Compared to Figure 1, the districts with elevated rates were closely linked to the areas with multiple urban medical centers, especially those along the MRT lines, implying that a high usage of ambulatory and emergency services for children's health might logically occur in these districts. Because the NHI program is characterized by its low co-payments and open access to providers without choice restrictions, it encourages those insured under the current Taiwan NHI system to seek care in these medical centers with minimal personal financial impact. Consequently, each person in Taiwan averages 14.2 clinic visits per year. In addition, some people may seek treatment of common diseases in hospitals or even tertiary medical centers, rather than clinics [24]. Further study is required to investigate the relationship between the identified spatial heterogeneity and the locations of medical centers in the study area.

Conclusion
In summary, the spatiotemporal analysis presented in this study identifies the temporal pattern of the health risks during and after ADS and analyzes them day-by-day by considering the spatial confounding factor. The study results clearly show significant and increased rates for respiratory clinic visits in the studied population of children over time in 5 of 7 days after ADS. The findings of this population-based study can provide governmental agencies with an important reference source in order to plan and implement policies that can help to both protect children from the possible adverse health effects of ADS and to provide care related to such health effects.

Author Contributions
Conceived and designed the experiments: HLY LCC. Analyzed the data: HLY LCC. Contributed reagents/materials/analysis tools: HLY LCC CHY. Wrote the paper: HLY LCC CHY. Provided empirical data: CHY HLY.