Ambient air pollutants concentration prediction during the COVID-19: A method based on transfer learning

2.1. Transfer learning

Given a labeled source domain $D_s={\left\{x_i,y_i\right\}}_{i=1}^n$ and an unlabeled target domain $D_t={\left\{x_j\right\}}_{j=n+1}^{n+m}$ or a labeled target domain with relatively little data $D_t={\left\{x_j,y_j\right\}}_{j=n+1}^{n+m}$, data distribution of these two domains, $P\left(x_s\right)$ and $P\left(x_t\right)$, are different, i.e., $P\left(x_s\right)\ne P\left(x_t\right)$. The TL is used to find the similarities between these two domains, thereby achieving the knowledge transferring from the source domain to the target domain to learn the labels [16]. According to the classification of learning methods, the existing TL method can be categorized into instance-based TL, feature-based TL, parameter-based TL, and relation-based TL [17], [18].

Instance-based TL focuses on how to select instances from the source domain that are useful for training in the target domain [19]. This type of TL is based on the precondition that there is a similar distribution between source and target domains. Then it can achieve the knowledge transformation by re-weighting samples in the source domain and applying the available information to the target domain. For example, Kim and Lee [20] proposed a new domain adversarial neural network, which can modify the original source domain data and convert it into an auxiliary target domain. By incorporating the attention mechanism into TL, He et al. [21] employed the source domain to create samples for the training of the target domain. Chen et al. [22] adopted the weight updating scheme to obtain valuable samples from the source domain, thereby reducing the effort on the source domain. In our study, the source domain refers to the observations of ambient air pollutant concentration before the COVID-19 epidemic and the target domain refers to the observations during COVID-19. Although there are sample distribution differences between these two domains, they are all time series data related to pollutant emissions. Therefore, the instance-based TL is suitable for the prediction task of this study and we will follow the instance-based transfer strategy to carry out APCP during the COVID-19 epidemic.

Feature-based TL focuses on how to find the common feature representation between the source domain and target domain, and then employs these features for knowledge transformation [23]. The feature-based TL emphasizes effective feature analysis like feature selection, mapping, and encoding. This type of TL method generally includes symmetric and asymmetric feature transformation [24]. The former aims to find valuable features across domains, and the latter reduces the domain difference by transforming the features of the source domain into the target domain. The commonly used transfer component analysis (TCA) is based on symmetric feature transformation, which can discover the representations of cross-domain features by minimizing the differences between marginal distributions [25]. Parameter-based TL assumes that the source domain and target domain share some model parameters or have the same prior distribution. This type of method aims to find these same model parameters or prior distributions to achieve knowledge transformation [26]. For example, the single-model knowledge transfer learns both the knowledge of the target domain and the transfer knowledge in the parameters of the pre-trained model to achieve the prediction task of the target domain [24]. Relation-based TL assumes that if two domains are similar, they will share a similar relationship. Specifically, this method uses the source domain to learn the logical relation network and then applies it to the target domain to achieve knowledge transformation [27]. However, this type of TL method is limited in practice because it involves the complex relational map between source and target domains.

By comparing different types of TL methods, this paper chooses the instance-based TL to obtain results of APCP during COVID-19. Considering that the existing TL methods for APCP have ignored the detailed domain analysis and comparison, this paper introduces a distribution description method, GMM, to obtain the distribution of each domain for further domain analysis and knowledge transformation. A detailed introduction to GMM is shown in the following subsection.

2.2. Gaussian mixture model

The TL method aims to reduce the distribution discrepancy between domains, and it mainly relies on the description of domains with appropriate distributions. For example, the distribution of atmospheric pollutant concentrations is uncertain and may change at any time, and it is therefore impossible to directly measure the difference between sample distributions. GMM is a popular method to explore the distribution structure of samples, which adopts Gaussian distribution to quantify the sample and decompose the sample into several Gaussian sub-models [28]. The domain distribution quantified by GMM can well support the discrepancy calculation between samples, and as a result, it provides a basis for distribution description and knowledge transfer procedure [29]. GMM has been widely used in TL research due to its excellent capabilities for sample distribution analysis [30]. In this study, we utilize this method to obtain the distribution description of historical observation. Specifically, GMM is a simple extension of the Gaussian distribution, which can be regarded as a mixed model composed of K Gaussian sub-distributions (namely, hidden variables). The distribution of GMM can be represented as:

where ${\alpha }_k$ is the probability that the observation data belongs to the kth sub-distribution, ${\alpha }_k\ge 0$ and ${\sum }_{k=1}^K{\alpha }_k=1$. $\varphi \left(x\left|{\theta }_k\right)\right)$ is the Gaussian distribution density function of the kth sub-distribution. And ${\theta }_k=\left({\mu }_k,{\Sigma }_k\right)$, where ${\mu }_k,{\Sigma }_k$ represent the mean and covariance matrix of the sample in the kth sub-distribution. $K,{\alpha }_k,{\mu }_k$ and ${\Sigma }_k$ are the parameters needed to be solved. The expectation maximization (EM) algorithm is the commonly used method for parameter training of GMM [31]. This is an iterative algorithm for estimating the maximum likelihood of the parameters of a probability model with hidden variables. Steps for updating parameters of GMM through EM iteration are as follows:

Step 1: Initialize the parameters.

Step 2: Calculate the probability of data j belongs to sub-distribution k, which is represented as $r_{jk}$ and calculated by:

Step 3: Calculate the model parameters for the new iteration.

Step 4: Repeat the above two steps until convergence. Thus, we can obtain the parameters of GMM.

After obtaining the domain distribution by GMM, we further focus on the distribution difference between domains. Commonly used distribution discrepancy measure methods include correlation alignment (CORAL) [32] and MMD [33]. The former realizes deep domain adaptation of knowledge by reducing the difference of covariance matrix between two domains. The latter is a popular unsupervised pattern recognition method, which can calculate the distribution difference between domains by matching appropriate feature representation and kernels [34]. Comparing these two methods, we find that MMD is more compatible with the domain distributions represented by GMM. In addition, the distance calculation of MMD can support parameter optimization to minimize the difference between domains. Hence, we employ MMD to determine the difference between the two distributions.

2.3. Shapley additive explanations and partial dependence plot

This subsection introduces two types of feature analysis methods to analyze the impact of influential factors on final results, which are helpful to reveal some valuable management implications. First, SHAP is used to rank the Shapley values of features, namely the average marginal contribution of features, to obtain the importance of variables [35], which can be represented by:

where $y_i$ represents the prediction result and $y_{base}$ represents the baseline for the prediction model (usually the mean of the target variables for all samples). $x_i$ represents the ith sample, and $x_{ij}$ represents the jth feature of the ith sample. $g\left(x_{ij}\right)$ represents the SHAP value of $x_{ij}$, that is, $g\left(x_{i1}\right)$ is the contribution value of the first feature to the final results $y_i$ in the ith sample. $g\left(x_{i1}\right)>0$ indicates that this feature shows a positive effect and $g\left(x_{i1}\right)<0$ represents the negative feature. That is, the results of SHAP generate both feature importance and the polarity of influence, positive or negative [36].

Second, PDP is another type of feature analysis method to further show how features affect the prediction results. That is, different from other feature importance analysis methods, which only focus on the feature importance or the polarity of influence, PDP can show the detailed influence of features on the prediction. This method can help to show how the prediction results vary with the changing of features by capturing the marginal effect of features [37]. Specifically, we can observe the relationship between prediction target and input features, such as whether it is linear, monotonic, or more complex. In addition to obtaining the relationship between the prediction target and input features, the PDP can also be adopted to analyze the relationships between different features.

3.1. Domain distribution comparison

Considering that the distribution of domains greatly influences knowledge transfer, we first explore the distribution comparison of the source domain and target domain. Despite the simple comparison of distribution in Fig. 1, this study employs a useful tool, the Q–Q plot, to show the specific distribution comparison between domains. This plot takes the percentile of each value in the sample data set as the abscissa value, and the percentile of the value in the reference data set as the vertical axis. It employs a 45-degree reference line to visualize distribution differences: if two sample sets come from a population with the same distribution, the sample points should fall near this reference line; on the contrary, the greater the departure from this reference line, the greater the evidence for the conclusion that the two data sets are distributed differently.

This study takes one of the important pollutant indicators, AQI, as an example to show the Q-Q plot in the source domain and target domain, as illustrated in Fig. 2. In this study, the source domain refers to the pollutant observation before the COVID-19 epidemic and the data span is from January 1, 2018, to January 23, 2020. The target domain refers to the pollutant observation during the COVID-19 epidemic and the selected data span is from January 24, 2020, to July 31, 2020. Detailed data introduction is in the following Section 4.1. In Fig. 2, we observe that both the distribution of points and the red line fitted by points deviate significantly from the blue reference line. We also obtain the $p$-value of the Kolmogorov–Smirnov test, which is a popular method to test whether a group of samples comes from a probability distribution or compare the distribution of the two groups of samples [38]. In this study, the obtained result of the Kolmogorov–Smirnov test is 3.66e−15. In the case of a confidence level of 5%, this result indicates that there is no sufficient reason to explain that the two datasets are subject to the same distribution, which also verifies the results shown in the Q-Q plot. That is, there are distribution differences between these two kinds of domains, and as a result, the traditional ML methods are not suitable for the prediction task. Therefore, this paper introduces an improved TL method combining the analysis of data distribution to achieve the prediction.

3.2. The improved TL prediction model based on GMM

4.1. Data collection

This paper collects monitor data for two years and seven months (ranging from 2018.01.01 to 2020.07.31) from the monitoring station in Wuhan, China. Websites for the monitor data collection include China National Environmental Monitoring Centre (CNEMC) (http://www.cnemc.cn) and National Climatic Data Center (NCDC). The collected data includes the concentration of main ambient air pollutants and some meteorological variables, which are commonly used variables in existing APCP studies. The collected ambient air pollutants/indicators include AQI, PM2.5, PM10, SO2, NO2, O3, and CO. These pollutants/indicators are the key factors used by monitoring stations to report air quality and are also the main reference information used in many existing studies on air pollutants. All the mentioned indicators in addition to AQI are measured by real-time concentration and 24-h moving mean (variables related to O3 include the real-time concentration, 8-h moving mean, 24-h moving mean, and the 24-h maximum of the 8-h sliding mean). The collected meteorological variables include some factors related to atmospheric pollutant concentrations such as temperature, relative humidity, windspeed, wind direction, air pressure, air pressure trend, and amount of precipitation. To simplify the data without losing important information, we collect data at three-h intervals starting from 2:00 every day. And the missing values can be filled by crawling data from adjacent monitoring sites at the same time. A detailed description of the collected data is shown in Table 1.

In this paper, AQI is taken as an example to carry out the prediction. That is, AQI is the dependent variable to be predicted, and other variables in Table 1 are independent variables. Note that AQI data is a kind of time series, the lagged AQI can also be applied as input variables in the prediction model. Considering that most AQI predictions are 3–10 days, we take the commonly used 7 days as the prediction cycle and select the AQI value at the same time 7 days earlier as the additional feature input. In this way, we obtain 7544 pieces of data from 2018.01.01 to 2020.07.31 and 22 features for the prediction task.

4.2. Prediction evaluation method and experiment setting

Given the collected dataset, this paper employs some commonly used methods to evaluate prediction performance, including mean square error (MSE), mean absolute error (MAE), explained variance score (EVS), and R2_score. The calculation of MSE is:

where N represents the number of samples. $A_i$ and $F_i$ represent the actual value and prediction value, respectively. The smaller the value of MSE, the better the fitting effect. MAE is another evaluation method by calculating the average of absolute error between the actual value and prediction value, and the function of MAE is:

The MAE is commonly used to assess the closeness of the predicted results to the actual value, and the smaller the value, the better the fitting effect. EVS is used to explain the variance scores of regression models, and the value of EVS is $\left[0,1\right]$. The closer the value of EVS is to 1, the more independent variables can explain variance changes of dependent variables. And the smaller the value, the worse the effect. R2_score refers to the judgment coefficient, which is also the variance score of the regression model. Its value range is $\left[0,1\right]$, and the larger the R2_score, the more independent variables can explain the variance change of dependent variables and vice versa.

We then introduce the experimental settings of this study. All experiments involved in this study are based on Python 3.7. The source codes of the main experiment in this study can refer to https://github.com/chenny1996/APCP-based-on-TL. To satisfy the predictive requirements, the proposed APCP model requires prospective validation, that is, the test set should be isolated from model tuning and forward in time. In this study, we divide the data set in Section 4.1 into a training set containing data from the first 25 months (i.e., source domain, 2018.01.01 to 2020.01.23) and a test set containing data from the last 6 months (i.e., target domain, 2020.01.24 to 2020.07.31). We first employ GMM to deal with the original domains and then obtain the new source domain by the proposed optimized model. Note that all the mentioned analysis of domain distribution in this study only contains the 22 input features without the predicted target AQI. Then this study trains the selected ML prediction models on the new source domain with these input features. Models estimated from the new source domain are then applied to the target domain for AQI prediction during the COVID-19 epidemic. All samples for training are normalized by the minimum–maximum method to ensure the comparable performance of different models.

Like other ML methods, the selected prediction models except LR have hyperparameters that need to be determined during model training, such as the penalty value in LASSO and the maximum depth in GBR. This study determines the optimal hyperparameters for each model by using grid search and 5-fold cross-validation. Then we train the selected models on the target domain and evaluate the prediction performance using the four aforementioned methods: MSE, MAE, EVS, and R2_score. To verify the validity of the proposed model and statistically compare different models, this study employs the n-out-of-n bootstrap sample based on the source domain to evaluate each model (including the proposed method and comparison models) 30 times. The n-out-of-n bootstrap is an emerging method for deriving test statistics under large samples [45], [46], and the comparison results with statistical significance can be obtained by combining this method with paired t-tests. The calculation of paired t-tests of n-out-of-n bootstrap is as follows:

Suppose that to compare the predictive performance of model A and model B, $N$ is the number of n-out-of-n bootstrap subsets, $P^A=\left[P_1^A,P_2^A,\dots ,P_N^A\right]$ and $P^B=\left[P_1^B,P_2^B,\dots ,P_N^B\right]$ are the obtained prediction performance of model A and model B on all subsets. $P_i=\left(P_i^A-P_i^B\right)$ represents the performance difference between model A and model B on the ith bootstrap subset. The paired t-test calculates the following t statistic for the null hypothesis that the mean difference $\overline{P}=\frac{1}{N}{\sum }_{i=1}^N{P}_i$ is equal to zero:

4.3. Results analysis

This study introduces the improved TL prediction model to obtain the results of APCP during COVID-19. We first obtain the sample distribution of the source domain and target domain by GMM. By employing the BIC method, we can obtain the number of clusters in GMM for the source domain and target domain, as shown in Fig. 4 and Fig. 5, respectively. The color-coded boxes in these two figures represent different covariances in GMM. The bar charts marked with * represent the best selection for each domain. We find that the optimal number of clustering of these two domains is 5, and the full covariance performs better.

The remaining parameters of GMM, i.e., parameters of sub-distributions and corresponding weights, can be obtained by EM optimization. Then we can obtain the sample distribution of the source domain and target domain represented by GMM. To visualize the obtained distribution of these samples, we select the first two variables, PM2.5 and PM2.5_24 h, as coordinates to visualize the clustering results of the predicted target AQI. The clustering result of AQI in the source domain is shown in Fig. 6. Due to the partial overlap of the sample distribution, the 3D graph in Fig. 6 still cannot show the results of each cluster. Five clusters of the source domain are separately shown in Fig. A.1, Fig. A.2, Fig. A.3, Fig. A.4, Fig. A.5. The clustering result of AQI in the target domain by GMM is shown in Fig. 7, and separate visualizations are shown in Fig. A.6, Fig. A.7, Fig. A.8, Fig. A.9, Fig. A.10, respectively. To clearly show the data distribution, the mentioned graphs are based on unnormalized data, and the following predictions are based on normalized data to compare different models.

Within the obtained sample distribution of the source domain and target domain, we then compare the similarity of these two domains and obtain the new source domain similar to the target domain by Eqs. (9), (10). And finally, training selected base models in Section 3.2.2 in the adjusted source domain and applying the well-trained model to the target domain. The results obtained by the proposed TL-based ML models are shown in Table 2. The numbers in bold represent the best, and the standard deviations of the 30 samples of the bootstrap test set are in parentheses, the same as below. In these prediction models, we observe that the prediction errors of all models are small and EVS and R2_score are close to 1, which verifies the effectiveness of the proposed model. Further analysis of the table reveals that LR and its variants, BR and LASSO, show little difference in results, while ENR generates the worst results. Among them, the GBR model performs the best, which can be attributed to the effect of model integration in GBR.

4.4. Comparison analysis

5.1. Results and discussion of SHAP analysis

Using the open-source package shap in Python, we can obtain the results of SHAP for this study. The top 20 important features obtained by SHAP are shown in Fig. 8. In this figure, the samples are represented by points with different colors. The redder the point color represents the greater the feature value, and the bluer the point color represents the smaller the feature value. The ordinate of Fig. 8 represents features and the abscissa is the SHAP value of features, i.e., the influence of features. In this figure, we find that some features show positive effects on prediction, and the representative features include PM10, PM2.5, and O3. And some of the features show no significant effects, that is, the increase or decrease of these types of features will not directly lead to the change of the influence on the results, or the corresponding change has no obvious rule. For example, as for the feature O3_24 h, we observe that when the feature value is small, the SHAP value is close to 0. But when the feature value is large, the changes of SHAP have no obvious rule. And as for some meteorological features, such as wind direction and windspeed, both the increase and decrease of feature value show no significant influence on SHAP.

To obtain a clearer understanding of how each feature affects the final results, Fig. 9 supplements the importance of features obtained by calculating the mean SHAP value. Combining these two figures, we can conclude: (1) Compared with meteorological features, pollutant-related features show more influence on the prediction of AQI; (2) Among pollutant-related features, PM10 and PM2.5 show more significant influence than other pollutants; (3) Some meteorological features that are generally considered important, like windspeed and amount of precipitation, show no significant effects on the prediction. While some meteorological features which are easily ignored, such as air pressure trends, have a great influence on the prediction results. The above results can help us understand the role of each feature in the prediction, and then the PDP is adopted to analyze how these features affect the final results.

5.2. Results and discussion of PDP analysis

This subsection selects four relatively important features in Fig. 9, PM10, PM2.5, O3, and air pressure trends, to analyze how they affect the prediction. Using the open-source package in Python, pdp, we obtain the PDPs of these four features, as shown in Fig. 10, Fig. 11, Fig. 12, Fig. 13. The $x$-axis of these figures represents the feature and the $y$-axis represents the change of predicted values. The blue shaded part represents the confidence interval. In Fig. 10, Fig. 11, Fig. 12, we observe a monotonically increasing linear relationship between the three pollutant-related features and the predicted target. This result indicates that the influence of these features on the predictive target is almost linear, and this influence increases with the increase of the feature value. These findings can provide guidance for managers to make decisions. For example, given the absolute contribution of PM10 and PM2.5 to the prediction, these two variables should be monitored more rigorously to accurately predict AQI values during COVID-19. In addition, the positive linear effects of these variables on prediction remind us to pay more attention to the major source of PM10 and PM2.5 shown in Table 1 and take corresponding control measures to alleviate air pollution. While in Fig. 13, we observe that the positive effect of air pressure trends on prediction is short-lived and disappears when the variable increases to a certain extent. This finding supplements Fig. 9 about how air pressure trends specifically affect the prediction.

The above analysis draws a conclusion about the relationship between each feature and the prediction target. Then we further focus on the analysis of different features to uncover the effect of controlling these features on other features. This study selects the feature with the greatest impact, i.e., PM10, and explores the relationship between PM10 and the other three main features. We employ the PDP to obtain the results, as shown in Fig. 14, Fig. 15, Fig. 16. In these figures, the $x$-axis and $y$-axis represent the features to be analyzed, and the $z$-axis represents the influence on the prediction results. In Fig. 14, we observe that when one of PM10 or PM2.5 is fixed, the value of another feature changes will also lead to the corresponding impact on the prediction results. In Fig. 15 and Fig. 16, we observe that when the values of PM10 are constant, the change of O3 and air pressure trends do not cause significant changes in the results. On the contrary, the change in PM10 can lead to significant changes in the results. These observations indicate that both PM2.5 and PM10 have a positive correlation with the predicted results, while O3 and air pressure trends show no significant positive effects on the prediction when fixing PM10. Therefore, we can conclude that controlling PM10 and PM2.5 has a greater impact on the prediction, and more attention should be paid to these two factors.