Insighting Drivers of Population Exposure to Ambient Ozone ( O 3 ) Concentrations across China Using a Spatiotemporal Causal Inference Method

: Ground-level ozone ( O 3 ) is a well-known atmospheric pollutant aside from particulate matter. China as a global populous country is facing serious surface O 3 pollution. To detect the complex spatiotemporal transformation of the population exposure to ambient O 3 pollution in China from 2005 to 2019, the Bayesian multi-stage spatiotemporal evolution hierarchy model was employed. To insight the drivers of the population exposure to ambient O 3 pollution in China, a Bayesian spatiotemporal LASSO regression model (BST-LASSO-RM) and a spatiotemporal propensity score matching (STPSM) were ﬁrstly applied; then, a spatiotemporal causal inference method integrating the BST-LASSO-RM and STPSM was presented. The results show that the spatial pattern of the annual population-weighted ground-level O 3 ( PWGLO 3 ) concentrations, representing population exposure to ambient O 3 , in China has transformed since 2014. Most regions (72.2%) experienced a decreasing trend in PWGLO 3 pollution in the early stage, but in the late stage, most areas (79.3%) underwent an increasing trend. Some drivers on PWGLO 3 concentrations have partial spatial spillover effects. The PWGLO 3 concentrations in a region can be driven by this region’s surrounding areas’ economic factors, wind speed, and PWGLO 3 concentrations. The major drivers with six local factors in 2005–2014 changed to ﬁve local factors and one spatial adjacent factor in 2015–2019. The driving of the trafﬁc and green factors have no spatial spillover effects. Three trafﬁc factors showed a negative driving effect in the early stage, but only one, bus ridership per capita (BRPC), retains the negative driving effect in the late stage. The factor with the maximum driving contribution is BRPC in the early stage, but PM 2.5 pollution in the late stage, and the corresponding driving contribution is 17.57%. Green area per capita and urban green coverage rates have positive driving effects. The driving effects of the climate factors intensiﬁed from the early to the later stage.


Introduction
Ground-level ozone (GLO 3 ) is a well-known major atmospheric pollutant threatening human health, vegetation and biodiversity worldwide [1][2][3].GLO 3 is also an important greenhouse gas contributing to climate change [4].Significant increasing trends of annual GLO 3 concentrations occurred in most global urban areas from the 1990s, particularly from 2005 to 2014, and annual GLO 3 concentrations increased at a fast rate in East Asia, North America, and Europe [5].China has established several air pollution control policies and programs [6,7] since 2016.Consequently, the concentrations of most air pollutants, including PM 2.5 , NO 2 and SO 2 , of China have decreased since 2013, but GLO 3 concentrations have shown an increasing trend [8].GLO 3 may have become the second most hazardous pollutant after PM 2.5 [9].Lyu et al. [10] pointed out that the ratio of the population living in regions with high PM 2.5 pollution decreased from 2015 to 2021; however, this ratio remained almost unchanged for the ozone and even increased.
Considering the complex GLO 3 concentration problems in China, a deep understanding of the spatiotemporal trends of O 3 concentrations in China can provide helpful evidence for environmental epidemiology and appropriate regional control strategies of GLO 3 .Some researchers have investigated the spatiotemporal characteristics and influencing factors of GLO 3 in a specific city or a limited spatial region in a short-term period, mainly using the in situ monitored O 3 concentrations data.For instance, An et al. [11] explored the temporal variations of O 3 concentrations in Nanjing, Eastern China, and Wang et al. [12] studied the spatial and temporal patterns of observed ground-level O 3 concentrations in China.The above research found O 3 concentrations to be highest in summer and lowest in winter, with significant diurnal variations.Zhang et al. [13] used the spatial-temporal kriging model to analyze the daily maximum 8 h ozone concentration data from 55 air quality monitoring stations in China's Pearl River Delta Metropolitan Region in 2015, revealing the O 3 concentrations have a spatial autocorrelation that undergoes temporal migration and conversion across the seasons.Tang et al. [14] and Gao et al. [15] underscored the influence of anthropogenic emissions, particularly in urban areas, on ozone levels.Liu et al. [16] suggested that weather patterns and regional interactions should be considered when designing ozone mitigation strategies.Lyu et al. [17] demonstrated the use of machine learning algorithms for predicting ozone levels, showing promising results.These findings emphasize the need for region-specific and season-specific strategies to effectively address the issue of ozone pollution in China.Some scholars have explored the space-time variations of the O 3 concentrations across China using simple descriptive statistical methods based on in situ monitored air quality data.Yang et al. [18], Xue et al. [19], and Li et al. [20] explored the temporal variation and spatial distribution of GLO 3 concentrations in China during [2013][2014][2015][2016][2017][2005][2006][2007][2008][2009][2010][2011][2012][2013][2014][2015][2016][2017], and 2015-2020, respectively.Based on multisource data and machine learning, Liu et al. [21] and Wei et al. [22] estimated and investigated briefly the temporal and spatial trends of the GLO 3 concentrations in China during 2005-2017 and 2013-2020.The spatial distribution of the population was not considered in the above research.However, the health exposure risk of GLO 3 concentrations is determined by GLO 3 concentrations and population spatial distribution jointly.Given this, our study focused on the population exposure to the GLO 3 concentrations, instead of semplice GLO 3 concentrations.
The formation regime or driver of GLO 3 concentrations is also a vital issue.A large number of studies have investigated the complicated relationship between GLO 3 production and main precursors, volatile organic compounds (VOC) and NO x , from the perspective of atmospheric physicochemistry [23][24][25][26][27].Some scholars have researched the spatiotemporal characteristics of GLO 3 formation regime of local areas of China, such as Beijing [28][29][30] and Shanghai [31,32].Li et al. [20] and Lu et al. [33] studied the spatiotemporal variations of GLO 3 formation regimes across China.Although it is well known that meteorological factors have an influence on O 3 concentrations, the specific influencing mechanism still needs to be further explored.The meteorological factors' influencing effects on GLO 3 concentrations were studied in some articles [34][35][36][37][38].These studies all found that temperature was the dominant driver of GLO 3 concentrations at the annual time scale.
Most of the existing research has focused mainly on the pure GLO 3 concentrations without considering the population spatial distribution, and the study phases of these studies were within a short time range.Additionally, to our knowledge, few studies have investigated the simultaneous influence of natural and anthropogenic factors on population exposure to GLO 3 concentrations across China.Given this, our paper utilised the advanced Bayesian spatiotemporal hierarchy model to investigate transitions of spatial distributions of the annual PWGLO 3 concentrations over China in recent 15 years, and then presented a spatiotemporal causal inference method for the multiple-factor analysis scenario, which was used to identify the potential drivers of population exposure to GLO 3 concentrations.The outcome variable in this study is population exposure to GLO 3 concentrations.In this study, population exposure to GLO 3 concentrations is represented by the annual population-weighted GLO 3 (PWGLO 3 ) concentrations calculated by the grid data of the GLO 3 concentrations and population distribution by the following formula,

Methods
where GLO 3,ikt and Pop ikt represent the GLO 3 concentrations and distributed population count of the k th grid with spatial resolution of 10 km × 10 km in the i th prefecture region.I ik is the total number of the 10 km grid of across the i th prefecture region.The grid data of the GLO 3 concentrations at a spatial resolution of 10 km × 10 km was collected from the CHAP dataset (https://weijing-rs.github.io/product.html(accessed on 6 June 2022)), which was estimated using the extended space-time extra-trees (STET) model [21] from big data including ground-based observations, remote sensing products, atmospheric reanalysis, and emission inventory [21].The GLO 3 estimates are highly consistent with the ground-based measurements with R 2 = 0.96 and the root-mean-square error (RMSE) of 8.64 µg/m 3 .Additionally, the GLO 3 estimates for the whole of China are highly consistent with surface measurements (R 2 = 0.87), and the values of the mean RMSE, the mean absolute error (MAE), and the mean relative error (MRE) over the entire domain, were 17.10 µg/m 3 , 11.29 µg/m 3 , and 18.38%, respectively [21].The grid data of the population density with the spatial resolution of 1 km × 1 km were collected from an open database, WorldPop data (https://www.worldpop.org/(accessed on 6 June 2022)).In consideration of the inconformity between WorldPop data and those reported in the statistical yearbook, the WorldPop data employed a ratio coefficient method to perform the correction based on the population data of the statistical yearbook.The population weight of each grid can be calculated by the total population of the grid divided by the total population of the prefecture-level region covering the spatial grid.Then, the PWGLO 3 concentrations of each prefecture-level region may be calculated with the GLO 3 concentrations multiplying by the population weight and zonal accumulation.

Influencing Factors
The influencing factors of the PWGLO 3 concentrations mainly came from three aspects [39] illustrated by Figure 1: social economics, traffic development, and natural environment.Each category driver contains a certain number of proxy variables.The social economics involved 7 proxy variables: population density (PD), gross domestic product per capita (GDPPC), tertiary industry proportion (TIP), real estate investment per capita (REIPC), industrial electricity consumption per capita (IECPC), electricity consumption of residents per capita (ECRPC), and total social electricity consumption per capita (TSECPC).It should be noted that, in China, the tertiary industry mainly refers to the service industry, while the secondary industry refers to the mining, manufacturing, electricity (heating power, fuel gas, and water) production and supply, and construction industries.The traffic development included 3 proxy variables: bus ridership per capita (BRPC), number of taxis per capita (NTPC), and number of buses per capita (NBPC).The data for the above 10 proxy variables were collected from the China city statistical yearbook in corresponding years.The natural environment contained 9 proxy variables: urban green coverage rate (UGCR), urban construction area per capita (UCAPC), ratio of urban construction area (RUCA), green area per capita (GAPC), difference of summer average temperature from that in 2004 (DSAT), annual PM 2.5 concentrations (APMC), annual cumulative duration of sunshine (ACDS), annual average wind speed (AAWS), and annual average relative humidity (AARH).It should be noted that the DSAT was selected in order to underline the rising temperature.As previously mentioned, VOC and NO x , the emissions of human activities that mainly include industry and traffic, are the common precursors of ozone.The three traffic factors in this paper include public traffic, bus, and taxi.Additionally, the factor pertaining to industry emissions lacked data availability.Consequently, the APMC, as a proxy variable of the emissions of the industry and private cars, was absorbed.Moreover, UGCR, UCAPC, RUCA, and GAPC, were collected from the China city statistical yearbook in corresponding years (https://data.cnki.net/trade/yearbook/Single/N2022040095?zcode=Z024 (accessed on 6 June 2022)).The raw grid data of annual average PM 2.5 and temperature with a spatial resolution of 10 km × 10 km were downloaded from the National Earth System Science Data Center (http://www.geodata.cn/data/(accessed on 6 June 2022)).The data of the three climate factors, ACDS, AAWS, and AARH, were collected from the China Meteorological Data Service Centre (http://data.cma.cn/en(accessed on 6 June 2022)).

Bayesian Multi-Stage Spatiotemporal Evolution Hierarchy Model (BMSSTEHM)
The BMSSTEHM presented by Li et al. [40] was employed to explore the spatiotemporal trends of the  concentrations from 2005 to 2019 in mainland China.The BMSSTEHM integrated the BSTHM [41] and piecewise regression model [42] and can investigate non-linear local trends by capturing the self-adaptive turning points of local

Bayesian Multi-Stage Spatiotemporal Evolution Hierarchy Model (BMSSTEHM)
The BMSSTEHM presented by Li et al. [40] was employed to explore the spatiotemporal trends of the PWGLO 3 concentrations from 2005 to 2019 in mainland China.The BMSSTEHM integrated the BSTHM [41] and piecewise regression model [42] and can investigate non-linear local trends by capturing the self-adaptive turning points of local trends and considering spatial correlations and heterogeneity.The mathematical forms of the BMSSTEHM are as follows: (2) where c it represents the observed annual average PWGLO 3 concentrations of the i-th prefecture-level region at year t.µ it and σ 2 c represent the mean and variance of the likelihood distribution.α is an intercept term.S i is the parameter of spatial fixed effect; (b 0 t + v t ) describes overall trends containing a linear b 0 and non-linear tendency v t , whose prior distribution adapted Gaussian distribution; b 1i , b 2i are local piecewise linear regression coefficients; a 1i is the first turning point of the i-th prefecture-level region.If there is one turning point, the spatiotemporal can be divided into the first (early) and second (later) stages.The corresponding linear variation parameters are k 1i and k 2i .G t,a 1i is the logistic function, as described in Equation (3).λ is a shape parameter that is generally assigned to a value greater than 10.T is the number of time nodes (T = 15 in this study).Linear regression is necessary only when there are more than two sample values, namely at least three sample values.a 1i is therefore limited to between 3 and T − 4.Then, it can guarantee that there are at least three sample values on both linear sides.a 1i is assigned to the prior uniform distribution and can be self-adaptively estimated according to the sample data.The term ε it is a Gaussian noise error whose prior distribution is assigned as a normal distribution N 0, σ 2 ε .The spatial relative magnitude of the PWGLO 3 concentrations of the i-th prefecture-level region in the two stages, denoted as SR 1i and SR 2i , quantifying the PWGLO 3 polluted level relative to the global overall level can be estimated.SR 1i and SR 2i can be expressed as follows: where t * 1 and t * 2 represent the temporal centre points of the early and later stages, respectively.In this study, M is the number of prefecture-level regions of mainland China, which is 338.The parameters, S i , b 1i , b 2i , and a 1i , are all considered simultaneously spatial structured and unstructured effects by assigning the prior of the BYM model [43].Therefore, the 'Queen' spatial adjacency structure sharing a common edge or a common point was adopted in this paper.ε 1 and ε 2 represent Gaussian error terms.
The Bayesian statistics calculations in this paper were implemented by WinBUGS 14.0 [44] based on the Markov chain Monte Carlo (MCMC) algorithm.The number of iterations was set to 200,000, 180,000 were the burn-in period, and 20,000 were the sampling number of the posterior distribution of parameters.The convergence of the Bayesian inferences in this study was monitored by standard autocorrelation plots and trace plots.

Bayesian Spatiotemporal LASSO Regression Model (BST-LASSO-RM)
For spatiotemporal process, the spatial endogenous and exogenous interaction effects should be considered [45][46][47].Considering possible multicollinearity problems due to the no little number of explanatory variables, this study presented a BST-LASSO-RM based on the Bayesian LASSO regression model [48] to overcome the problem of multicollinearity, considering the spatial endogenous and exogenous interaction effects.The regression parameters of the BST-LASSO-RM were assigned to Laplace priors [49].The BST-LASSO-RM can acquire a more stable estimation and automatically provide interval estimates for all parameters, including the error variance [48].The mathematical structure of the BST-LASSO-RM can be expressed as follows: where µ it and σ 2 y represent the mean and variance of the normal likelihood distribution of the observed PWGLO 3 concentrations.φ is an intercept term.W i presents the spatial adjacent matrix of the i-th prefecture-level city and adopts Queen contiguity that includes common vertices and edges.W i µ it and W i x j,i,t represent the spatial adjacent average PWGLO 3 concentrations and the j-th explaining variable around the i-th prefecture-level region in the t-th year.ρ is the parameter of the spatial endogenous effect of the annual PWGLO 3 concentrations, W i µ it .β j is the regression parameters of local influencing factors, x j,i,t., and ϑ j is the parameter of spatial exogenous effects of spatial adjacent influencing factors, W i x j,i,t ; k is the number of the explanatory variables, ξ it represents the random error, β, ρ, and θ represent the estimate of the regression parameters β, ρ and ϑ. µ and X are the matrices of µ it and x j,i,t , respectively.λ β , λ ρ , and λ ϑ are the shrinkage parameters.The Bayesian calculations of the BST-LASSO-RM inferences were implemented mainly by the Python package, PyMC3.

A Spatiotemporal Propensity Score Matching (STPSM) Method
Similarly, to take into account the spatial endogenous and exogenous interaction effects, a STPSM method was presented based on the PSM method [50] in this study.This method can be adopted to recognise the causality between the PWGLO 3 concentrations, Y, and the spatial adjacent PWGLO 3 concentrations, WY, the local and spatial adjacent influencing factors, X and WX.Matched sets of treated and untreated groups that share a similar value of the propensity score can be formed by the STPSM, and then the average treatment effect for the treated, , is the probability of treatment assignment, Z i = 1, conditional on other observed covariates, X i .The propensity score of one factor can be obtained by logistic regression as follows: Considering that all the variations in this study are continuous numeric variables, the sample will be divided into treated and untreated groups according to whether the target variate is greater than its median.STPSM is then employed to implement 1:1 matching, in which pairs of treated and untreated groups are formed with similar values of the propensity score.

A Spatiotemporal Causal Inference Method
This study proposed a method for recognizing spatiotemporal causality and causeeffect, integrating STPSM and BST-LASSO-RM.As shown in Figure 2, this method consists of three steps.First, the BST-LASSO-RM was employed to recognise the significant correlated, Λ (1) , or not significant correlated variables, Λ (0) , with the outcome variable, Y, from all independent variables, X, WY, and WX.
where Y i and y i are the observed value and Z standardized score, Y mean and Y std are the corresponding mean and standard deviation.Specifically, the representations of the W i µ it , x j,i,t and W i x j,i,t are the same as the above.And the three type of variables can be identified as the significant correlated variable if the posterior probability of the corresponding regression coefficients, P(ρ > 0|Data) or P(ρ < 0|Data), P(β j > 0|Data) or P(β j |Data) < 0, P(ϑ j > 0|Data) or P(ϑ j |Data) < 0, are greater than 1 − α, e.g., 0.95, and vice versa.Second, the STPSM method was used to recognise the causality of the significant correlated independent variables, Λ (1) including WY, X, and WX, to outcome variable, Y. Let us suppose there are K(K ≤ k) significant correlated variables, Λ j (j = 1, 2, . . ., K), and then follow the steps below to recognise the causality of each of the Λ  (i) Let Λ (1) 1,i,t ; (ii) Then, the observed outcome, Y it , can be divided into two groups, they can be called treatment and control group, according to the value of Λ (1) k,i,t ; (iv) Using STPSM to match the sample data between the treatment and control groups; (v) Adopt t-test to identify the significant difference between the treatment and control groups and consequently to recognise the causality of Λ M .
The identified causal and not causal variables were denoted as Λ (11) and Λ (10) l (l = 1, 2, . . ., K − M), it identifies that all significant correlated influencing factors are causal factors to explained variable, Y, if K − M = 0. Through the above steps, possible three types of causal effects, spatial endogenous causal effect from WY, spatial exogenous effect from WX, and local variational causal effect from X, can be identified.Third, again employing the BST-LASSO-RM to estimate the quantitative cause-effects, δ Λ (11)

The Spatial Pattern and Its Transformation of the PWGLO 3
The spatiotemporal process can be generally divided into two phases called early and late stages.The overall spatial distributions of the PWGLO 3 concentrations in the two phases can be identified by the BMSSTEHM (Figure 3).The results show that the overall spatial pattern of the PWGLO 3 concentrations transformed from the early to late stage.Specifically, the spatial distribution of the PWGLO 3 concentrations in the early stage is disperser than that in the late stage.The high PWGLO 3 polluted regions transferred from the eastern and southern coastal areas in the early stage into the northern and North China Plain regions, including Hebei, Shanxi, Shandong, Henan, Jiangsu, and Anhui, in the late stage.The composition of the number of prefecture-level regions with the five classifications of the PWGLO 3 levels changed between the two stages.The numbers of the high PWGLO 3 polluted prefecture-level regions with spatial relative magnitude (SRM) of 1.05~1.20 and greater than 1.20 changed from 110 (32.5%) and 13 (3.8%) in the early stage to 94 (27.8%) and 68 (20.1%) in the late stage.The number of the medium PWGLO 3 polluted prefecture-level regions decreased from 96 (28.4%) in the early stage to 64 (18.9%) in the late stage.The number of the lightly PWGLO 3 polluted prefecture-level regions with SRM of less than 0.95 remained stable during the study period, 2005-2019.
Furthermore, the SRM of PWGLO 3 concentrations over the southeast coast areas of China decreased from 1.05~1.20 in the early stage to less than 1.05 in the late stage, while that over the North China Plain regions increased to greater than 1.20.The PWGLO 3 concentrations SRM over the northern regions, such as Shanxi, Inner Mongolia, Gansu and Ningxia, increased from less than 1.05 in the early stage to 1.05~1.20 and greater than 1.20 in the late stage.

The Local Trends of the 𝑃𝑊𝐺𝐿𝑂 Concentrations at the Sub-Provincial Level
The annual change in the  concentrations in the two stages at the prefecture level can also be estimated by the parameters,  and   , of the BMSSTEHM.Research shows that the  production is a photochemical oxidation process determined highly by meteorological conditions and precursors, such as VOC and  .Consequently, variations in the annual mean  concentrations could be modulated by many factors, especially climate factors, such as air temperature and sun radiation, which can elevate or lower annual mean  concentrations.To validate the reliability of the annual temporal trends estimated by the BMSSTEHM, the scatters and BMSSTEHM fitted       In the late stage, the spatial pattern of the local trends of the  concentrations formed a distinct spatial structure (Figure 6).The  concentrations in the most sub-provincial areas (79.3%) showed an increased trend.In particular, some cities (dark red coloured areas in Figure 6) of central China, such as Ningxia, Shanxi, Henan, Anhui, and south Gansu, experienced a strong increased trend of greater than 4.00 μg/m per year.The median, lower, and upper quartiles of all the prefecture-level regions with increasing local trends are 2.82, 1.19, and 4.38 μg/m per year.The local trends of the areas located in Tibet, Xinjiang, Guangxi, and Guangdong converted from upward in the early stage to downward trend in the late stage.In the late stage, the spatial pattern of the local trends of the PWGLO 3 concentrations formed a distinct spatial structure (Figure 6).The PWGLO 3 concentrations in the most sub-provincial areas (79.3%) showed an increased trend.In particular, some cities (dark red coloured areas in Figure 6) of central China, such as Ningxia, Shanxi, Henan, Anhui, and south Gansu, experienced a strong increased trend of greater than 4.00 µg/m 3 per year.The median, lower, and upper quartiles of all the prefecture-level regions with increasing local trends are 2.82, 1.19, and 4.38 µg/m 3 per year.The local trends of the areas located in Tibet, Xinjiang, Guangxi, and Guangdong converted from upward in the early stage to downward trend in the late stage.Secondly, recognising the causality of the significant influence factors on the  concentrations.As mentioned in Section 2, the STPSM method was adopted to recognise causality.Considering that all the variables are numeric types for each significant influence variable, the prefecture-level regions can be bounded by the median to be  Secondly, recognising the causality of the significant influence factors on the PWGLO 3 concentrations.As mentioned in Section 2, the STPSM method was adopted to recognise causality.Considering that all the variables are numeric types for each significant influence variable, the prefecture-level regions can be by the median to be divided into the treatment group and the control group.Then, the prefecture-level regions in the treatment group and control group were matched by the STPSM.Subsequently, the relative risk (RR) and risk difference (RD) of each selected significant influence variable on the explanatory variable, the PWGLO 3 concentrations, can be (Table 1).Furthermore, the RR and RD can be calculated by the ratio and difference of the average PWGLO 3 concentrations in the matched treatment group and control group.The results show that most of the selected significant influence variables possessed statistically significant (p < 0.05) causality with the PWGLO  Figure 7 shows the normalised regression coefficients of the BST-LASSO-RM between the PWGLO 3 concentrations and the driving factors in the two stages.The results indicate that the structure of the drivers of the PWGLO 3 concentrations has changed from the early stage to the late stage.Concretely, the eight factors, GDPPC, TSECPC, BRPC, NTPC, DSAT, ACDS, ARRH, and three spatial adjacent factors, W•PWGLO 3 , W•GDPPC, W•TSECPC, had driving effects in both stages.However, the driving direction of the NTPC factor changed from negative to positive.This is possibly due to the fact that the use efficiency of taxis was reduced from the early to the late stage due to the increase in private cars.Moreover, GAPC and NBPC had driving effects in 2005-2014, excluding 2015-2019; however, TIP, W•TIP, ECRPC, W•ECRPC, UGCR, W•APMC, AAWS, and W•AAWS were on the contrary.

Drivers of the 𝑃𝑊𝐺𝐿𝑂 Concentrations
finding is that GAPC has a positive driving effect in only one stage.The reason for this could be that GAPC decreased from the early stage to the late stage.The negative driving effects in 2005-2014 for BRPC, NTPC, and NBPC, and 2015-2019 for BRPC, indicate that the development of public transportation, especially bus ridership, can alleviate PWGLO 3 Six local factors, BRPC, APMC, TSECPC, GDPPC, NBPC, and GAPC, were the major drivers in 2005-2014, and the corresponding explanatory powers, the percentages of the absolute value of the normalised regression coefficients, were 13.56%, 11.95%, 9.82%, 8.72%, 8.46%, and 8.27%.Five local factors and one spatial adjacent factor, APMC, ACDS, TSECPC, BRPC, AARH, and W•APMC, the major drivers with explanatory powers of 17.57%, 9.75%, 9.71%, 7.05%, 6.28%, and 8.94%, in 2015-2019.

Implications of the Results
The results of the transitions of the overall spatial pattern and local trends of the annual PWGLO 3 concentrations can clarify the spatiotemporal evolution mechanism of the PWGLO 3 concentrations in China.This indicates that the emphatically concerned regions should focus on the areas with higher levels and increasing trends of PWGLO 3 concentrations (i.e., the red and dark red areas in Figures 3B and 5).The identification of the potential drivers regarding the prefecture level of the PWGLO 3 concentrations in China can help policy-makers design more effective policies to control PWGLO 3 concentrations.The annual PWGLO 3 concentrations can be reduced by decreasing the factors with positive driving effects and increasing the factors with negative driving effects.Specifically, although local and spatial adjacent social economic factors (GDPPC, W•GDPPC, TIP, W•TIP, ECRPC, W•ECRPC, TSECPC, and W•TSECPC) have positive effects, these factors cannot be easily decreased due to the development of the economy.Nevertheless, the annual PWGLO 3 concentrations can be lowered by decreasing NTPC and increasing BRPC simultaneously.The increase in UGCR can result in the growth of the annual PWGLO 3 concentrations, and green coverage may provide many human health benefits.This finding suggests that the vegetation, with more VOC emissions as precursors of ozone, should be cut down in the urban area.The results estimated by the spatiotemporal causal inference method show that the spatial adjacent PWGLO 3 , i.e., W•PWGLO 3 , has positive driving effects in the two stages, 2005-2014 and 2015-2019, this states that the PWGLO 3 concentrations have a significant spatial causal spillover effect.Namely, the local PWGLO 3 concentrations in a prefecture-level region may also be aggravated by its spatial adjacent PWGLO 3 concentrations.Meanwhile, considering that the AMPC and W•APMPC have positive driving effects, this indicates that surface ozone pollution may be synergistically controlled with local and circumjacent PM 2.5 pollution, i.e., synergistic governance.The climate drivers, such as temperature, sunshine, wind speed, and relative humidity, can also not be changed; however, the related results may help us understand the variation of PWGLO 3 concentrations in the background of global climate change.

Discussion
This paper focused on the problem of population exposure to surface ozone concentrations, measured by the annual PWGLO 3 concentrations, at a prefecture-level scale across China mainland during 2005-2019.In this study, an advanced method called BMSSTEHM was used to carefully detect the space-time trends the annual PWGLO 3 concentrations.Furthermore, to explore the drivers of the PWGLO 3 concentrations, this study transcended the constraints of traditional correlation analyses by presenting a spatiotemporal causal inference method that discerns causality relationship and estimates the cause-effect.This method combines STPSM and BST-LASSO-RM, both of which are also proposed in this paper.The methods in this paper provide a methodological reference for other similar studies.And the spatiotemporal causal inference method can be applied potentially in analysing remote sensing data.Furthermore, the results of the paper can help in the creation of diverse regional policies to reduce population exposure to GLO 3 concentrations.
The spatial distribution of PWGLO 3 concentrations in mainland China has shifted since 2015.It is important to carry out in-depth research on why the spatial pattern of PWGLO 3 concentrations on China's mainland has changed.The spatial structure of the local trends of the annual PWGLO 3 concentrations at prefecture-level in China's mainland showed differences in the two stages.The results indicated that the annual PWGLO 3 concentrations in the early stage maintained a steady state, i.e., most prefecture-level regions experienced a low local trend.But in the later period, the regions with high PWGLO 3 concentrations level experienced simultaneously high local trends.This means that the PWGLO 3 concentrations of the hot spot regions will become more severe.Additionally, the areas with high levels of PWGLO 3 concentrations have the tendency to expand outwards.
The transformation reasons for the static spatial pattern and local trends of the annual PWGLO 3 concentrations are very intricate.Generally, the reasons can be classified as anthropogenic and meteorological influences.Some researchers [51][52][53] have concluded that meteorological factors were not negligible but anthropogenic factors were dominant.This study argues that the economic development and remaining higher PM 2.5 concentrations jointly lead to the transformation of spatiotemporal trends of Chinese PWGLO 3 concentrations.According to Wei et al.'s newest study [54], the PM 2.5 concentrations across China decreased significantly from 2015-2020, but those in North China still remained at a higher level; in contrast, those in South China reached a lower level.In particular, the Pearl River Delta region of China has carried out pioneering pollution control measures since 2014 (https://www.mee.gov.cn/ywgz/xcjy/shxc/201909/t20190909_733028.shtml(accessed on 6 January 2023)), been "withdrawn" from the three key regions of air pollution prevention and control in China.The results of this study indicated that the PM 2.5 as a precursor of ozone has a maximal driving effect on PWGLO 3 in the late stage.Moreover, the rapid economic development of North China in the late stage may be also another reason.
The formation mechanism of the PWGLO 3 concentrations is very complex.Our study has preliminarily explored the drivers of the annual PWGLO 3 concentrations at the prefecture-level scale.The driving pattern of the annual PWGLO 3 concentrations in China's mainland showed different features in the two phases.The four local and spatial adjacent economic factors, GDPPC, W•GDPPC, TIP, W•TIP, ECRPC, W•ECRPC, TSECPC, and W•TSECPC, all exhibited positive driving effects on the PWGLO 3 concentrations in the late stage.This indicates that the PM 2.5 concentrations were reduced by the shrinkage of the secondary industry due to environmental regulations; however, the development of the tertiary industry led to increasing GLO 3 concentrations.Moreover, not only the local economic factor can drive the PWGLO 3 concentrations, but also of the corresponding surrounding areas also can.This implies that the process of adjusting the policies aimed to alleviating surface ozone pollution must consider regional coordination.
In terms of traffic development factors, increasing bus ridership can promote the reduction in the PWGLO 3 concentrations in the two stages, and the driving effects of the BRPC in the early and later stages came close.It is interesting that the increase in NTPC and NBPC can lower the 3 concentrations in the early stage; however, in the later stage, NBPC had no significant driving effect and the NTPC's driving direction changed from negative to positive.The reason for the above phenomenon might be that travel demand within the city has increased rapidly due to the of the tertiary industry.But the transport efficiency of the taxi is analogous to that of the private car, and less than of the bus.The low transport efficiency of the taxi may lead to high per capita exhaust gas emission, a precursor of ozone pollution.This suggests that public traffic should be energetically developed to improve public transport volumes and then reduce the GLO 3 concentrations.
This study found that GAPC and UGCR had positive driving effects in the early and later stages, respectively.This result is consistent with the conclusion of related environmental chemistry research [55].Furthermore, the GLO 3 is mainly generated by photochemical reactions of NO x and VOCs that were released from anthropogenic activities [56], while the VOCs are an emission from vegetation [55].Hence, green areas or green covers become risk factors for the GLO 3 concentrations.Our study quantified the relationship between the annual PWGLO 3 concentrations and GAPC and UGCR.Eucalyptus, Populus, Cunninghamia, etc., are plants with high VOC emissions; however, thickets and grasslands have low VOC emissions [57].This result suggests that plants with high VOC emissions should be decreased in the urban green coverage.
For the factors of traffic development and green cover, the corresponding driving effects only occurred in the local area.That is, the factors of traffic development and green cover in the spatial adjacent areas have not driving effects on the local area's PWGLO 3 concentrations.This finding provided important evidence or reference to control surface ozone pollution by changing the local factors of the traffic development and green cover.
From the view of ecological studies, this paper also demonstrated the conclusion that temperature and PM 2.5 are important drivers on the GLO 3 concentrations.We evaluated the quantitative driving effects of the local DSAT and APMC, and spatial adjacent PM 2.5 , W•APMC, on annual PWGLO 3 concentrations at the prefectural scale, which is different from an environmental chemistry study.This can provide valuable evidence for policymaking for mitigating PWGLO 3 concentrations in China and other global regions.
This paper contributes by providing not only the concrete statistical results of the spatiotemporal trends of the PWGLO 3 concentrations at a sub-provincial scale, but also a frame for recognising causality and estimating causal effects for the situation with multifactors, and investigating the drivers of the PWGLO 3 concentrations in mainland China.However, this study also has some limitations.First, the study period was not comprehensive enough, and the statistical results would be more robust if the study period was long enough.Second, the spatial scale was adopted with sub-province, but this can still be refined, such as the county scale.The above limitations will be emphatically focused on.

Conclusions
First, this study proposed a spatiotemporal causal inference method combining STPSM and BST-LASSO-RM, both of which are also proposed in this paper, and this method can be applied potentially in analyzing remote sensing spatiotemporal data.Second, the spatial distribution of PWGLO 3 concentrations in mainland China has shifted since 2015.Third, while early-stage annual PWGLO 3 concentrations remained stable across most prefecture-level regions in China's mainland, the later stage the regions with high PWGLO 3 concentrations also exhibited high local trends.Fourth, a shift in the drivers of PWGLO 3 concentrations occurred from the early to late stages, the factors like NTPC changing their driving direction.Notably, the development of public transportation, especially bus ridership, can mitigate PWGLO 3 concentrations, with six local factors being the primary drivers in 2005-2014 and five local factors, along with one spatial adjacent factor, dominating in 2015-2019.Fifth, the annual PWGLO 3 concentrations can be lowered by decreasing NTPC and increasing BRPC simultaneously.Synergistic governance with local and adjacent PM 2.5 pollution can help control surface ozone pollution in the context

22 Figure 1 .
Figure 1.The schematic diagram between outcome variable (dependent variable) and determinants and its proxy variables.

Figure 1 .
Figure 1.The schematic diagram between outcome variable (dependent variable) and determinants and its proxy variables.

Figure 2 .
Figure 2. Flow diagram of the spatiotemporal causal inference framework.Figure 2. Flow diagram of the spatiotemporal causal inference framework.
(vi)  Repeat the above five steps, (i)-(v), to recognise the causality of other significant correlated independent variables, Λ(1)2 , . . ., Λ Figure S1 shows the temporal series of the spatial distribution of PWGLO 3 annual concentrations in mainland China from 2005 to 2019.This reveals that the spatial pattern of the PWGLO 3 annual concentrations in mainland China has changed since 2014.The number of prefecture-level regions, with the annual PWGLO 3 concentrations greater than 80 µg/m 3 , increased from 172 in 2014 to 244 in 2019.The PWGLO 3 concentrations in the North China Plain region have become even more serious since 2014.In other words, different spatial structures of the PWGLO 3 concentrations were formed during 2005-2014 and 2015-2019.

Figure 3 .
Figure 3.The overall spatial pattern of the  concentrations in the early stage (A) and latestage (B) and the transformation (C) of the composition of the five classes of  levels between the two stages.

Figure 3 .
Figure 3.The overall spatial pattern of the PWGLO 3 concentrations in the early stage (A) and latestage (B) and the transformation (C) of the composition of the five classes of PWGLO 3 levels between the two stages.

3. 3 .
The Local Trends of the PWGLO 3 Concentrations at the Sub-Provincial Level The annual change in the PWGLO 3 concentrations in the two stages at the prefecture level can also be estimated by the parameters, b 1i and b 1i + b 2i , of the BMSSTEHM.Research shows that the GLO 3 production is a photochemical oxidation process determined highly by meteorological conditions and precursors, such as VOC and NO x .Consequently, variations in the annual mean PWGLO 3 concentrations could be modulated by many factors, especially climate factors, such as air temperature and sun radiation, which can elevate or lower annual mean GLO 3 concentrations.To validate the reliability of the annual temporal trends estimated by the BMSSTEHM, the scatters and BMSSTEHM fitted polylines of the annual PWGLO 3 concentrations in the example of 12 prefecture-level regions illustrated in Figure 4.The 12 prefecture-level regions were randomly selected while ensuring coverage in all directions of China.The results show that the fitted polylines can capture different linear trends in the two stages, although the annual PWGLO 3 concentrations in some years were saltatory or abnormal, possibly due to unstable climate factors.It should also be noted that the parameters of the BMSSTEHM, b 1i and b 1i + b 2i , measured the average annual change within the corresponding two stages; thus, the mutations in some years had little influence on the two parameters.Remote Sens. 2023, 15, x FOR PEER REVIEW 12 of 22 capture different linear trends in the two stages, although the annual  concentrations in some years were saltatory or abnormal, possibly due to unstable climate factors.It should also be noted that the parameters of the BMSSTEHM,  and   , measured the average annual change within the corresponding two stages; thus, the mutations in some years had little influence on the two parameters.

Figure 4 .
Figure 4.The scatters and Bayesian multi-stage spatiotemporal evolution hierarchy model fitted polylines of the annual population-weighted ozone ( ) concentrations in the example of 12 prefecture-level regions in China from 2005 to 2019.

Figure 5
Figure 5 illustrates the local trends of the  concentrations in the early stages.Generally, most prefecture-level regions (72.2%) experienced a decreased local trend of the  concentrations in the early stages.The corresponding median, lower, and upper quartiles of the annual change are −0.36,−0.53, and −0.19 μg/m per year.The other 27.8% of the prefecture-level regions are distributed dispersedly across China's mainland.The median, lower, and upper quartiles of the annual increasing changes are 0.20, 0.08,

Figure 4 .
Figure 4.The scatters and Bayesian multi-stage spatiotemporal evolution hierarchy model fitted polylines of the annual population-weighted ozone (O 3 ) concentrations in the example of 12 prefecture-level regions in China from 2005 to 2019.

Figure 5 22 and
Figure 5 illustrates the local trends of the PWGLO 3 concentrations in the early stages.Generally, most prefecture-level regions (72.2%) experienced a decreased local trend of the PWGLO 3 concentrations in the early stages.The corresponding median, lower, and upper quartiles of the annual change are −0.36,−0.53, and −0.19 µg/m 3 per year.The other 27.8% of the prefecture-level regions are distributed dispersedly across China's mainland.The median, lower, and upper quartiles of the annual increasing changes are 0.20, 0.08, and 0.64 µg/m 3 per year.Some sub-provincial regions located in Guangxi, Guangdong, and Tibet possessed high PWGLO 3 polluted levels and increased local trends concurrently in the early stage.moteSens. 2023, 15, x FOR PEER REVIEW 13 of 22

Figure 5 .
Figure 5.The local annual change of the  concentrations in China at the sub-provincial scale in the early stage.

3 Figure 5 .
Figure 5.The local annual change of the PWGLO 3 concentrations in China at the sub-provincial scale in the early stage.

Figure 6 .
Figure 6.The local annual change in the  concentrations of China at a sub-provincial scale in the late stage.

3 Figure 6 .
Figure 6.The local annual change in the PWGLO 3 concentrations of China at a sub-provincial scale in the late stage.

Figure 7 .
Figure 7. Normalised regression results of the Bayesian LASSO regression model between the  concentrations and the significant driving factors in the two stages.

Figure 7 .
Figure 7. Normalised regression results of the Bayesian LASSO regression model between the PWGLO 3 concentrations and the significant driving factors in the two stages.
2.1.Variable and Data 2.1.1.Population Exposure to GLO 3 . Results 3.1.Descriptive Statistics Overall, Chinese ground-level O 3 concentrations maintain a slight downward trend from 2005 to 2014, but an increasing tendency arose after 2014.The median of ground-level O 3 annual concentrations decreased from 82.73 in 2005 to 81.57µg/m 3 in 2014, but increased to 87.84 µg/m 3 in 2019.The PWGLO 3 annual concentrations maintained relatively steady from 2005 (81.80 µg/m 3 ) to 2014 (81.06 µg/m 3 ), an increasing trend occurred from 2014 to 2019 (89.92 µg/m 3 ).The heterogeneity of the PWGLO 3 concentrations at the prefectural level in China showed an overall rising trend.The corresponding coefficient of variation (standard deviation divided by mean) increased from 14.16% in 2005 to 15.60% in 2014, and to 17.48% in 2019.

Table 1 .
The estimated relative risk (RR) and risk difference (RD) of the significant influence factors in the two periods, 2005-2014 and 2015-2019, based on the STPSM eliminating confounding effects.

Table 2 .
Non-normalised regression results of Bayesian LASSO regression model between the PWGLO 3 concentrations and the significant driving factors in the two stages.