A Two-Stage Method to Estimate the Contribution of Road Traffic to PM2.5 Concentrations in Beijing, China

Background: Fine particulate matters with aerodynamic diameters smaller than 2.5 micrometers (PM2.5) have been a critical environmental problem in China due to the rapid road vehicle growth in recent years. To date, most methods available to estimate traffic contributions to ambient PM2.5 concentration are often hampered by the need for collecting data on traffic volume, vehicle type and emission profile. Objective: To develop a simplified and indirect method to estimate the contribution of traffic to PM2.5 concentration in Beijing, China. Methods: Hourly PM2.5 concentration data, daily meteorological data and geographic information were collected at 35 air quality monitoring (AQM) stations in Beijing between 2013 and 2014. Based on the PM2.5 concentrations of different AQM station types, a two-stage method comprising a dispersion model and generalized additive mixed model (GAMM) was developed to estimate separately the traffic and non-traffic contributions to daily PM2.5 concentration. The geographical trend of PM2.5 concentrations was investigated using generalized linear mixed model. The temporal trend of PM2.5 and non-linear relationship between PM2.5 and meteorological conditions were assessed using GAMM. Results: The medians of daily PM2.5 concentrations during 2013–2014 at 35 AQM stations in Beijing ranged from 40 to 92 μg/m3. There was a significant increasing trend of PM2.5 concentration from north to south. The contributions of road traffic to daily PM2.5 concentrations ranged from 17.2% to 37.3% with an average 30%. The greatest contribution was found at AQM stations near busy roads. On average, the contribution of road traffic at urban stations was 14% higher than that at rural stations. Conclusions: Traffic emissions account for a substantial share of daily total PM2.5 concentrations in Beijing. Our two-stage method is a useful and convenient tool in ecological and epidemiological studies to estimate the traffic contribution to PM2.5 concentrations when there is limited information on vehicle number and types and emission profile.


Introduction
According to a growing body of epidemiological evidence traffic-related air pollution has been shown to have adverse health impacts. Fine particulate matters with aerodynamic diameters smaller than 2.5 micrometers (PM 2.5 ) pose great public health hazards, including higher risks of respiratory diseases, impaired lung function, asthma attacks, cardiovascular diseases, and potentially also premature death [1].
The particulates generated from combustion are more harmful than those generated from other processes, and traffic emissions are estimated to account for up to 50% of combustion-generated particulates in urban areas in developing countries [2]. According to the Ministry of Environmental Protection of China, traffic emissions have become the main source of air pollution in Beijing [3]. Among all air pollutants, PM 2.5 is of special importance in China due to the rapidly growing number of road vehicles in recent years. By collecting and analyzing aerosol samples of PM 2.5 and PM 10 both in summer and winter seasons at different traffic, industrial and residential areas in Beijing, a multisite study found that industrial and motor vehicle emissions, together with coal burning, were the major contributors to the air-borne particulate pollution in Beijing [4].
Although the Beijing Environmental Protection Bureau started monitoring air pollution in 1984, monitoring of PM 2.5 only started in 2006. Prior to that, PM 2.5 was mainly used for air pollution research purposes [5]. As a result of increasing demand from the public, since October 2012, Beijing has increased its number of fixed air quality monitoring (AQM) stations from 27 to 35 across the entire municipal area. In addition to carbon dioxide, sulfur dioxide, nitrogen dioxide, ozone and PM 10 , PM 2.5 has also been included in the air quality evaluations of these AQM stations. A study found that, while burning of coal for power plants is a major source of air pollution across China, vehicle emissions are one of the biggest sources of PM 2.5 in Beijing, with greater impact than soil dust, fossil fuel combustion, biomass burning and some industrial sources [6]. Although previous studies have clearly shown that the contribution of traffic emissions to total air pollution varies largely with time and space, they were unable to characterize the spatiotemporal features of the traffic-related PM 2.5 because of limited information on location and time period for air sample collection [5].
Chemical mass balanced receptor models and source-oriented chemical transport models have been used to estimate the contributions of various sources to PM 2.5 , but most of them require the knowledge of the chemical profile of both the emissions of the sources and the air samples of the receptors (i.e., the impacted locations) [7,8]. Although other models such as principal component and factor analyses do not require a priori knowledge of the source profile, application of these models yielded controversial results. For example, the estimated motor-vehicle contribution to PM 2.5 ranged from 6% in Beijing, China to 53% in Barcelona, Spain [9].
Although traffic emission is the principal source of intra-urban concentration of PM 2.5 , one reason that the direct measurement of motor-vehicle emission may not be feasible in epidemiological studies is that it is usually not possible to track all the vehicles and measure corresponding components of the traffic-pollutant mix in the whole study area [10]. As a result, different surrogates of traffic-related pollution have been used to assess the contribution of road traffic to ambient air pollution. In epidemiological studies, the commonly used surrogate models include geostatistical interpolation [11], land-use regression [12], dispersion [13] and hybrid [14] models. Hybrid models combine personal activity of residents in the study area and exposure data, and incorporate various measurements, therefore better quantify the contribution of traffic on air pollution, against a background concentration of specific regions. However, none of the models has an ideal surrogate to access the emissions from all sources over time and space, posing a significant challenge in disentangling the contribution of road traffic from other sources.
To improve the assessment of traffic-related contributions to PM 2.5 , a promising method is the deployment of a large number of AQM stations in places where concentrations of PM 2.5 are expected to be highly variable, and with available information on temporal and spatial factors [15]. The intensive air quality data that we collected from 35 AQM stations in Beijing, one of the most populous cities in the world, between 2013 and 2014, provided us a unique opportunity to achieve this purpose. In our paper, we presented a two-stage method using dispersion model and generalized additive mixed model (GAMM) to estimate the contribution of road traffic to PM 2.5 concentrations in Beijing. We used different types of the AQM stations (described in Material and Methods section) to distinguish the emission sources of PM 2.5 , adjusted for the location of these stations, traffic density and meteorological conditions. In the first stage, a Community Multi-scale Air Quality (CMAQ) based model was built to estimate the contribution of road vehicle emission to PM 2.5 as a result of dispersion and decay in the areas represented by background stations [16]. In the second stage, a GAMM with meteorological and geographic data was developed to estimate the non-traffic contribution to PM 2.5 at the rest stations. The traffic contribution to PM 2.5 was then calculated by subtracting the total PM 2.5 concentration with non-traffic concentration. The study was approved by the Institutional Review Board of Karolinska Institutet, Sweden.

Data Collection
Hourly concentrations of PM 2.5 were collected from 35 AQM stations in Beijing from 1 January 2013 to 31 December 2014. The AQM stations were installed by the Beijing Municipal Environmental Protection Bureau. The aim of these stations was to assess the air quality under different conditions from the most polluted area with high density of traffic to the least polluted rural areas in Beijing. Thus, air pollution concentrations of these stations vary largely from each other due to the variation of their distances to pollution sources, e.g., traffic emissions and industrial emissions. The distribution of the AQM stations was shown in Figure 1. These stations scattered from the very south to the north of Beijing, from the central urban areas to countryside, covering most of the spatial regions and typical land types. Geographic information of these stations was attained from College of Resources and Environment, University of Chinese Academy of Sciences. According to the Ambient Air Quality Standards and Technical Regulation on Ambient Quality Index of China, 24-hour concentrations of PM 2.5 and individual air quality index (IAQI) were reported hourly from these stations [17]. The air quality has been classified by Chinese Environmental Protection Agency into six categories, i.e., "Good", "Moderate", "Unhealthy for Sensitive Groups", "Unhealthy", "Very Unhealthy" and "Hazardous" [17,18]. Duplicated records were first removed from the dataset, and the records with empty or 0 value were treated as missing. The missing rate was 9% and no apparent trend was found for the missing values. In total, 553,877 PM 2.5 concentration records were collected from the 35 monitoring stations in 730 days in 2013 and 2014. Values greater than 10 times the 75% percentile or smaller than one-tenth of the 25% percentile of all the records were treated as abnormal values and only included in sensitivity analyses.
Daily meteorological data were collected from National Meteorological Information Center of China in the same period, including air temperature, atmospheric pressure, wind speed, wind direction, volume of rainfall and hours of daylight. Five-minute traffic volume and speed data per 30 minutes for four days from eight crossroads in core districts in Beijing were collected by the College of Resources and Environment, University of Chinese Academy of Sciences. The traffic density of the monitoring stations in these districts was characterized by an inverse function of mean road vehicle speed on the main roads [19].

Fitting Spatial Trend of PM2.5 Concentration
Historical data and previous findings showed that air pollution was often heavier in the southern part than the northern part of Beijing [20], therefore a three-level generalized linear mixed model (GLMM) was fitted between the geographical Y coordinates (i.e., distance from an AQM station to the southern boundary of Beijing ) of the AQM stations in a rectangular coordinate system and the log transformed PM2.5 concentrations (logPM2.5). The Y coordinates were used as an independent variable, whereas calendar days and hours of each calendar day were include as random effects in the model. Because background stations are less but still affected by traffic pollution, and non-traffic portion of PM2.5 pollutants is more geographically stable, fitting a regional non-traffic trend in the study area that takes advantage of the background stations is plausible. The final traffic contribution could be calculated by subtracting the non-traffic portion from the total observed concentration. The two-stage method is described in detail below:

Fitting Spatial Trend of PM 2.5 Concentration
Historical data and previous findings showed that air pollution was often heavier in the southern part than the northern part of Beijing [20], therefore a three-level generalized linear mixed model (GLMM) was fitted between the geographical Y coordinates (i.e., distance from an AQM station to the southern boundary of Beijing ) of the AQM stations in a rectangular coordinate system and the log transformed PM 2.5 concentrations (logPM 2.5 ). The Y coordinates were used as an independent variable, whereas calendar days and hours of each calendar day were include as random effects in the model. Because background stations are less but still affected by traffic pollution, and non-traffic portion of PM 2.5 pollutants is more geographically stable, fitting a regional non-traffic trend in the study area that takes advantage of the background stations is plausible. The final traffic contribution could be calculated by subtracting the non-traffic portion from the total observed concentration. The two-stage method is described in detail below: Badaling, located far away from both urban areas and industrial areas and had few direct traffic and industrial emissions. The air pollution at these stations is mainly from dispersed pollutants, and the PM 2.5 concentration of these stations can be regarded as the background pollution concentration in each region. The five traffic stations include Dongsihuan, Nansanhuan, Qianmen, Xizhimenbei and Yongdingmen, which are less than 10 meters away from the main roads of Beijing, where the PM 2.5 concentration mainly derives from traffic emissions. The two industrial stations include Liulihe and Yufa which are located at the southern boundary between Beijing and Heibei Province where the PM 2.5 concentration is mainly caused by local industrial emission and dispersion. For traffic stations, the PM 2.5 pollutants were mainly from vehicle emissions. The total PM 2.5 concentration of the five stations was considered as a surrogate of the PM 2.5 from traffic emissions. Two industrial stations close to the southern boundary of Beijing are located near to an industrial area of Hebei Province. The PM 2.5 concentrations of these two stations are therefore treated as a surrogate of the industrial emissions.
Based on the Hybrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT) model, backward trajectories were used to track the transport corridors that are regarded as a "region of influence" i.e., the five traffic stations and two industrial stations in our study [21]. Both traffic and industrial stations were considered as PM 2.5 sources of the six background stations because no other major PM 2.5 contributors were found near the background stations. Because dispersion processes are largely additive, PM 2.5 pollution at every station is supposed to be consisted of remaining from daily deposition and dispersion from different emission source points, i.e., traffic, factories and other sources (such as household cooking and coal burning) [22,23]. For the background stations, PM 2.5 contribution other than traffic and industrial dispersed is considered as station-specific background PM 2.5 concentration.
Since distance also played an important role for pollutant dispersion, the inversed value of distance from source stations (i.e., traffic stations and industrial stations) to the receptor stations (i.e., background station) was put as a weight of dispersion factor.
According to the Community Multiscale Air Quality (CMAQ) model, all emissions are assumed to be instantaneously well-mixed and have their own atmospheric lifetime [24]. Therefore a strong daily dependence is expected on consecutive days. We assume that in the condition of wind, PM 2.5 can partly linger for at least one day [25]. Analogously, we built a dispersion model as shown in model (1), in which the PM 2.5 concentration was presented as a summation of traffic dispersion, industrial dispersion and the remaining from daily deposition. Because the pollution carried by wind had a strong positive relationship with the wind flux, a power function was used to fit dispersion effect. The daily deposition of pollution with interaction of wind was fitted by the exponential function: (1) In model (1):Ĉ pptq denotes the expected PM 2.5 concentration at station p on day t. C ppt´1q denotes the observed PM 2.5 concentration on day t-1; D ind p represents the average distance from station p to industrial stations; C indptq denotes the observed PM 2.5 concentration of industrials stations on day t; and D tra f f ic p represents the average distance from station p to traffic stations; C tra f f icptq denotes the observed PM 2.5 concentration of traffic stations on day t;Ŵ indptq denotes the summation of valid flux of wind from industrial stations andŴ tra f f icptq means the summation of valid flux of wind from traffic stations on day t; W avg is the average wind speed of the year; W ptq is the maximum wind speed on day t; and k 1 ,¨¨¨, k 5 are the parameters to be estimated by Levenberg-Marquardt and global minimum algorithm till their convergence [26].
In model (1), k 1ˆCppt´1q describes the residual concentration on last day's pollution; illustrates the traffic PM 2.5 concentration from traffic stations by dispersion. The sum of these three components is allowed to decay with increasing wind by the factor e´k 5ˆWptq .
According to model (1), if W ptq " 0, which means there were no wind/dispersion at all, the model (1) reduces to its simplest format:Ĉ pptq " k 1ˆCppt´1q (2) To estimate the parameters in model (1), we obtained following data: daily PM 2.5 concentration of the six background stations C pptq , (p = 1, . . . , 6); daily PM 2.5 concentration of the two industrial stations C jptq .(j = 1, 2); daily PM 2.5 concentration of the five traffic stations C kptq .(k = 1, . . . , 5); distance from each background station to each industrial station D ind jp ; distance from each background station to each traffic station D tra f f ic kp ; daily maximum wind speed W ptq ; direction of daily maximum wind á m w , which is given in 16 compass points clockwise from the Y coordinate. Let á m w denote the unit vector of wind direction, and θ be the degree of direction from source station to receptor station p deviating from the Y coordinate, then: Let á R p denote the direction vector from the centroid of source stations to the receptor station p, then the summation of valid wind fluxŴ ptq to station p is given as: In our projection, we limited the minimum pollution brought by wind to nonnegative value, thus: Based on model (1), the daily traffic contribution to PM 2.5 at background stations can be calculated as: where T pptq % is estimated percentage of daily traffic contribution to total PM 2.5 concentration at background stations. Meanwhile, the expected daily non-traffic contribution NT pptq˚c an be calculated as: A GAMM was fitted between log transformed daily non-traffic PM 2.5 concentration logNT pptq and Y coordinates (Y p ) for the background stations. Because there were apparent nonlinear relationship between daily PM 2.5 concentration and day (t) (Figure 2a), humidity, temperature and atmospheric pressure (atmos) (Figure 2c), we used B-spline penalized by discrete penalties as additive smoothing function in the GAMM. Ten knots per year were set for day, five for humid, five for temperature and four for atmos, respectively. Numbers of knots were determined by minimizing Akaike Information Criterion (AIC) [27]. Besides the Y p , our preliminary analyses suggested linear associations between daily PM 2.5 concentration and wind speed (Wind (t) ), rain volume (Rain (t) ) and hours of daylight (Light (t) ) ( Figure 2d), they were included as covariates in the GAMM in addition to Y p . Day of week (DOW (t) ) and direction of daily maximum wind speed (Max_wind_dir (t) ) were included as factor variables in the model. In addition, considering the intra-cluster correlation of PM 2.5 concentration within stations, we included a random effect for stations in the model (Figure 2b). The selection of explanatory variables was also decided by a top-down rule [28]. The model was run by stepwise approach and generalized cross-validation (GCV) criterion [29]. The final GAMM is: Rain ptq`β5ˆM ax _ wind _ dir ptq`β6ˆD OW t s pt, k " 10 per yearq`s ptemperature ptq , k " 5q s phumid ptq , k " 5q`s patmos ptq , k " 4q`µˆZ p where logpNT pptq q˚is expected log transformed non-traffic PM 2.5 concentration; βs are parameters to be estimated; sp.qs are additive smoothing functions which illustrate the effects of day, temperature and humidity on non-traffic concentrations; Z p is a random intercept for station p. PM2.5 concentration within stations, we included a random effect for stations in the model ( Figure  2b). The selection of explanatory variables was also decided by a top-down rule [28]. The model was run by stepwise approach and generalized cross-validation (GCV) criterion [29]. The final GAMM is:  Log transformed non-traffic PM2.5 concentrations at non-background station q, were then predicted using the GAMM fitted in model (8). The estimated contribution of road traffic to PM2.5 contribution at non-background station q, ( ) % , was calculated as observed PM2.5 concentration deducted by estimated non-traffic PM2.5 concentration: Log transformed non-traffic PM 2.5 concentrations at non-background station q, logpNT qptq q˚, were then predicted using the GAMM fitted in model (8). The estimated contribution of road traffic to PM 2.5 contribution at non-background station q, T qptq %, was calculated as observed PM 2.5 concentration deducted by estimated non-traffic PM 2.5 concentration: The whole process of the method is shown in Figure 3. The parameters for dispersion model were estimated in software 1stOpt [26]. GLMM was fitted in Stata 13.1 and GAMM was fitted in R 3.2.2 using mgcv package.
Int. J. Environ. Res. Public Health 2016, 13, 124 8 The whole process of the method is shown in Figure 3. The parameters for dispersion model were estimated in software 1stOpt [26]. GLMM was fitted in Stata 13.1 and GAMM was fitted in R 3.2.2 using mgcv package.  There was a significant linear relationship between Y coordinates and log transformed PM2.5 concentrations both in all stations and in background stations (Figure 4), supporting our assumption that PM2.5 concentration followed an exponential decline function on distance. The Y coordinates could explain more than 80% variation of log transformed annual average PM2.5 concentrations in all stations. The closer a station was to the south border of the southern industrial area, the heavier the pollution level it had.
The optimal estimation of the parameters and fitness of the model was shown in Table 3. The dispersion model can explain more than 60% variation of the daily PM2.5 concentration of the background stations. The unexplained variation might on the other hand be due to temporal trend and meteorological conditions and was modeled in the GAMM later.  There was a significant linear relationship between Y coordinates and log transformed PM 2.5 concentrations both in all stations and in background stations (Figure 4), supporting our assumption that PM 2.5 concentration followed an exponential decline function on distance. The Y coordinates could explain more than 80% variation of log transformed annual average PM 2.5 concentrations in all stations. The closer a station was to the south border of the southern industrial area, the heavier the pollution level it had.

Results
The optimal estimation of the parameters and fitness of the model was shown in Table 3. The dispersion model can explain more than 60% variation of the daily PM 2.5 concentration of the background stations. The unexplained variation might on the other hand be due to temporal trend and meteorological conditions and was modeled in the GAMM later.        Based on Equation (6), the road traffic contribution to PM 2.5 concentration of the background stations is shown in Table 4. The contributions ranged from 17.2% in Yungang to 25.3% in Zhiwuyuan. The estimations of parameters and the approximate test of smoothing of GAMM are shown in Tables 5 and 6. All coefficients of the linear components and the smooth terms are significant at α = 0.05 level. The result is also in line with the fact that increasing pollution dilution was expected to be associated with greater wind speed and rain volume. According to Yu et al. [30], average PM 2.5 concentration during the days with wind speed higher than 2 m/s was 13% lower than those during the days with weaker wind. Average PM 2.5 concentration during the rainy days was 21% lower than those during the days without rain. But it is interesting that hours of daylight were negatively associated with the PM 2.5 concentration. This may be partly due to low dispersion rate during days with fewer daylight hours (usually in hazy and cloudy days) and accelerated accumulation of pollutants. The partial regression smooth plots (Figure 5b-e) and normal Q-Q plot of Pearson residual (Figure 5f) showed a good fit of GAMM. Based on Equation (9), the traffic contribution to PM 2.5 concentration of other stations is shown in Table 7. The absolute and relative contributions of road traffic to PM 2.5 concentrations of all stations were summarized in Figure 6. The average annual contribution of road traffic to PM 2.5 concentration ranged from 17.2% to 37.3% with a mean contribution 30%. The highest contribution was found in busy road areas, and the contribution in traffic-related stations is about 14% higher than those in rural areas.
Because there were no PM 2.5 values lower than one-tenth of the 25% percentile and only 5% values were higher than 10 times the 75% percentile, the estimated contributions changed little when including these abnormal values in sensitivity analysis (results not shown).

Discussion
Exhaust emissions due to road traffic are known to make a large contribution to total PM 2.5 concentrations in urban areas [32][33][34][35] and exposure to PM 2.5 from vehicular emissions has been demonstrated to have a negative impact on human health [36][37][38][39][40]. An improved understanding of the traffic-related contribution to PM 2.5 is therefore vital for conducting source apportionment and health effect studies. Due to rapid economic and industrial development and urbanization in the past few decades, energy consumption and the number of motor vehicles are rapidly escalating in China [41]. As the capital of China, Beijing has witnessed a devastating increase in air pollution in the past decades. To develop effective PM 2.5 reduction strategies, major sources of PM 2.5 and contributions from each source need to be understood thoroughly. A recent study claimed that vehicles had limited contribution to atmospheric particulate pollution in Beijing [42], and had since caused the public to question the governmental policy in limiting car use. The study presented PM 2.5 concentrations in all seasons in Beijing and concluded that vehicle emissions accounted for less than 4% of the total PM 2.5 [42], much smaller than the previous estimates of the Chinese Environmental Protection Agency or as reported by other studies [43][44][45][46]. Other studies using the same data sources suggested however that vehicle contribution to PM 2.5 in Beijing could vary between 10% and 50% [47,48].
Quantifying traffic-related contribution to PM 2.5 requires the compilation of detailed traffic data according to time and space, including, for example, traffic counts, vehicle types, travel speeds, fuel types, and emission controls [9]. Receptor models and air-quality dispersion models have been used to assess the contribution of different types of sources, including motor vehicles, to ambient pollution in urban and rural areas [49]. Traditionally, source apportionment estimation methods [50] such as chemical mass balance (CMB) [51] or positive matrix factorization (PMF) have been applied to analyze the contribution of pollutant source. Air mass trajectory analysis is also a useful tool for detecting the direction and location of sources for various air pollutants as a PM 2.5 forecast model [52]. However, these models heavily rely on the accuracy of source profile information. Some other models were also commonly used, mainly including source apportionment model [53], land use regression model and Gaussian dispersion model [54][55][56]. However, the limited numbers of roadside monitors have made it difficult to catch the geographical variation in motor-vehicle emissions. Resource requirements for collecting these data can be prohibitive and have led to the use of source-oriented dispersion based models [57], meteorological-chemical transport based models [58] and observation-based statistical models [59].
In our study, we developed a two-stage method to estimate the traffic-related contribution to PM 2.5 concentration that utilized the air-quality data from different types of AQM stations. This method combined atmospheric chemistry dispersion model and statistical GAMM model, and simplified the mathematical algorithm by omitting the detailed traffic-related information, e.g., types, number and density of vehicles, and incorporated the temporal trend of PM 2.5 concentration in a more precise way. We collected hourly PM 2.5 data at 35 monitoring stations to estimate the road traffic contributions to PM 2.5 concentrations. The results revealed that 17.2%-37.3% of PM 2.5 might be attributable to traffic emissions. Compared to the results released by Beijing Municipal Environmental Protection Bureau (22%-30%) [60], our reported contribution is higher and may partly be due to the rapid increase of traffic volume and decrease of industrial and coal burning emissions in recent years in Beijing [61].
Usually, the estimation of traffic-related emission relies on the analysis of road side measurements correcting for background concentrations [62]. In our study, we carefully defined the components of PM 2.5 concentration of background stations from two major sources, i.e., traffic emission and industrial sources. Considering the complex components of the traffic related PM 2.5 source at the traffic stations and industrial stations, relative to the background stations, we modeled the non-traffic PM 2.5 concentration for all stations using GAMM. The results from previous studies using particulate matter source apportioning and Comprehensive Air Quality Model with Extensions (CAMx) revealed that the maximum level of uncertainty for secondary production was low (6%), hence the application of an additive linear relationship was considered reliable [63,64]. In our dispersion model, the coefficients k 3 , k 4 , k 5 determine the precision of the estimated traffic contribution to PM 2.5 . We made simulation using different k 3 , k 4 , k 5 settings for the purpose of sensitivity analysis. The results showed that a 20% deviation in k 3 , k 4 , or k 5 would result in <7% change in the estimated traffic contribution. It indicated that our dispersion model was robust regarding the variation of the estimates of different parameters.
In order to avoid over-fitting or under-fitting, frequent in GAMM, we used penalized B-splines (P-splines). The P-spline approach controls the coefficients of the smooth function for which a certain penalty term is specified. In this approach, the crucial point is the selection of smoothing parameter. We tested the residual of the model and the scatter plots showed a clear homogeneity around smoothing curves with no specific trend (Figure 5b-e). In our model, the geographical variations were efficiently explained by Y coordinators. A few meteorological variables were selected in the models as previously suggested [65][66][67].
Our study has several strengths. First, most of the previous researches were performed in the United States or Europe, while reliable information from Africa, Asia and South America is lacking. Our study provides important evidence to fill in this information gap and offers an opportunity to develop enhanced methods for quantification of the contribution of traffic emission to air pollution. Second, the two-stage method predicted the background pollution instead of traffic emissions directly. In this case, the residual of the first dispersion model could be further decomposed in the GAMM and the unknown non-linear relationships and temporal autocorrelation were modeled using smoothing functions. Third, although existing dispersion models can give an approximate estimation of traffic emissions based on a big database, they need rich information in terms of vehicle types and fuels, traffic stop-and-go-driving situations, average speed and traffic density, etc. [68]. Moreover, the advanced Gaussian dispersion model also requires more complicated 3-dimensional meteorological and location information, making it unfeasible to adapt in less developed countries and regions. Our simplified dispersion model, on the other hand, needs less traffic and geographic data and applies simpler estimation algorithm, and therefore increases flexibility and feasibility of usage. In such context, it is a convenient tool on operational basis for estimating traffic contribution to PM 2.5 over a region with moderate number of AQM stations. Lastly, because of the limited number of AQM stations available, previous estimates of traffic contribution to PM 2.5 were mainly based on GAM that might not precisely reflect the variation between stations and correlation within stations in areas with various land use types [69]. The results of such studies were consequently very sensitive to the location of monitoring stations. However, the use of widespread AQM stations and intensive air quality data collected in our study made it possible to involve the different type of stations as a random factor in the mixed effect model that may sufficiently reflect the variation of contribution over a wide region.
Our study also had some limitations. Given the complexity of pollution sources and dynamic dispersion mechanisms, our simplified dispersion model only took into account industrial and traffic emissions, whereas it combined all other pollution sources as a whole. As a result, our method might have led to an overestimation of the traffic contribution. Although we examined the influence of daily average vehicle speed on PM 2.5 concentrations at five traffic stations and found no statistically significant association, this variable was not included in the GAMM since such information was not available for other stations. Finally, we did not consider some indirect sources from vehicles, such as tire type and asphalt roads that may also increase PM 2.5 concentration [70]. Future efforts are needed to compare methods using direct traffic emission measurements with our simplified indirect method. We also admit that the predictability of our models is not high and the accuracy of the estimated contributions needs to be assessed by further studies.

Conclusions
We developed a two-stage method to estimate the traffic contribution to daily PM 2.5 concentrations in Beijing using hourly PM 2.5 concentration data, daily meteorological data and geographic information collected at 35 AQM stations in Beijing between 2013 and 2014. Our results showed that traffic emissions accounted for a substantial share of total PM 2.5 concentrations, ranging from 17% at rural stations to 37% at stations close to busy roads. Our estimates were not only comparable to reports from the Beijing Municipal Environmental Protection Bureau but also reflected the spatial and temporal trends of traffic contribution in a large area. Lacking complete direct measurements of traffic emissions throughout the study area, this method fully utilized the characteristics of different station types. Our method is a useful and feasible tool in ecological and epidemiological studies to estimate the burden of PM 2.5 derived from road traffic when there is no sufficient traffic-related information.