This study uses the publicly available hourly PM_2.5 data consisting of 2 years of EPA and PA measurements from October 2019 to September 2021 [10, 19]. The hourly averages were then converted to 24-h averages for this analysis. We collected meteorological data from five nearby stations of the National Oceanic and Atmospheric Administration (NOAA), [20] with a distance of each nearest EPA, PA sensor, and NOAA station. The meteorological variables include temperature, relative humidity (RH), wind speed (WS), and wind direction (WD). The information on these EPA, PA sampling can be extracted from Figure 1 and from supplementary of [14], the information on NOAA sites can be extracted from Supplementary Table S1. For consistency and validation of results with [12, 14], the PM_2.5 data range was set to be [1,70] μg/m³. Furthermore, the monitoring locations of EPA and PA are plotted on the map with the population density, and housing units around the sampling locations in Figure 1. The total population, housing units, and median housing income were calculated in the census blocks, as defined by the U.S. Census Bureau [38]. The PA sensors are located in urban areas with higher populations and incomes.

The Pearson correlation coefficient was calculated to quantify the linear association between pairs of EPA, PA, and NOAA sensors using

where x_i and y_i are the individual sample points, and x̄ and ȳ are the means of the variables x and y, respectively. In the context of PM_2.5 measurements and meteorological data, x_i and y_i represent data from either PM_2.5 or meteorological variables.

To assess the significance of measurement errors between EPA and PA, we conducted Bland–Altman analysis [21]. Additionally, we performed Levene's test of equality of variances to determine whether the variability on the observed data from both EPA and PA sensors was statistically different [22].

Recognizing that the correction models using only RH and temperature do not uniformly account for the contribution of all sources to the PA data, particularly at the short-term component of PM_2.5 [14], we plan to investigate the source components impacting both short-term and baseline components of PM_2.5. To further examine this, we have separated the data into short-term and baseline components. The short-term component includes high-frequency data that is influenced by local anthropogenic sources such as traffic and short-term weather events. The baseline component, on the other hand, includes low-frequency data that are related to seasonal changes in weather, and changes in emission rates over time [23–25]. We ignored the medium-term component of both EPA and PA from the analysis as the medium-term component of raw PA data matches with raw EPA data.

The PM_2.5 time series data are separated into short-term and baseline components using the Kolmogorov-Zurbenko (KZ) filter technique [14, 23]. The KZ filter is a low-pass filter produced through repeated iterations of a moving average with parameters moving window (m), and iterations (p) also known as KZ_{m, p}. For details on the KZ filter formulation for air sensor data, please refer to [14].

The baseline part of PM_2.5(t) time series denoted as PM_{2.5, B}(t) and baseline part of meteorological variables' time series M_i(t) denoted as M_Bi are obtained by

The short-term part of PM_2.5(t) time series denoted as PM_{2.5, S}(t) and short-term part of meteorological time series M_i(t) denoted as M_Si(t) are obtained by

By separating the data into short-term and baseline components, we can analyze and examine how each component contributes to the overall variance of the time series data for both EPA and PA PM_2.5 [26]. The relative contributions of temporal components are obtained as follows:

where Var(i(t)) is variance of short-term, or baseline component, and Var(PM_2.5(t)) is variance of total PM_2.5 time series.

The short-term and baseline components of PM_2.5 can be combined with short-term and baseline components of meteorology to quantify the effect of meteorology and relatively estimate the effect of anthropogenic activities on PM_2.5 data [27–29]. The PM_2.5 data can be approximated as short-term and baseline PM_2.5 measurements as

The multiple linear regression (MLR) models for short-term and baseline components of PM_2.5 with short-term and baseline meteorology and anthropogenic activities can be written as

where,

M_Si(t) and M_Bi(t) are time series of the i^th meteorological variable for short-term and baseline components, respectively, and a₀, b₀, a_i, and b_i are regression model parameters to be estimated using a step-wise algorithm in MLR model. The residuals ϵ_S(t), ϵ_B(t) represent changes in PM_2.5 concentrations that cannot be attributed to meteorological variables present in the model and are mainly due to anthropological activities in the short-term and baseline components, respectively [27, 28]. To estimate the impact of meteorology and anthropogenic impact on both short-term and baseline PM_2.5(t), we built models considering PM_2.5(t) data as the response variable and meteorological data from nearby NOAA sensor as the predictor variable for each EPA sensor and PA sensor. We used the variance inflation factor (VIF) to assess the multicollinearity between the explanatory variables [30].

MLR models can only quantify the overall impact of meteorology on PM_2.5 measurements of both EPA and PA networks in short-term and baseline components. However, the question of which predictor most influences the data of both networks has no trivial answer due to the presence of correlated predictors. Correlation analysis is often used to examine the relationship between two variables. However, when there are many predictors, correlation analysis is not the best method to use. Here, we use the LMG measure proposed by Lindeman, Merenda, and Gold [31] and popularized by [32] to determine the relative importance of predictors.

The LMG measure uses sequential R², but it accounts for the dependence on orderings by averaging over all possible orderings. According to [33, 34], the variance decomposition for a linear model with k predictors can be defined as

and

where ν_j and ν_k are variances of each predictor M_i(t), and ρ_jk is the covariance of predictor j = 1, 2, 3, 4 with k = j+1, ..., 4.The R² for a model with predictors in set S is given as

The additional R² adding set X to a model with the predictors in set S is given by

where S and X are disjoint sets of predictors.

where r denotes permutations, r = 1, 2, ..., p!; seqR2(xk|r) denotes the sequential sum of squares for the predictors x_k in the ordering of the predictors in the r-th permutation.

The LMG measure for the k-th predictor x_k based on sequential sums of squares from all possible (p!) orderings for p predictors is given by

For example, for three explanatory variables (p=4), there are 24 different orderings (4!) and six different estimations (sequential sum of squares) for each explanatory variable. The relative importance of each explanatory variable is the mean of the six estimations. We applied the LMG measure, defined in Equation 16 on short-term and baseline components of PM_2.5 in Equation 9 and Equation 10 to get the relative contribution of each meteorological predictor on short-term and baseline PM_2.5, respectively.

For the entire workflow (reading and organizing data, descriptive analysis, and data analyses), we used the R software (R: A Language and Environment for Statistical Computing) (version 4.2.3), along with the following libraries in our coding: readxl, dplyr, tidyr, ggplot2, car, qqplotr, kza, stats, relaimpo. The “relaimpo” package was developed by [35], which can calculate the relative importance of predictor variables in multiple regression using the LMG measure and bootstrap confidence intervals.

This analysis utilizes PM_2.5 time series data from 5 EPA, 9 PA sensors, and 5 NOAA monitors located in Cook County, IL, from October 2019 to September 2021. The distribution of PM_2.5 data at each of the EPA and PA sensors over the entire analysis period is presented in Table 1. The overall distribution of PA sensors is broader compared to EPA monitors with higher mean PM_2.5 concentrations. Furthermore, to understand the overall linear relationship of EPA with PA, we applied the Pearson correlation Equation 1 to each pair of EPA and PA sensors. The correlation analysis was conducted for all possible combinations: PA with PA, EPA with EPA, and PA with EPA. The results of the correlation analysis are presented in Table 2. The PA sensor network shows correlations within the PA network with correlation coefficients ranging from 0.81 to 0.90 and with EPA coefficients ranging from 0.59 to 0.72. The correlation coefficient within the EPA network ranges from 0.51 to 0.67. We calculated the correlation coefficient of meteorological variables, as shown in Supplementary Table 2. Relative humidity (RH), temperature, wind speed (WS), and wind direction (WD) are all correlated with each other. Specifically, T and RH are negatively correlated, and WS is also negatively correlated with RH and T.

To assess the significance of measurement errors between EPA and PA, we conducted Bland–Altman analysis as shown in Supplementary Figure 1. The Bland–Altman analysis on EPA E2 and PA P6, the closest and comparable sensors pair, shows that measurements of E2 tend to report values that are lower than those reported by P6 and have significant measurement error. This is also reported in many earlier studies that PA sensors overestimate the measurements [12, 14]. Additionally, we applied Levene's test of equality of variances to assess whether the variances of the data from the EPA and PA sensors were statistically different. The results showed a p < 0.001 and an F-test value of 225.43. This indicates that the measurement variances from the EPA (E2) and PA (P6) sensors differ significantly. The difference in variances suggests that the measurements are inconsistent between the two sensor types.

To investigate the source components influencing the short-term and baseline fluctuations of low-cost sensor PM_2.5 data and to compare it with PM_2.5 data from reference monitors of EPA, we use the KZ filtering approach to separate the short-term and baseline components of the PM_2.5 time series at each selected EPA monitor and PA sensor, as well as meteorological variables including RH, temperature, WS, and WD from a nearby NOAA station, following Equations 2–5. The summary of KZ filtered short-term as baseline components is presented in Tables 3, 4. For illustration of temporal variations, one combination of EPA and PA datasets (E2 and nearby PA sensor P6), the short-term and baseline components are shown in Figure 2. The total PM_2.5 time series in Figure 2A for EPA sensor E2 has a range from 0 to 30μg/m³, whereas raw data from PA sensor P6 have an almost double range from 0 to 60μg/m³, which is also observed in the standard deviation of PA sensor data compared to EPA. Similarly, after decomposing PM_2.5 time series into short-term components as shown in Figure 2B and Table 3, the standard deviation is double in PA data compared to EPA. In the baseline component of PM_2.5 in Figure 2C, Table 4, the standard deviations in both datasets are similar, suggesting that the data have more variation in the short-term for PA sensors compared to EPA.

Our analysis of the time series decomposition revealed differences in the short-term and baseline components in both networks. By measuring the variation in each temporal component, we quantified the proportion of variation in each component relative to the total variation of the EPA and PA PM_2.5 time series data using Equation 6. The results of the relative contributions of short-term and baseline components to total data are presented in Figure 3A, Supplementary Table S3A. In Figure 3A, it can be observed that both PA and EPA PM_2.5 data have a greater relative contribution in the short-term component to the total variance, which is approximately 60–86%. However, comparing the two networks, the PA sensors have a relatively higher contribution than EPA to short-term variations. The relative contribution of the short-term component to total variance is greater in PA data, with a narrow distribution, indicating that the short-term component variance is largely uniform in the PA network and may be due to the capture of a local source that is independent of the sensor's location, as observed in correlations within the PA network in Table 2. On the other hand, the EPA sensors exhibit a broader variation in their short-term component, indicating that they capture local sources based on their location. The relative contribution of baseline to total variance is approximately 8–37%, but when comparing the two networks, the relative contribution of EPA sensors is more than that of PA sensors in the baseline component. The higher relative contribution of EPA to the baseline could be due to the impacts of meteorology as meteorology contributes to baseline trends of air quality data.

To understand and quantify the effect of meteorological conditions on the PM_2.5 data from PA and to compare with the data from EPA, we used MLR models with the stepwise forward selection algorithm. We included short-term and baseline temperature, RH, WS, and WD as meteorological predictors and short-term and baseline PM_2.5 data as response variables in our analysis. In the final model, only variables that were significant according to the stepwise forward selection algorithm were included. This technique involves adding variables one at a time based on their p-value and determines the optimal set of parameters for the model. The model performance was compared using the R² values. From the previous section, we noted that the relative contributions of short-term components were higher in PA sensors as shown in Figure 3A but looking at MLR models, meteorology has a greater impact on PA sensors' short-term variations as observed in Figure 3B. This implies that higher variations were indeed due to weather in the short-term component of low-cost PA sensors.

In Figure 3B, it can be observed that the short-term component of the PA sensors has R² ranging from 0.33 to 0.42. On average, this is 11% more R² than EPA sensors, as seen in Supplementary Table S3A. Likewise, the baseline component of PA sensors in Figure 3B, Supplementary S3A has higher R² ranging from 0.23 to 0.67. Despite having lesser relative contributions in the baseline component, it still has an average of 18% more R² than EPA sensors. This shows that weather is a higher contributor to the variance of PM_2.5 in both short-term and baseline components in low-cost sensors compared to reference monitors. The short-term and baseline components of all PA sensor PM_2.5 data have higher R² with meteorological parameters, indicating that weather influences PA sensors but is less responsive to anthropogenic emissions from traffic and other sources compared to EPA reference monitors.

We also used the variance inflation factor (VIF) to assess multicollinearity among the explanatory variables, as a complement to the LMG method. The VIF analysis indicates that when the model is fitted to the PM_2.5 time series without separating it into short-term and baseline components, all weather variables, including temperature, relative humidity (RH), wind speed (WS), and wind direction (WD), have VIF values below 2. However, after applying the KZ filter, the baseline component models show VIF values of T greater than 4 and WS greater than 2, while RH and WD remain below 2. In contrast, the short-term component has VIF values below 2 for all weather variables, suggesting that multicollinearity is present in the baseline component.

We utilized multiple linear regression models to determine that low-cost sensors are more sensitive to weather parameters compared to reference monitors, in both the short-term and baseline components. However, it is not possible to determine the individual influence of each meteorological predictor using MLR analysis due to their correlation with each other. Therefore, we used the LMG measure to determine the relative importance of each predictor in both short-term and baseline PM_2.5. The output of LMG(x_k) is partial R² for the variable that adds up to 1 for all predictors x_k, for k = 1, 2, ...., n. For our study case, k = 1, 2, 3, 4, for RH, temperature, WS, and WD. The LMG(x_k) measure was calculated using Equation 16, applied on short-term (Equation 9) and baseline (Equation 10) PM_2.5, and the results of LMG(x_k) measure for PM_2.5 time series, short-term, and baseline components of PM_2.5 are summarized and presented in Tables 5–7, Figures 4A, 4B, respectively. Based on the LMG measure, we found that wind speed (WS) is the most influential factor for both PA and EPA time series before the data are broken down into short-term and baseline components. However, if we exclude WS, relative humidity (RH) becomes the most influential factor. Once the time series is broken down into components, WS emerges as the most influential factor in almost half of the PA sensors (P1, P4, P5, and P6) in the short-term components, while the temperature is the most important factor in the short-term component of all EPA sensors except E1. In the baseline, RH consistently remains an important factor for PA sensors, except for P6 and P8. There are varying responses of meteorological factors in the baseline and short-term of both datasets.