Evaluation of Machine Learning Models for Estimating PM_2.5 Concentrations across Malaysia

Abstract

Southeast Asia (SEA) is a hotspot region for atmospheric pollution and haze conditions, due to extensive forest, agricultural and peat fires. This study aims to estimate the PM_2.5 concentrations across Malaysia using machine-learning (ML) models like Random Forest (RF) and Support Vector Regression (SVR), based on satellite AOD (aerosol optical depth) observations, ground measured air pollutants (NO₂, SO₂, CO, O₃) and meteorological parameters (air temperature, relative humidity, wind speed and direction). The estimated PM_2.5 concentrations for a two-year period (2018–2019) are evaluated against measurements performed at 65 air-quality monitoring stations located at urban, industrial, suburban and rural sites. PM_2.5 concentrations varied widely between the stations, with higher values (mean of 24.2 ± 21.6 µg m⁻³) at urban/industrial stations and lower (mean of 21.3 ± 18.4 µg m⁻³) at suburban/rural sites. Furthermore, pronounced seasonal variability in PM_2.5 is recorded across Malaysia, with highest concentrations during the dry season (June–September). Seven models were developed for PM_2.5 predictions, i.e., separately for urban/industrial and suburban/rural sites, for the four dominant seasons (dry, wet and two inter-monsoon), and an overall model, which displayed accuracies in the order of R² = 0.46–0.76. The validation analysis reveals that the RF model (R² = 0.53–0.76) exhibits slightly better performance than SVR, except for the overall model. This is the first study conducted in Malaysia for PM_2.5 estimations at a national scale combining satellite aerosol retrievals with ground-based pollutants, meteorological factors and ML techniques. The satisfactory prediction of PM_2.5 concentrations across Malaysia allows a continuous monitoring of the pollution levels at remote areas with absence of measurement networks.

1. Introduction

Air pollution has become an acute environmental and health issue in developing countries during the last decades due to intense industrialization and urbanization processes [1,2,3]. It is estimated that about 7 million people die every year worldwide because of exposure to fine particulate pollution < 2.5 μm (PM_2.5), while about 91% of the world’s population live in areas with PM_2.5 concentrations above the allowable limits of 10–20 μg m⁻³ [4]. Southeast Asia (SEA) is not an exception to high pollution levels and experiences persistent haze conditions, especially during the dry season (June to September) due to extensive forest, agricultural and peat fires [5,6,7,8,9]. Malaysia is located in the main pathway of the SEA pollution outflow [10] that escalates the pollution levels due to trans-boundary aerosol transport. Air quality in Malaysia is considered rather degraded, as the annual mean PM_2.5 concentration is about 20 μg m⁻³ in 2019 [11], thus exceeding the limit set by WHO. Local sources such as traffic and industrial emissions, as well as biomass burning significantly contribute to the local pollution and aerosol loading [12,13,14]. According to the Department of Statistics Malaysia [15], chronic lower respiratory diseases (CLRD) is the fifth leading cause of death in Malaysia that was increased from 2.0% to 2.6% between the years 2017 and 2018. In fact, chronic exposure to PM_2.5 has significant impacts on human health by causing asthma, chronic obstructive pulmonary disease, lung cancer, cardiovascular and neurotoxic effects [16,17]. Therefore, continuous monitoring of the levels and mapping of spatial distribution of PM_2.5 is especially important for taking appropriate actions to maintain a good air quality over Malaysia [18,19].

Nowadays, many satellites provide reliable AOD (aerosol optical depth) products that can be used through various techniques for PM estimations from space, such as the Advanced Very High Resolution Radiometer (AVHRR) [20,21], Multiangle Imaging Spectroradiometer (MISR) [22], Medium Resolution Imaging Spectrometer (MERIS) [23], Spinning Enhanced Visible and Infrared Imager (SEVIRI) [24] and Moderate Resolution Imaging Spectroradiometer (MODIS) [25,26,27,28]. Among the satellite data products, AOD provided by MODIS sensor on board Terra and Aqua satellites has been widely used for many applications due to its high retrieval accuracy over land and near daily global coverage at 10 km and 3 km spatial resolution [26,27,29]. However, its usability in SEA is still limited due to large missing data series as a result of heavy and extended cloudiness [18,30,31,32]. On the other hand, geostationary satellites provide high temporal resolution (~15 min) data, thus limiting the problem of cloudiness. In the past few years, geostationary satellites such as Geostationary Operational Environmental Satellite (GOES) [33], Geostationary Ocean Color Imager (GOCI) [34], Fengyun-4 [35] and Himawari-8 [36,37,38] are available and provide continuous data over the SEA region.

Significant progress has been made in developing and establishing various techniques for estimating PM concentrations from space at local, regional and global scales. The progress in usage of linear and multi-linear statistics, regression-based, machine learning and hybrid models for estimations of PM_2.5 and PM₁₀ concentrations during the last decades is reviewed in recent works [1,39,40]. In the early 2000s, most researchers predicted PM using only AOD, by means of simple linear regression techniques [41,42]. Later on, more advanced techniques were developed to incorporate AOD and other important parameters that may influence PM distribution spatially and temporarily, starting from multiple linear regressions [18,43,44,45,46], chemical transport models (CTM) [47,48,49,50], mixed effect models (MEM) [51,52,53,54], artificial neural networks (ANN) [18,55,56,57,58], geographic weighted regression (GWR) [59,60,61,62] and generalized additive models (GAM) [63,64,65,66]. These techniques were used to capture the non-linear relationships that exist between the variables. Consequently, complex techniques have been developed by combining two or more statistical techniques; for instance, merging MEM and GWR [67] or incorporating MEM into GAM [68]. Nowadays, Machine Learning (ML) techniques such as deep neural network (DNN), support vector regression (SVR) and random forest (RF) enable to capture the complex relationships between parameters, exhibiting greater performance in estimating PM_2.5 [69,70] and are increasingly used in air quality studies [71,72,73,74]. Furthermore, some studies have incorporated meteorological factors with land use variables to predict the spatial and temporal variation of aeolian erosion and PM [75,76,77,78,79]. However, only few studies have incorporated air pollutant concentrations for PM_2.5 estimations [65,80,81]. Song et al. [65] explored the use of generalized additive model to estimate PM_2.5 concentrations in the Xi’an City, China (3581 square kilometres) using a combination of air pollutants (SO₂, CO, NO₂, and O₃), AOD and meteorological variables. The model was found to explain ~70% of the variance in PM_2.5 concentrations, with CO concentration and AOD represented most of the variation. The influence of air pollutants on the seasonal variability of PM_2.5 in an urban-industrial environment in Malaysia was investigated by [81] for one year. This study concluded that only gases (CO, NO₂, NO and SO₂) significantly affected the PM_2.5 mass, but not the meteorological factors (rainfall, wind speed and wind direction). Based on the prediction model of [80], temperature, Normalized Difference Vegetation Index (NDVI), humidity and residential area were found to be important parameters affecting the spatial variation of PM_2.5 in Jakarta, Indonesia. The same study [80], also showed that several parameters (PM₁₀, NO₂, SO₂, UV, rainfall, land use and NDVI) influenced the distribution of PM_2.5 in Taipei, Taiwan. More studies are therefore needed to find out the role of gases and meteorology in affecting the spatial and seasonal patterns of PM_2.5, even using meteorological normalization techniques in order to exclude the effect of changing meteorology on PM concentrations trends [82]. Recently, ML approaches (random forest regression models) were implemented to predict the large reductions in air pollutants, i.e., PM₁₀, NO₂, O₃, during the COVID-19 lockdown period [83]. Furthermore, PM₁₀, NO₂ and carbonaceous aerosols (organic carbon, elemental carbon) were also used in ML techniques (Lasso, Random Forest, AdaBoost, Support Vector Machine and Partials Least squares) for analysing air pollutants at street canyons [84]. Studies dealing with PM estimations in Malaysia are rather limited [18,30]. Shaziayani et al. [85] has reviewed PM₁₀ modelling studies in Malaysia, and only four studies used ML techniques in predicting PM₁₀. On the other hand, PM_2.5 studies are even fewer and most of them in Malaysia have been performed at small spatial scales [86,87].

The current study is the first in Malaysia and one of the very few works conducted worldwide aiming to estimate the PM_2.5 concentrations at a large (national) scale using pollution gases, AOD and meteorological factors based on machine-learning techniques. In order to extend the spatial coverage to the whole country (both Peninsular and Island Malaysia) and aiming to improve the accuracy of PM_2.5 estimates, this study integrates hourly AOD products from Himawari-8 satellite sensor, along with meteorological parameters and gaseous pollutants using machine learning techniques, i.e., random forest (RF) and Support Vector Regression (SVR). The models were developed separately for urban/industrial, suburban/rural sites and for the four dominant seasons in order to better represent the spatial (between sites) and temporal (between seasons) variation of the PM_2.5 concentrations. Variable importance analysis was conducted to identify the primary parameters that affect the PM_2.5 concentrations in Malaysia in a way to develop regression models for estimations of PM_2.5. The performance of RF and SVR was evaluated at different seasons and locations against measured PM_2.5 concentrations. The results will assist in representing the spatial and temporal evolution of PM_2.5 in Malaysia and for establishing measures in a way to improve air quality across the country.

Section 2 briefly describes the study area; Section 3 refers to the dataset that was used as variables in the prediction models; Section 4 refers to the PM_2.5 measurements across Malaysia. The results and model evaluation are included in Section 5, while Section 6 summarises the conclusions.

2. Study Area

Malaysia is one of the developed countries in SEA region with a rapid urbanization rate since 1970 [88]. Consequently, air pollution has become one of the serious environmental and human health concerns across the country [5,89], particularly in urban, industrialised and congested traffic areas such as Klang Valley (in the west coast of Peninsular Malaysia), Johor Bahru (southern tip of the Peninsula) and Georgetown, Penang (north of Peninsula) [90,91]. Air quality deteriorates at several parts of Peninsular Malaysia and in Borneo Island during the dry season mainly due to trans-boundary haze from neighbouring countries and regional/local forest fires. The concentrations of aerosols and air pollutants display a distinct seasonality, influenced by local meteorological conditions, i.e., rainfall, wind speed, relative humidity (RH) and temperature [86], being lower during the monsoon rainy season (November–March). In this study, 65 air quality monitoring stations (Figure 1) distributed across the Peninsular and Island Malaysia (Labuan, Sabah and Sarawak) were used to analyse the air pollution levels (PM_2.5, SO₂, NO₂, CO and O₃ concentrations). The stations are representative of industrial (7 stations), urban (10 stations), suburban (36 stations) and rural (12 stations) areas.

3. Dataset

The dataset used in this study consist of PM_2.5 and air pollutant concentrations, along with meteorological parameters from ground stations, as well as AOD data from Himawari-8 satellite.

3.1. Ground Measurements

PM_2.5 concentrations in Malaysia have been measured since April 2017. The Department of Environment, Malaysia (DOE) increased the number of air quality monitoring stations across the country from 52 to 65 stations in 2017 (Figure 1). Furthermore, these stations also measure meteorological parameters (e.g., ambient temperature, TEMP; RH, wind speed, WS; wind direction, WD) and gaseous pollutants (e.g., nitrogen dioxide, NO₂; carbon monoxide, CO; sulphur dioxide, SO₂; ozone, O₃). The stations are strategically distributed to represent urban, industrial, suburban and rural areas [92]. PM_2.5 measurements were performed via the TEOM 1405DF, which is a continuous dichotomous ambient air monitoring system with two Filter Dynamics Measurements Systems [93], able to measure PM_2.5 and PM₁₀. SO₂, NO_2, CO and O₃ are measured using Thermo Scientific model 43i, model 42i, model 48i and model 49i, respectively [93,94]. RH and TEMP were recorded using a Climatronic AIO 2 Weather Sensor (Climatronic Corporation) [95]. All the ground data were obtained on an hourly basis covering the period from January 2018 to December 2019.

All air quality and meteorological measurements went through quality assurance and quality control (QA/QC) procedures. Instruments for the detection of gases were manually calibrated once a fortnight. Flow verification for PM₁₀ and PM_2.5 measurements using TEOM was conducted once a month. The data removal during the QC check was predominantly due to insufficient measurements and instrument failure, while some perturbed data were also excluded as outliers in a second-level QC check [95].

3.2. Satellite Data

Himawari-8 is a geostationary satellite operated by the Japan Meteorological Agency. It was launched on 7 October 2014 and carries the Advanced Himawari Imager (AHI) sensor, which is equipped with 16 bands from visible to infrared [36]. Himawari-8 releases AOD products at two levels, namely Level 2 (10 min temporal) and Level 3 (hourly and daily), which have been used for various applications including estimation of PM [96,97,98,99], dust detection [100] and aerosol data assimilation [101]. The L3 product is an improved version of the L2 AOD product that minimized cloud contamination [102] and has a 5 km spatial resolution. Himawari-8 AOD at 500 nm is associated with quality assurance levels namely “very good”, “good”, “marginal” and “no confident (or no retrieval)” [99]. In this study, only the “very good” L3 AOD₅₀₀ retrievals were considered for PM_2.5 estimations, downloaded from the Japan Aerospace Exploration Agency (JAXA) website: available online: http://www.eorc.jaxa.jp/ptree/index.html (accessed on 10 May 2021) for the period January 2018–December 2019. Recently, Himawari-8 AODs were used to estimate the PM_2.5 concentrations over Hubei province, China [103]. Application of the Himawari-8 L2 AOD data over Malaysia revealed an overestimation by 24.2% [104], while the L3 AOD products displayed a better agreement with the Aerosol Robotic Network (AERONET) AODs with a coefficient of determination R² = 0.81, root mean square error (RMSE) of 0.13 and an overall overestimation of only 1% [92]. Moreover, Himawari-8 AODs presented a good agreement with AERONET AODs in China (R² = 0.41–0.83; RMSE = 0.18–0.31) [105], Southeast Asia (R² = 0.64; RMSE = 0.28), East Asia (R² = 0.83; RMSE = 0.14) [106], Korea (R² = 0.69; RMSE = 0.19) and Beijing-Tianjin-Hebei (R² = 0.76; RMSE = 0.36) [107].

4. PM_2.5 Estimation

The overall methodology used for the PM_2.5 estimations over Malaysia is illustrated in Figure 2. The model inputs consist of hourly AOD, SO₂, NO₂, CO, O₃, WS, WD, TEMP and RH values. Hourly AOD data from Himawari-8 were extracted at 5 × 5 km over the air quality monitoring stations and temporally collocated with ground measurements. Wind direction was used in the model because wind blowing from a highly polluted area can influence air quality in other downwind places. Wind speed enables to accelerate pollutants travelling from other places but also contributes to the dilution processes at local level [108]. On the other hand, temperature can trigger biogenic emissions, photochemical reactions and secondary aerosol formation over the region [109] and also control the temperature inversions, which can trap the pollutants near the surface [45,108,110]. Finally, RH may affect the hygroscopic growth of particles and enhance the aerosol scattering [111,112,113].

In this study, we utilized and evaluated two ML models, namely, Support Vector Regression (SVR) and Random Forest (RF) to estimate PM_2.5 concentrations at the 65 air quality monitoring stations across Malaysia. The input variables to the SVR and RF models are selected based on our previous study [18,31] and other literature [114,115]. In total, we developed 7 different models to represent the spatial and seasonal effects on PM_2.5 distributions. Model 1 considers all data from the 65 stations, but other models, i.e., models 2 and 3, only represent urban/industrial and suburban sites, respectively, while models 4 to 7 represent different seasons (wet, dry and two inter monsoon). The models are described in the subsequent sections.

4.1. Machine Learning (ML) for PM_2.5 Estimation

Nowadays, several studies have used machine learning (ML) techniques, including RF model, aiming to increase the accuracy in prediction of PM concentrations, since these models are flexible in nonlinear approaches [72,116,117,118,119]. The SVR and RF techniques were particularly selected in this study to achieve more accurate PM_2.5 estimations using satellite-derived AOD, meteorological parameters and gaseous pollutants as predictor variables. The Classification And REgression Training (CARET) package was used to perform the RF and SVR modelling. The data splitting, pre-processing, model tuning and variable importance analysis were executed in a R environment. SVR depends on the kernel function and due to its excellent generalization capability, it is able to minimize the overfitting [120], and therefore, it has been used for PM estimations [69,108,121]. SVR can fit the errors within a certain threshold by finding an appropriate boundary line (between hyperplane) to suit the data. The flexibility of SVR depends on the selection of the parameter such as kernel function, cost function and epsilon value. There are four types of kernel functions namely linear, polynomial, sigmoid and radial basis function (RBF) that were used in this study for capturing the non-linear dynamics [69], whilst cost function was used to avoid any overfitting of the data, as small cost value leads to large margin (or wide boundary line) and causes overfitting in the model. The epsilon value controls the number of support vectors used to develop the regression function, while the smaller epsilon value indicates an optimum accuracy. Initially, the SVR parameters are selected based on trial-and-error values, but we found that the default values (Supplementary Materials Table S1) included in the R package “e1071” provided the most promising results.

Random Forest (RF) is a tree-based ML technique proposed by [122]. Theoretically, RF model is an ensemble of multiple decision trees and uses the majority vote/decision of the trees as the final RF model [123,124]. The algorithm becomes more robust when more decision trees are constructed. RF randomly selects parameters in order to develop each tree, and therefore, it reflects potentially complex effects of predictors on the prediction [125]. The purpose of selecting random predictors instead of all predictors is to reduce the correlation between trees in order to make them disparate [126]. Thus, the variance of the RF prediction can reduce any overfittings. The number of decision trees can be modified to reduce the training time according to a required accuracy and computing capability [127]. In this study, RF model was run using the “Random Forest” in the R package. Since, RF is a non-parametric algorithm, here we only set the two most important parameters—although RF can have more parameters—which are mtry and ntree. Parameter mtry is a number of predictors sampled for splitting at each node while ntree is a number of trees in the forest. If mtry value is too small, it might be none of significant parameters included in the subset, and the insignificant parameters would be selected for a split. Therefore, the trees have poor predictive ability [126]. In this study, we set the mtry = 3 (as default: mtry = p/3, where p is the number of parameters used in the model), while ntree is set as 500 in the model. For tuning, we only tune mtry because the CARET package has automatic tuning for mtry only. Therefore, in this study, for ntree, we used the default value [128,129]. The results were obtained based on best mtry tuning accuracy. It should be mentioned that a limitation of this study is that models were not broadly optimized.

4.2. Model Validation

The total number of the matching dataset, covering all parameters at the 65 stations in Malaysia from 2018 to 2019, is 13,376. The matching data were randomly partitioned at a fraction of 70% for model calibration (model development) and 30% for model validation. In the model development, a sample based 10-fold cross validation technique was used, where the calibration data were randomly divided into 10 subsets; at any single moment, one subset was used for validation and the remaining subsets were used for calibrating the model. The average value of the results of the 10 subsets was adopted as the model accuracy. The sample based 10 cross validation (CV) performed validation with matchup sample from both spatial and temporal dimension. This is a commonly used CV-based technique to reveal the overall predictive ability of PM_2.5 estimation models [130]. Then, the final models were validated using the 30% of the remaining data. Statistical indicators such as the coefficient of determination (R²), Root Mean Square Error (RMSE), mean bias error (MBE) and Nash-Sutcliffe Efficiency (NSE) were used to evaluate the accuracy of the models. The NSE is a normalized statistic, which can determine the magnitude of the residual variance to the measured data variance and indicates how well the measured PM_2.5 versus estimated PM_2.5 data fits the 1:1 line (best fit line).

$$RMSE=\sqrt{\frac{{\displaystyle \sum _{i=1}^N{\left(I_m-I_c\right)}^2}}{N_t}}$$

(1)

$$MBE\hspace{0.17em}=\frac{{\displaystyle \sum _{i=1}^N\left(I_m-I_c\right)}}{N_t}$$

(2)

$$NSE=1-\frac{{\displaystyle \sum _{i=1}^N{\left(I_m-I_c\right)}^2}}{{\displaystyle \sum _{i=1}^N{\left(I_m-\overline{I_m}\right)}^2}}$$

(3)

where I_m and I_c are the measured and estimated PM_2.5 concentrations, respectively, $\overline{I_m}$ is the average of the measured PM_2.5 and N the total number of measurements. NSE = 1 indicates an ideal model performance, while NSE = 0 shows model predictions as accurate as the mean of the observed data. Lower RMSE and MBE values correspond to better performance and to lower biases from the ML models.

4.3. Variable Importance

Variable importance statistic was used to analyse the contribution of each variable in PM_2.5 estimations. Since SVR is a kernel-based model, and we do not know the concrete form of its nonlinear mapping function, and the weight vector (ω) cannot be computed directly [128], it is complicated to analyse the variables importance statistic. On the other hand, there are two famous measures for RF, which are mean decrease accuracy (MDA) and mean decrease Gini (MDG). The MDG is based on Gini importance which measures the average gain by splits of a given variable, whilst MDA is based on out of bag (OOB) samples. In RF model, each tree is grown based on a bootstrap sample of the training data, and those data that were not used in the bootstrap sample are known as out of bag (OOB) samples [126]. The MDA measures the accuracy of the model losses by permuting each variable. This technique is considered as most efficient variable importance for random forest [129,131], and it was preferred in this study as less bias compared to Gini importance. The higher percentage value of the variable importance indicates higher influence of the corresponding variable to PM_2.5 estimations. In R script, we used the “varImp” function, which can automatically scale the importance scores in values between 0 and 100.

5. Results and Discussion

5.1. Descriptive Statistics

The descriptive statistics of the measured variables that are used in SVR and RF models for all stations in Malaysia are summarized in

Table 1. Statistical values for the measured parameters in all air-pollution monitoring sites.

	PM2.5 (µg m−3)	AOD	SO2 (ppb)	NO2 (ppb)	O3 (ppb)	CO (ppm)	WS (ms−1)	RH (%)	TEMP (°C)
Mean	21.86	0.69	1.2	5.23	25.2	0.60	1.73	66.58	30.70
Median	17.07	0.46	1.0	3.81	27.3	0.56	1.60	66.85	30.88
Stdev	19.15	0.68	0.9	6.1	15.1	0.28	1.02	10.35	2.31

. The columnar AOD₅₀₀ over the Malaysian sites during 2018–2019 exhibits a mean of 0.69, which is above the median value 0.46 due to many episodic aerosol events with AODs above 2, representing thick smoke plumes from extensive fires in Indonesia and Indochina [10,132,133]. The measured PM_2.5 concentrations at the 65 examined sites follow a similar distribution with a higher mean (21.9 µg m⁻³) than median (17.1 µg m⁻³) and a maximum value of 230 µg m⁻³ (Table 1, Figure 3). These PM_2.5 levels are similar to those reported at several sites in Southeast Asia [86,134]. NO₂ and CO exhibit means of 5.2 ppb and 0.6 ppm, respectively, while tropospheric O₃ levels (25.2 ppb) are considered rather high with deleterious effects on human health [135,136].

The box-whisker plots for AOD₅₀₀ and air pollutants, separated for urban/industrial and suburban/rural sites and for the four seasons, are shown in Figure 3. The analysis showed that the seasonal-mean PM_2.5 at the urban/industrial sites displayed higher values in all seasons. During the dry season, the maximum PM_2.5 levels were found to be 230.3 μg m⁻³ for the urban/industrial and 226.3 μg m⁻³ for the suburban sites, with means of 31.26 µg m⁻³ and 26.38 µg m⁻³, respectively. This is attributed to the prevailing southwest wind carrying biomass-burning aerosols from Indonesia due to extensive forest fires in this season [90]. However, the seasonal mean PM_2.5 concentrations do not notably differ in the other seasons (lying between 17.80 µg m⁻³ and 22.33 µg m⁻³ for both urban/industrial and suburban/rural sites), indicating a year-long PM_2.5 laden atmosphere across Malaysia. Regarding the Himawari-8 AOD, it is higher during the dry season, while slightly lower mean values (~0.6 to 0.8) were observed in the other seasons, with marginal differences between urban/industrial and suburban/rural sites (Figure 3b). The seasonal and site variations of the columnar AODs are in agreement with ground PM_2.5 concentrations, indicating that the industrial and traffic emissions are the main pollution sources for urban centres and surrounding areas. In general, heavy precipitation during the wet season only marginally reduced the aerosol levels since severe pollution episodes with PM_2.5 > 100 µg m⁻³ and AODs > 3 were also present. However, the columnar AOD displayed a very low correlation (R² = 0.09) with the surface PM_2.5, indicating (i) a significant aerosol loading aloft and (ii) different sources and temporal variability between surface PM_2.5 and AODs [137]. In addition, the mean concentrations of CO, NO₂ and O₃ (Figure 3c–e) are higher over the urban/industrial areas compared to suburban/rural sites in each season. Motor vehicle and power plants emissions are the major contributors to CO and NO₂ concentrations in Malaysia with about 95.7% and 66%, respectively [11], thus explaining the higher NO₂ and CO levels in urban/industrial areas. The stronger correlation was found between CO and PM_2.5 (R² = 0.33), revealing that the particulate pollution in Malaysia is mostly related to local sources of fossil fuel and biofuel combustion, which enhance CO emissions [86,138]. NO₂, which is mostly related to vehicular emissions, was negligibly associated with PM_2.5 concentrations (R² = 0.1), The overall mean SO₂ concentration was found to be 1.2 ppb, with slightly larger levels in the dry season, exhibiting means of 1.5 ppb and 1.3 ppb for the urban/industrial and suburban sites, respectively.

5.2. Models for PM_2.5 Estimation

In this study, seven models were developed for PM_2.5 estimations in Malaysia using ML techniques, namely, SVR and RF. These models were developed in order to better capture the remarkable spatial (between stations of different characteristics) and temporal (between seasons) variations of PM_2.5 and to examine the model’s capability in representing the levels and evolution of PM_2.5. The developed models are:

• Model 1: Overall model;
• Model 2: Spatial model (urban/industrial);
• Model 3: Spatial model (suburban/rural);
• Model 4: Temporal model (dry season);
• Model 5: Temporal model (wet season);
• Model 6: Temporal model (inter-monsoon, April–May);
• Model 7: Temporal model (inter-monsoon, October).

The models’ inputs consist of AOD, SO₂, NO₂, CO, O₃, WS, WD, TEMP and RH, with the total number of the matching samples for all variables (overall model) to be 13,376 at the 65 monitoring stations from January 2018 to December 2019. The scatter plots between measured and predicted PM_2.5 concentrations via the SVR and RF models are shown in Figure 4 and Figure 5 for the validation of the 7 developed models, while Supplementary Materials Table S2 summarizes the statistical indicators of R², RMSE, MBE and NSE, as well as the importance ranking, for each input variable at the 7 developed models. The statistical indicators correspond to the averaged values of each model for the given number of data (N).

The evaluation of the calibration datasets for the overall model 1 showed that the SVR model yielded comparable accuracy (R² = 0.69; RMSE = 10.62 μg m⁻³, NSE = 0.679 and MBE = −1.392) to RF (R² = 0.66; RMSE = 11.28 μg m⁻³, NSE = 0.939 and low MBE = 0.066) (Supplementary Materials Table S2). The validation dataset displayed also small differences between the two models and comparable statistics with the calibration datasets i.e., R² = 0.66 and RMSE = 12.11 μg m⁻³ for SVR and R² = 0.62, RMSE = 11.40 μg m⁻³ for RF (Figure 4).

Furthermore, other models for estimating PM_2.5 were also developed by splitting the entire datasets initially into two categories namely urban/industrial (3357 data points) and suburban/rural (N = 9798). The RF model performed slightly better compared to SVR for urban/industrial (model 2) and suburban (model 3) sites, with R² = 0.76, RMSE = 11.47 μg m⁻³, NSE = 0.735 and R² = 0.67, RMSE = 12.47 μg m⁻³, NSE = 0.661, respectively. Regarding the calibration datasets, the NSE values of all the RF models were usually above 0.9, indicating an accurate model’s performance, compared to the SVR models, which exhibited NSE values in the range of 0.55 to 0.79 for the calibration datasets (Supplementary Materials Table S2). The model validation revealed that RF models performed slightly better than the SVR models for both locations. The validation of SVR and RF models for suburban sites was considered satisfactory with R² = 0.61, RMSE = 11.53 μg m⁻³ and R² = 0.64, RMSE= 10.76 μg m⁻³, respectively (Figure 4b,c). Both models performed better at urban/industrial sites compared to the suburban sites, while all models underestimated the large PM_2.5 concentrations (PM_2.5 > 60 μg m⁻³). For PM_2.5 below 50 µg m⁻³, where the vast majority of the data points lie, the underestimated and overestimated data points are almost equal for all the models and the regression line tends to coincide with the 1-1 line (Figure 4).

In addition, four temporal models were also developed for each season in Malaysia, namely, dry season (June–September; Figure 5a,b), wet season (November–March; Figure 5c,d) and two inter-monsoon seasons (April–May; Figure 5e,f and October; Figure 5g,h). The total number of data for the temporal models is 4223, 5797, 2441 and 915 for the dry season, wet season, April-May and October, respectively. These datasets were randomly divided into calibration and validation groups, and the results of the statistical evaluations of the seasonal PM_2.5 predictions are included in Supplementary Materials Table S2 (models 4–7). The estimation accuracy of PM_2.5 concentrations via the SVR and RF models varied between the seasons, while it was relatively lower in the wet and inter-monsoon (April-May) seasons for both models. For instance, the R² values for the calibration dataset of the SVR model were 0.81, 0.62, 0.58 and 0.74 for the dry, wet seasons, April-May and October, respectively. The RF models exhibited slightly lower RMSE values compared to the SVR models, as 13.61 μg m⁻³, 9.99 μg m⁻³, 8.36 μg m⁻³ and 12.41 μg m⁻³ against 14.0, 10.9, 9.38 and 12.46 μg m⁻³ for dry, wet, April-May and October, respectively. Furthermore, RF models displayed higher NSE values compared to SVR, while the PM_2.5 underestimations from the SVR models (negative MBE values) are eliminated by using the RF models (Supplementary Materials Table S2; Figure 5).

The statistical evaluators of the developed models in Malaysia are mostly comparable to those found from multi-variate models including AOD and several meteorological parameters (Temp, RH, WS, Dew point, mixing height) for PM_2.5 estimations in Indian cities [139]. More recently, [39] developed a deep neural network consisted of recurrent layers for extracting the relationship between high-resolution (1 km) MODIS observations and PM₁₀, PM_2.5 concentrations in Tehran, Iran. The PM_2.5 and PM₁₀ estimations resulted in RMSE values of 11.66 μg m⁻³ and 23.79 μg m⁻³, respectively, comparable to the current results and previous RMSE values for PM₁₀ estimations in Malaysia (11.61 μg m⁻³; [18]) and Delhi (18.99 μg m⁻³; [140]). MODIS-MAIAC AODs and columnar water vapor (CWV), along with meteorological parameters and land-use data, were included in a linear mixed effect model (LME) and a RF model for daily PM_2.5 estimations at high spatial resolution (1 km × 1 km) over the Indo-Gangetic Plains (IGP), India [118]. The RF model exhibited higher accuracy with R² = 0.87 and relative RMSE of 24.5%, compared to LME [118]. The spatial distributions of the R² (~0.6 to 0.9) and RMSE (~20 to 40 μg m⁻³) values from PM_2.5 estimations across the IGP [118] were mostly comparable to those observed over Malaysia.

The hourly time series of the measured and predicted PM_2.5 concentrations via the SVR model, separately for the station characteristics and seasons and for the overall model, are shown in Figure 6. The results verify the good performance of the SVR model in predicting the PM_2.5 concentrations across Malaysia—RF performed slightly better with very similar results—also revealing an underestimation at the highest PM_2.5 values. However, the model’s underestimation in representing PM_2.5 peaks is not systematic, and in many cases, models reproduce satisfactorily the high PM_2.5 concentrations, even overestimating them (Figure 6).

The residual analysis for the model validation datasets revealed that the frequency of residuals approached the normal distribution peaking around zero for all the models (Supplementary Materials Figure S1). However, the frequency distributions were slightly shifted towards negative values, as the highest model underestimations may reach to −60%, but for very rare cases. In the vast majority of the cases, the predictions were quite accurate, indicating that the used ML models are satisfactory for PM_2.5 estimations in Malaysia.

Atmosphere is a complex system and composed by various substances like air molecules and solid and liquid particles of various sizes, shapes and chemical composition [141,142]. Therefore, combining many auxiliary data such as meteorological factors, aerosols, gases and land use allows for a better estimation of PM_2.5. To date, most PM prediction studies found that inclusion of meteorological factors has improved the PM estimations, because each meteorology parameter may modulate the PM concentrations in a different way [43,45,71]. Nowadays, ML techniques and RF models are frequently used in estimating PM concentrations at many regions around the world [39,116,119]. The statistical indicators from the models’ calibration and validation in this study (Supplementary Materials Table S2; Figure 4 and Figure 5) are mostly comparable to those presented in other studies using various methods and ML approaches for estimates of PM concentrations around the world (

Table 2. Results from previous studies using machine learning techniques for PM estimations from space.

Author	Study Area	Input Data			Method	R/R2	Accuracy
		Source of AOD	Other Parameters	Output
[143]	China	MODIS AOD 10 km (Terra and Aqua)	RH, AT, WS, SP, PBLH, NDVI, population and road data	PM2.5	Geoi-DBN	Sample based CV R2 = 0.88 Site based CV R2 = 0.82	Sample based CV RMSE = 13.03 μg m−3 Site based CV RMSE = 16.42 μg m−3
[146]	China	MODIS AOD 3 km (Terra and Aqua)	Lat, long, month, RH, AT, WS, SP, PBLH	PM2.5	GRNN	R2 = 0.89	RMSE = 16.51 μg m−3
[147]	China	MAIAC AOD 1 km	AT, AP, evaporation, precipitation, RH, sunshine duration and WS	PM2.5	GW-GBM	Exclude missing AOD R2 = 0.74 Include missing AOD R2 = 0.76	Exclude missing AOD RMSE = 24.3 μg m−3 Include missing AOD RMSE = 23.0 μg m−3
[128]	Cincinnati, OH, USA	MODIS AOD 3 km (Terra and Aqua)	Visibility, PBLH, TEMP, RH, total and rate precipitation, P, WS, WD, land cover, roadways, green space, spatiotemporal convolution layer	PM2.5	RF	Overall R2 = 0.90 Spatial R2 = 0.87 Temporal R2 = 0.84	Overall RMSE = 2.45 μg m−3 Spatial RMSE = 2.83 μg m−3 Temporal RMSE = 3.13 μg m−3
[148]	British Columbia, Canada	MODIS AOD 3 km (Terra)	LST, humidity, vapour, NDVI, albedo from MODIS product. PBLH, WS. Elevation from SRTM	PM2.5	MLR BRNN SVM LASSO MARS RF XGBoost Cubist	R2 = 0.22 R2 = 0.31 R2 = 0.30 R2 = 0.24 R2 = 0.31 R2 = 0.49 R2 = 0.46 R2 = 0.48	RMSE = 3.24 μg m−3 RMSE = 3.04 μg m−3 RMSE = 3.13 μg m−3 RMSE = 3.20 μg m−3 RMSE = 3.05 μg m−3 RMSE = 2.67 μg m−3 RMSE = 2.71 μg m−3 RMSE = 2.64 μg m−3
[69]	BTH, China	MODIS AOD 10 km (Terra and Aqua)	PBLH, TEMP, SLP, humidity, WD and WS	PM2.5	OR Rpart RF SVM	R = 0.73–0.76 R = 0.68–0.83 R = 0.69–0.84 R = 0.77–0.88	RMSE = 36.92–42.48 μg m−3 RMSE = 35.42–46.20 μg m−3 RMSE = 36.34–44.59 μg m−3 RMSE = 29.50–38.32 μg m−3
[115]	East coast peninsular Malaysia	-	AT, RH, WS, GR, MSLP, rainfall, CO, NO2, and SO2	PM10	MLR MLP RBF	R2 = 0.594–0.706 R2 = 0.691–0.794 R2 = 0.827–0.929	VIF = 1.077–1.926 RMSE = 8.49–9.57 μg m−3 RMSE = 9.19–4.08 μg m−3
[141]	BTH, China	MODIS AOD 10 km (Aqua)	AT, RH, WS, WD and P	PM2.5	MLR MARS SVR RSRF	R2 = 0.733 R2 = 0.776 R2 = 0.850 R2 = 0.843	RMSE = 33.016 μg m−3 RMSE = 30.180 μg m−3 RMSE = 24.745 μg m−3 RMSE = 25.320 μg m−3
[149]	Wuhan, China	Himawari-8 AOD L3	MODIS NDVI, RH, AT, WS, SP, PBLH, DEM	PM2.5	DL	R2 = 0.850	RMSE = 9.303 μg m−3
[144]	China	MAIAC AOD 1 km	TEMP, total precipitation, evaporation, PBLH, RH, SP, WS, WD, MODIS Land use Cover, NDVI, DEM	PM2.5	RF STRF	R2 = 0.98 R2 = 0.98	RMSE = 6.40 μg m−3 RMSE = 5.57 μg m−3
[150]	Shenzhen, China	MAIAC AOD 1 km	EWS and RH	PM2.5	RF IRF	R2 = 0.88 R2 = 0.91	RMSE = 4.3 μg m−3 RMSE = 3.66 μg m−3
[145]	East Asia (Eastern China, Korean Peninsular and Japan)	GOCI, GEOS-Chem	NDVI, urban ratio, DEM, precipitation, AT, ST, dew point temperature, RH, max WS, visibility, PBLH, SP, solar radiation, road density, population density	PM10 PM2.5	RF	R2 = 0.88 R2 = 0.90	RMSE = 26.9 μg m−3 RMSE = 15.77 μg m−3
[70]	Greater London	MAIAC AOD 1 km	Population density, cloudiness, barometric pressure, WD, WS, dew point temperature, land use variable (type, distance to water, airport, PBLH, NDVI, traffic count, elevation etc)	PM2.5	GBM RF Deep NN KNN ensemble model	Overall model R2 = 0.826 R2 = 0.830 R2 = 0.793 R2 = 0.791 R2 = 0.828	Overall model RMSE = 4.331 μg m−3 RMSE = 4.278 μg m−3 RMSE = 4.728 μg m−3 RMSE = 4.721 μg m−3 RMSE = 4.231 μg m−3
[151]	Guwahati, India	−	CO, NO2, SO2, AT, RH, WS, rainfall	PM10	MLR MLP CART	R2 = 0.61–0.68 R2 = 0.64–0.69 R2 = 0.52–0.63	RMSE =29.31–31.99 μg m−3 RMSE =31.02–31.74 μg m−3 RMSE = 39.98–41.24 μg m−3

List of abbreviations: 1. Study area: BTH (Beijing-Tianjin-Hebei). 2. Sensor: GOCI (Geostationary Ocean Color Imager), MAIAC (Multiangle Implementation of Atmospheric Correction), MODIS (Moderate Resolution Imaging Spectroradiometer), SRTM (Shuttle Radar Topography Mission), 3. Parameter: AT (air temperature), CO (carbon monoxide), EWS (extreme wind speed), GR (global radiation), Lat (Latitude), Long (longitude), LST (land surface temperature), MSLP (mean sea level pressure), NDVI (Normalized difference vegetation index), NO₂ (nitrogen dioxide), P (pressure), PBLH (planetary boundary layer height), RH (relative humidity), SF (surface temperature), SLP (sea level pressure), SO₂ (sulphur dioxide), SP (surface pressure), ST (surface temperature), TEMP (temperature), WD (wind direction), WS (wind speed). 4. Method: BRNN (Bayesian Regularized Neural Networks), CART (Classification and Regression Trees), Cubist (rule based tree model), DL (Deep Learning), DNN (deep neural network), GBM (Gradient Boosting Machine), Geoi-DBN (Geo-intelligent Deep Belief Network), GRNN (generalized regression neural network), GW-GBM (Geographically- Weighted Gradient Boosting Machine), IRF (improved random forest), KNN (k-nearest neighbour), LASSO (Least Absolute Shrinkage and Selection Operator), MARS (Multivariate Adaptive Regression Splines), MLP (multilayer perceptron), MLR (Multiple Linear Regression), OR (orthogonal regression), RBF (radial basis function), RF (random forest), Rpart (regression tree), RSRF (hybrid remote sensing and random forest), STRF (space-time random forest), SVM (Support Vector Machines), SVR (support vector regression), XGBoost (eXtreme Gradient Boosting). 5. Accuracy: RMSE (root means square error), VIF (variance inflation factor).

). It should be noted that this study obtained reasonable results at national scale without including land use information compared to previous works [128,143,144,145]. Besides that, the validation techniques may be also different, as for instance [143] used sample and site based 10 CV in order to assess the spatial performance, whilst our study only used sample based 10 CV since it can be used to reflect the overall predictive ability [130].

Determining the strength of the correlation between PM_2.5 and all the parameters used for its prediction is very important because it can indirectly portray the pollution process and the source of the pollution. The results of the variable importance analysis for the RF model have been included in Supplementary Materials Table S2 and are shown in Figure 7. For the overall model 1, CO is the highest contributor to the PM_2.5 estimations, similar to the other models, as discussed above, and is followed by AOD, O₃, NO₂, SO₂ and the meteorological parameters. This is because both PM_2.5 and CO are originated from common sources in Malaysia like biomass burning and traffic which enhance the CO emissions [86,138]. Besides that, CO may stay in the atmosphere for a long period (weeks or months), being able to get transported in high concentrations from biomass-burning areas in Indonesia [152]. Parameters with the least importance in PM_2.5 estimations are RH, WD, TEMP and finally the WS with a zero score (Figure 7a). Similar to the overall model (model 1), the contributions of the meteorological parameters were relatively weak in the spatial models as well, which also revealed CO, AOD and O₃ as the most important variables (Figure 7b,c). CO remains the most important predictor in the seasonal models as well, with minimal contributions from the meteorological variables (Figure 7d–g). The WS and WD have minimal contributions, in agreement with [81], who found that both parameters were not significantly correlated with PM_2.5 in the Klang Valley region in Peninsular Malaysia. Although AOD usually exhibited a high correlation with PM_2.5 [153], in our case, there was not a direct association with PM_2.5 concentrations, implying complicated pollution conditions in the vertical layer over Malaysia [141]. Generally, PM_2.5 is affected more by gaseous pollutants and not so much by the columnar AOD (missing values due to cloud cover and elevated aerosol layers) and meteorological conditions since Malaysia has rather uniform weather conditions throughout the year. Therefore, influence of the meteorological parameters is minimal towards seasonally changing PM_2.5 concentrations. However, previous studies in Malaysia showed that the meteorological parameters affected the coarse particles, e.g., PM₁₀ concentrations [81], indicating a meteorological-dependent character of the coarse-mode aerosols. In our previous study [18], we found that the estimations of PM₁₀ concentrations based on satellite AODs were significantly improved after inclusion of the meteorological parameters in the model.

Similar to the current results, [154] also found that the parameters with the highest importance in predicting PM_2.5 concentrations are CO, NO₂, SO₂ and AOD. Inclusion of the pollutant gases improved the performance of their RF model from R² = 0.69, RMSE = 41.63 μg m⁻³ to R² = 0.81, RMSE = 32.74 μg m⁻³, in a similar way as in the present study. These results were also in agreement with [65] who found that CO was the most important variable that explained 20.65% of the variation in estimated PM_2.5 concentrations in Xi’an, China, exhibiting a strong correlation with AOD. Furthermore, [81] studied the PM_2.5 composition in Klang Valley, Malaysia, and concluded that CO, NO₂, NO and SO₂ mostly affected the PM_2.5 concentrations.

This is the first study conducted in Malaysia aiming to estimate the PM_2.5 concentrations based on machine-learning techniques. The satisfactory accuracy of the estimates, despite the biases and challenges in representing PM_2.5 pollution episodes, is especially important for the development of models aiming to systematically monitor PM_2.5 over the country, especially at remote areas with unavailability of measurements. However, the only 65 operational stations are still insufficient to cover the whole Malaysian territory with an area of 330,290 km². Establishing more air quality monitoring stations is very costly, and a certain station is only capable to satisfactorily represent the pollutant concentrations within a radius of about 15 km [155]. Alternatively, remote sensing data encourages more studies to be conducted on atmospheric particulates and air quality, since satellite technology considers AOD as a key predictor of PM over a large area [156]. This would also help in evaluating the influence of the local/regional emissions from anthropogenic activities against those attributed to natural causes or long-range transported aerosols, mostly smoke, from Indonesia and other parts in northern Indochina.

6. Conclusions

The current study developed new machine-learning models, namely, Random Forest (RF) and Support Vector Regression (SVR), to estimate PM_2.5 concentrations across Malaysia for the first time covering the years 2018 and 2019. Satellite (Himawari-8) AOD, ground-measured air pollutants (NO₂, CO, SO₂, O₃) and meteorological parameters (temperature, relative humidity, wind speed and direction) were used as input variables. Due to the high spatial (between stations with different characteristics like urban/industrial and suburban/rural) and temporal (between seasons) evolution of the PM_2.5 levels across Malaysia, seven sub-models were developed separately for the different sites (urban/industrial, suburban/rural) and seasons (dry, wet and two inter-monsoons (April-May and October)), and one overall model. Of the available dataset, 70% was randomly selected for the model calibration, and the remaining 30% for the model validation. The PM_2.5 predictions of each model are compared to those measured at 65 air pollution monitoring stations, using standard statistical estimators.

For the overall model, SVR calibration performed slightly better than RF with R² = 0.69 and RMSE = 10.62 μg m⁻³ against measured PM_2.5 concentrations. Whilst for the spatial models, the RF validation performed slightly better than SVR, with statistical indicators of R² = 0.76, RMSE = 11.47 μg m⁻³ for urban/industrial, and R² = 0.64, RMSE = 10.76 μg m⁻³ for the suburban/rural sites. Therefore, both RF and SVR models displayed slightly higher performance for PM_2.5 estimations at urban/industrial sites with higher levels of AOD and air pollutants. Furthermore, the estimation accuracy of SVR and RF models was lower in the wet (November–March) and inter-monsoon (April–May) seasons compared to the dry (June–September) season. Based on the model accuracy and variable importance analysis, CO was always the most influential predictor variable for PM_2.5 estimations in Malaysia, followed by AOD, O₃, NO₂, SO₂ and meteorological parameters but with different order depending on the dataset and model. An important finding was the very weak correlation and contribution of the meteorological variables to PM_2.5 estimations. Furthermore, very low correlation was found between PM_2.5 and columnar AOD, indicating that surface pollution followed a different temporal pattern than AOD and the presence of a significant aerosol layer aloft due to transported smoke plumes from wildfires in southeast Asia. The current results showed that the use of machine-learning techniques for PM_2.5 estimations over Malaysia was promising as these models can satisfactorily represent the values and temporal evolution of PM_2.5 concentrations over both urban/industrial and suburban/rural sites although the underestimation of the highest PM_2.5 levels. In a next step, gaseous pollutants from satellite remote sensing observations will be included in ML approaches in order to estimate PM_2.5 concentrations over large areas, aiming to cover the whole Malaysian territory.