Particulate matter of less than 2.5 μm(PM_2.5) is a critical issues in modern society, and the World Health Organization WHO (2018) has estimated that around 7 million people die every year from exposure to fine particles in polluted air. Numerous studies have been conducted to show the associations between exposure to ambient PM_2.5 and adverse health effects. Pui et al. (2014) reviewed various aspects of PM_2.5, including its measurement, source apportionment, visibility, and health effects, and mitigation, and Burnett et al. (2014) developed a fine particulate mass-based relative risk function. In addition, Song et al. (2017) estimated the health burden attributable to PM_2.5 based on three years of observed data.

South Korea is one of the worst countries in terms of severe air pollution, and many studies have focused on the PM_2.5 in South Korea. Choi et al. (2012) examined the characteristics, sources, and distributions of PM_2.5 and carbonaceous species in Incheon, South Korea, and Ryou et al. (2018) summarized the findings of PM source apportionment studies on South Korea. More recently, Bae et al. (2020) estimated long-term foreign contributions to the PM_2.5 concentrations in South Korea with a set of air quality simulations.

Various statistical models have been applied to the data to predict PM_2.5 level. Ordieres et al. (2005) compared three different topologies of neural networks for predicting average PM_2.5 concentrations: multilayer perceptron (MLP), radial basis function (RBF) and square MLP. Dong et al. (2009) proposed a model based on hidden semi-Markov models for high PM_2.5 concentration value prediction, and Qiao et al. (2019) developed a model based on wavelet transform–stacked autoencoder–long short term memory (LSTM).

These models perform well in the relatively usual event, but they work poorly to predict extremely high levels of PM_2.5. Since extreme events are rarely occurred, few data are used to learn the extreme patterns of the PM_2.5 level. Therefore, conventional statistical models does not perform well in predicting extreme events, if there is no assumption for such events in the models.

Various models have been developed to predict extreme values to overcome this issue. D’Amico et al. (2015) advanced the generalized Pareto distribution to model the probability distribution function’s tail to predict wind speed in Alaska. Quintela-del-Rı and Francisco-Fernández (2011) proposed a nonparametric functional data analysis to estimate ozone data in the UK, and Schaumburg (2012) combined nonparametric quantile regression with extreme value theory. However, few studies have been conducted to predict extreme PM2.5 levels (Qin et al., 2015).

In this paper, we adapt the three-stage model proposed by Wang and Li (2013) to predict extreme concentrations of PM_2.5 in Seoul, South Korea. The three-stage model extends the model by Wang et al. (2012a) that relaxed the assumptions of linear quantile functions of Y and tail equivalency across covariates, x. Wang and Li (2013) integrated quantile regression and extreme value theory by estimating intermediate conditional quantiles using quantile regression and extrapolating these estimates to tails based on the extreme value theory.

In this paper, we modify the three-stage model with lasso regression (Tibshirani, 1996) to improve the prediction performance of the extreme concentration of PM_2.5. Lasso regression performs variable selection, and provides a sparse solution. Therefore, even when numerous predictors are included in the initial model stage, penalized regression methods, such as lasso, provide a parsimonious linear model. There are various penalized quantile regression models (Wu and Liu, 2009; Alhamzawi et al., 2012; Wang et al., 2012b), and we consider the lasso quantile regression in this study. We further combine the three-stage models with the lasso penalty obtained from the various extreme quantile values, and generate a single prediction value. Finally, a binary classifier based on the proposed model is constructed.

The rest of the paper is organized as follows. The data description and exploratory data analysis are provided in Section 2. In Section 3, we propose a binary classifier based on a three-stage model with lasso regression, and an algorithm for implementation is presented in Section 4. In Section 5, we apply the proposed method to the PM_2.5 in Seoul, South Korea, and validate the results. Finally, the concluding remarks are presented in Section 6.

Table 1 lists the variables used in this study. The meteorological variables, hourly precipitation, temperature, wind speed, and humidity, were collected from the South Korea Meteorological Administration. Several studies indicate that the concentration level of PM_2.5 is affected by the meteorological factors. For example, Zhang et al. (2015) confirmed the critical role of meteorological parameters in air pollution formation, and Zhang et al. (2018) studied the influences of critical meteorological parameters, such as wind and precipitation, on PM concentrations. Therefore, we also included these variables in this study.

Air pollutant data were obtained from AIR KOREA, affiliated with the Korea Environment Corporation. The measurements include the hourly PM_2.5, PM₁₀, SO₂, NO₂, CO, and O₃ content. Stracquadanio et al. (2007) studied the correlations between PM_2.5 and gases (benzene, O₃, SO₂, NO₂, and CO). Song et al. (2015) used the generalized additive model to determine the statistical relationships between PM_2.5 concentrations and other pollutants, including SO₂, NO₂, CO, and O₃. Based on these studies, we also consider these pollutant variables to be covariates.

The variables in Table 1 are measured in 25 districts in Seoul, South Korea (Figure 1). Differences existed in places where the meteorological variables and air pollutant variables were measured, but the distance between the stations is relatively short.

All variables were obtained from January 1, 2015, to May 29, 2020, but in this study, we consider only the spring season (March, April, and May) every year, when the extremely high levels of PM_2.5 are observed. In addition, for the fast computation, we used four-hour average values for all variables. Therefore, the number of time points is 3298 for each district. For example, the four hour average values of PM_2.5 in the Gangnam district are presented in Figure 2.

Missing values occur in each measurements, thus, we imputed the missing points as follows. First, linear imputation was applied if length of sequential missing points was less or equal to four. Second, if the length of sequential missing points was greater than four and less than or equal to six, the sequence was imputed by the value of the nearest region within 5 km. Finally, the rows were omitted from the data if the previous steps did not impute the points.

To observe the variable correlations, we computed the correlation coefficients. Figure 3 displays the correlation in the Gangnam district, where PM₁₀, SO₂, and CO have a positive association of more than 0.5 with PM_2.5, a relatively strong association. However, others have weak relationships with less than 0.5. We also observed relatively strong associations between O₃, NO₂, CO, and wind speed, implying multicollinearity between variables. For the other districts, we observed a similar variable correlation.

The Ministry of Environment in South Korea categorizes the concentrations of PM_2.5 levels between 0 and 15 μg/m³ as ‘good’, between 16 and 35 as ‘normal’, between 36 and 75 as ‘bad’, and more than 76 as ‘very bad’ (Table 2). In this study, we focus on the forecasting ‘very bad’ cases in Seoul that causes severe health effects. The number of time points that categorized based on the PM_2.5 levels are presented in Table 3 for 25 districts. ‘Very bad’ cases are rarely observed for all districts, implying that conventional mean based statistical models may not work well in predicting extremely high PM_2.5 levels.

The proposed method is based on the three-stage model by Wang and Li (2013) with some modifications. We consider the response variable Y and covariate vector X = (X₁, . . . , X_p)^T with the X₁ = 1, and assume that we observe a random sample {(y_i, x_i), i = 1, . . . , n} of the random vector (Y, X).

The goal is to estimate the extremely high conditional quantiles for τ_n → 1 as n→∞,

where F_Y (·|x) is the conditional cumulative distribution function of Y.

Similar to Wang and Li (2013), we assume that F_Y (·|x) is in the maximum domain of attraction of an extreme value distribution G_γ(x₎, where γ(x) > 0 is the extreme value index. For the random sample Z₁, . . . , Z_n from F_Y (·|x), constants a_n(x) > 0 and b_n(x) ∈ ℝ exist such that,

where Z_(n) is the largest order statistic of the samples.

Although the three-stage model uses conventional quantile regression as a base model, we propose regularized quantile regression using the lasso penalty to overcome the multicollinearity problem and improve the prediction performance.

The entire procedure consists of four steps. First, the power transformation parameter λ of response variable Y is estimated. Secondly, the conditional intermediate quantiles are fitted to the transformed response variable. Then, the extreme quantile is estimated by extrapolating the intermediate quantile estimates and we ensemble the estimates using a simple average. Finally, the result is divided according to the threshold into two groups: ‘not very bad’ and ‘very bad’. The following section details the description of the procedure.

3.1. Power transformation

In quantile regression, conditional quantiles of Y are assumed to be linear in x at the tails. To relax this linearity assumption, Wang and Li (2013) considered the power transformation of Y. The power-transformed quantile regression model is defined as follows,

where τ ∈ [1 – ∊, 1], where ∊ is a small positive constant, and θ(τ) is the τ-th quantile regression coefficient.

The power transformation parameter λ in (3.4) is estimated as follows,

where Rn(t,λ;τ)=1/n Σi=1nI(xi≤t)[τ-I{Λλ(yi)-xiTθ^LASSO(τ;λ)≤0}], and I is the indicator function. The estimated lasso coefficient, θ̂^LASSO(τ; λ), is computed as follows,

where ρ_τ(x) = x{τ – I(x < 0)} is the τ^th quantile loss function, and ν is a penalty parameter in lasso regression estimated at each λ through using the cross-validation method (Tibshirani, 1996; Wu and Liu, 2009). Wang and Li (2013) suggested using the upper quantile level τ = 1 – ε with small positive constant ε, and we set τ = 0.95.

3.2. Estimating intermediate quantiles

In this step, we estimate the intermediate quantiles Q_Y (τ_j|x) for τ_j = j/(n + 1), j = 1, . . . ,m, with m = n – [n^η]. Parameter η is set as 0.1, as suggested by Wang and Li (2013), and [x] denotes the integer part of x.

The estimate of the intermediate quantile is

where λ̂ is a power transformation parameter estimated from (3.5), and

3.3. Extrapolating intermediate quantiles to tails

Next, we extrapolated the intermediate quantile estimates to the extreme tails. For τ_n → 1, we estimated Q_Y (τ_n|x) as follows,

where k = k_n →∞and k/n → 0.

The extreme value index γ̂_k(x) is estimated using the followings,

The selection of k is a crucial part of the three-stage model, and we applied a selection procedure, a modified version of Section 3.3 of Wang and Li (2013), by adding a lasso penalty term.

Then, we ensemble the extreme quantile results, Q̂_Y (τ_n|x), for , using a simple average,

where contains extreme quantile levels, and we set = {0.950, 0.955, . . . , 0.990, 0.995, 0.999}. | | denotes the number of elements in the set, which is equal to 11.

3.4. Binary classification procedure

In the last step, the estimated extreme quantile, Q̂_ensemble is classified into two classes based on the threshold value. This study focuses on the prediction of ‘very bad’ cases of PM_2.5; thus, we set the threshold as 76μg/m³, which is defined in Table 2.

We observed that {(y_i, x_i), i = 1, . . . , n}, where y_i is PM_2.5, and x_i = (x_1,i, . . . , x_9,i)^T is a covariate vector on the ith time. We used nine variables as covariates: PM₁₀, SO₂, NO₂, CO, O₃, precipitation, temperature, wind speed, and humidity.

Note that, we have used predicted explanatory variables. Therefore, proposed model can be written as,

where f (·) is proposed three-stage model and g(·) is the forecast model based on the exponential smoothing (ETS) algorithm (Hyndman et al., 2008) computed as follows,

where α is the smoothing parameter estimated by maximizing the likelihood.

The number of time points in the data set is 3, 298. We have few ‘very bad’ cases in the overall data; thus, a small test set contains few ‘very bad’ cases. Therefore, to validate the performance in predicting extreme cases, we have considered splitting the data into training and test sets at a ratio of 1:2. In this study, n₁ = 1,099 observations from March 1, 2015, to March 31, 2017, were used for the training data, and the test data are n₂ = 2,199 observations obtained from April 1, 2017, to May 29, 2020.

We constructed the proposed model using training data, and validated the model using test data. The training data are denoted as {(yitr,xitr),i = 1, . . . , n₁}, and the test data are {(yite,xite),i = 1, . . . , n₂}.

For each district in Seoul, we ran the following algorithm.

5.1. Evaluation metrics

The results were summarized using the sensitivity, specificity, positive predictive value, negative predictive value, and the F-score, which are statistical performance measures for a binary classification test. Based on the confusion matrix in Table 4, the following measures are defined,

• Sensitivity measures the proportion of correctly identified positives (the proportion of ‘very bad’ time points that are correctly identified as ‘very bad’),

  Sensitivity=TPTP+FN.

• Specificity measures the proportion of correctly identified negatives (the proportion of ‘not very bad’ time points that are correctly identified as ‘not very bad’)

  Specificity=TNTN+FP.

• The positive predictive value (PPV) is calculated as follows

  PPV=TPTP+FP.

• The negative predictive value (NPV) is calculated as follows

  NPV=TNTN+FN.

• The F-score measures overall accuracy is calculated as follows,

  F-score=(1+β2) TP(1+β2) TP+β2FN+FP=(1+β2)×PPV×Sensitivity(β2×PPV)+Sensitivity.

• The F-score is distributed from 0 to 1 and controlled by β, and it is chosen such that sensitivity is considered β times as important as the PPV. Considering our purpose is to accurately predict a ‘very bad’ event and the rareness of the event, we set β = 2 (Sasaki, 2007).

5.2. Comparison methods

Three conventional classification models are applied to evaluate the relative performance of the proposed model. As a representative non-parametric ensemble algorithm, the random forest (RF) is known to be robust to outliers and performs well in many classification problems. The MLP and LSTM models are artificial neural networks which are widely used in various applications. Significantly, variants of the LSTM have been used in predicting of PM_2.5 because it is appropriate for time series data.

The details of each model are summurized in Table 5. The hyperparameters are selected by trial and error. For example, we considered 400, 600, and 800 number of trees, and choose the 800 as the optimal value.

The commonly used threshold value which is used to classify the observation as positive is 0.5. However, the value may be inappropriate, considering the extremeness of the ‘very bad’ events in the data (Zou et al., 2016). Therefore, we considered threshold values from 0.01 to 1 in the test data and choose the result that provides the best F-score (Lakshmi and Prasad, 2014).

5.3. Results

The results from three districts, Guro, Yangcheon, and Yeongdeungpo, are presented in Table 6. The proposed three-stage model with lasso ensemble (TSLE) method provides the highest PPVs and comparable sensitivity values, and works best according to the F-score. Although the sensitivity values for the MLP are higher than those for the other models on for Yangcheon and Yeongdeungpo, the MLP has the lowest PPVs, which indicates that a false alarm frequently occurs. This outcome implies that the proposed method has relatively balanced performance for predicting ‘very bad’ events. However, all models perform well in terms of the NPV and specificity, because of the extremely imbalanced data.

In Figure 4, we plot the F-scores from the three-stage model with lasso regression (TSL) before ensemble. The specific quantile result may provide the best result but there is inconsistent. For example, TSL prediction when τ_n = 0.990 has the highest F-score in Guro, whereas the highest value occurs when τ_n = 0.999 in Yangcheon. Therefore, the ensemble technique solves this selection problem, and offers superior performance overall.

We obtained similar results for all 25 districts (not shown), and we plotted the average measures of all districts in Figure 5. Overall, the proposed TSLE method works best for all five measures.

As the lasso regression performs the variable selection, some variables are not selected in the regression model. Therefore, we present the proportion of the selection for each variable in Table 7. Temperature, humidity, and PM₁₀ were selected in most models, and NO₂, O₃ and SO₂ were rarely selected.

In this paper, we consider the prediction of extreme values of PM_2.5 in Seoul, South Korea. Compared to the conventional mean-based models, the proposed method is based on the quantile regression with the extreme value theory. Therefore, the proposed model predicts the extremely high values of PM_2.5 especially well. Moreover, we added the lasso penalty term to the quantile model, so it performs variable selections. Based on the statistical measures of the performance of a binary classification test, such as the sensitivity, PPV, and F-score, the proposed method works well for 25 districts in Seoul, South Korea.

We expected that the proposed model performance would improve by adding meteorological variables in China, which significantly affect the atmospheric conditions in South Korea. In addition, the data file and R code for implementation are provided at https://github.com/SaeSimcheon/Extreme-PM2.5-prediction.