Neural Models for Imputation of Missing Ozone Data in Air-Quality Datasets

Research Article Neural Models for Imputation of Missing Ozone Data in Air-Quality Datasets Ángel Arroyo ,1Álvaro Herrero,1Verónica Tricio,2 Emilio Corchado,3and MichaBWofniak4 1Department of Civil Engineering, University of Burgos, Burgos, Spain 2Department of Physics, University of Burgos, Burgos, Spain 3Departamento de Inform´ atica y Autom´ atica, University of Salamanca, Salamanca, Spain 4Department of Systems and Computer Networks, Wrocław University of Science and Technology, Wrocław, Poland Correspondence should be addressed to ´ Angel Arroyo; aarroyo[email protected] Received 5 December 2017; Accepted 31 January 2018; Published 8 March 2018 Academic Editor: Eloy Irigoyen Copyright©2018´ Angel Arroyo et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Ozone is one of the pollutants with most negative effects on human health and in general on the biosphere. Many data-acquisition networks collect data about ozone values in both urban and background areas. Usually, these data are incomplete or corrupt and the imputation of the missing values is a priority in order to obtain complete datasets, solving the uncertainty and vagueness of existing problems to manage complexity. In the present paper, multiple-regression techniques and Artificial Neural Network models are applied to approximate the absent ozone values from five explanatory variables containing air-quality information. To compare the different imputation methods, real-life data from six data-acquisition stations from the region of Castilla y Le´ on (Spain) are gathered in different ways and then analyzed. The results obtained in the estimation of the missing values by applying these techniques and models are compared, analyzing the possible causes of the given response. 1. Introduction and Related Work The ozone (O3) is an odorless, colorless, and highly reactive gas composed of three oxygen atoms. It is formed both in the Earth’s upper atmosphere (stratospheric ozone) and at ground level (tropospheric ozone). It can be “good” or “bad” for people’s health and for the environment, depending on its concentration levels and location in the atmosphere [1]. Stratospheric O3is formed naturally through the interaction of solar UltraViolet (UV) radiation with molecular oxygen (O2). Ground-level or “bad” ozone is not emitted directly into the air. In the 1950s, hydrocarbons and nitrogen oxides (NO𝑥) were identified as the two key chemical precursors of photochemical smog and its concomitant high concentrations of O3and other photochemical oxidant [2]. The majority of ground-level O3is formed from the photochemical oxidation of Volatile Organic Compounds (VOCs) in the presence of NO and other NO𝑥. Significant sources of VOCs are chemical plants, gasoline pumps, oilbased paints, autobody shops, and print shops. NO𝑥result primarily from high temperature combustion, and its most significant sources are power plants, industrial furnaces and boilers, and motor vehicles [3]. 1.1. Importance of Ozone. The O3exposition can cause damage in different ways. In the stratosphere, reduced O3levels asaresultofO 3layer depletion mean less protection from the sun’s rays and more exposure to UltraViolet B (shortwave) rays (UVB) radiation at the Earth’s surface [4]. The effects on human health of the O3layer depletion have been much analyzed, increasing the amount of UVB that reaches the Earth’s surface. UVB causes nonmelanoma skin cancer and plays a major role in malignant melanoma development. In addition, UVB has been linked to the development of certain cataracts, negative effects in patients with asthma, and Hindawi Complexity Volume 2018, Article ID 7238015, 14 pages https://doi.org/10.1155/2018/7238015 2Complexity other chronic respiratory disease. With respect to groundlevel O3, and its effects on human health, breathing O3can triggeravarietyofhealthproblems.Peoplewithasthmaand other chronic respiratory disease are a large and growing segment of the population and are also known to be especially susceptible to the effects of O3exposure. On days with high levels of O3, people with asthma tend to experience increased respiratory symptoms [3]. The layer O3depletion has also negative effects on the process of the development of plants, effects on the marine ecosystems like a direct reduction in phytoplankton production, negative effects on materials like biopolymers, and so forth. Tropospheric O3does not provide theprotectivefunctionthatitfulfillsinthestratosphere,being high reactivity. Its strong oxidizing capacity, when its levels rise above the natural background, can cause adverse effects in materials (derived from its corrosive effects), on vegetation and ecosystems. The present work focuses on tropospheric O3,whichis ariskfortheairquality[3].GiventheincreaseinO 3levels in the troposphere, it is currently considered one of the most important atmospheric pollutants. 1.2. Ozone Level Monitoring. Around the world there are numerous data-acquisition networks for the measurement of O3 levels and other pollutants, which consist of many stations in different locations where different sensors measure corresponding magnitudes. These network stations acquire data at periodic intervals of time (periods between ten and fifteen minutes are the most frequent ones) but frequently appear missing or corrupted data. In Europe, data are considered as corrupted when not meeting the Council Decision 97/101/EC of January 27, 1997 [5], which establish a reciprocal exchange of information and data from networks and individual stations measuring ambient air pollution within the Member States. Some of these networks provide information about the validity of the data, indicating through codes if the data is correct, it has not been possible to acquire, or it is corrupt, but in other occasions this type of information is not provided while the data are still missing. Some reasons for such failures havebeenpinpointed[6],namely,adamagedcable,thelossof proper electrical grounding, half-melted frost or snow on the dome, communications failure, and so forth. Some of these causes are temporary and may disappear spontaneously, but other ones require the intervention of a maintenance task force, and therefore errors persist for different periods of time.Theabsenceofvaliddatamayalsobeduetoreasons such as the following: mishandling of samples, low signalto-noise ratio, measurement error, nonresponse, or deleted aberrant value [7]. This is a problem for the analysis of the information coming from the measurement networks, and the imputation of these missing data [8] is necessary. Any of the variables acquired in network stations may suffer from the problem of the absence of data. If many data variables are omitted or corrupted in the same record, the whole sample must be withdrawn, when some models are applied [9], for subsequent tasks such as control, classification, forecast. Alternatively, if data for the same pollutant are missing in several adjacent rows, removing that variable may also be an alternative solution. In conclusion, having a complete set of data is necessary to perform a reliable study and to apply some models that cannot deal with missing data. 1.3. Missing Values and Related Work. The standard classification of missing data phenomenon [10] includes different situations: (i) Missing Completely At Random (MCAR), when the probability of an instance (case) having a missing value for a variable does not depend on either the known values or the missing data. (ii) Missing At Random (MAR), when the probability of an instance having a missing value for a variable may depend on the known values but not on the value of the missing data itself. (iii) Not Missing At Random (NMAR), when the probability of an instance having a missing value for a variable could depend on the value of that variable. As previous authors have pointed out, the complexity varies between these patterns of missing data [11]. Usually, in thecaseofair-qualitydata,missingvaluesareassociatedwith MAR or MCAR. The circumstances that may interfere with the acquisition of the data are many and not easily predictable [12]. To solve the missing data problem, a wide variety of differentmethodshavebeenapplieduptonow[8,10,13]. These imputation methods (IMs) are usually classified as follows: (i) Single imputation (SI): the method fills in one value for each missing one [12]. (ii) Multiple imputation (MI): multiple simulated values are generated at the same time [14]. The univariate and multivariate imputation methods differ in which the approximation of the missing values of the variable under study are calculated from the rest of the values of the very same variable (univariate) or using values of the rest of the variables (multivariate) [12]. With the aim of reducing the complexity of other MI applied methods [11], the present paper focuses on single and multivariate imputation for the O3magnitude in air pollution datasets. To do so, multiple-regression (linear and nonlinear) techniques together with Artificial Neural Networks (ANN) are applied to real-life datasets obtained from public airquality networks. Up to now, different Artificial-Intelligence (AI) techniques have been applied for imputation of missing data. In [7] imputation methods based on six different techniques are compared: K-Nearest Neighbors (KNN), Fuzzy K-Means (FKM), Singular Value Decomposition, Bayesian Principal Component Analysis (bPCA) and Multiple Imputations by Chained Equations. These methods are applied to four datasets split into two groups of various sizes: small datasets (Iris and E. coli) and large datasets (breast cancers 1 and 2). bPCA and FKM appeared to be the most robust imputation methods in the tested conditions. In [15] the accuracy of different imputation methods is evaluated: MissForest (MF) and Multiple Imputation based Complexity 3 on Expectation-Maximization (MIEM), along with two other imputation methods: Sequential Hot-Deck and Multiple Imputation based on Logistic Regression (MILR). The models are applied over fourteen binary datasets, with a range of missing data rates between 5% and 50%. The results from 10fold Cross-Validation (CV) show that the performance of the imputation methods varies substantially between different classifiers and at different rates of missing values. Although many imputation methods have been proposed up to now, scant attention has been paid to validate ANN for such a task, taking advantage of their regression capability [16]. Among these previous studies, ANN have been applied for the estimation of lost values in [17], where the main goal is identifying Learning Disabilities (LD) in children at early stages. In [18], authors proposed a SI approach relying on a Multilayer Perceptron (MLP) whose training is conducted with different learning rules, and a MI approach based on the combination of MLP and KNN. 24 real and simulated datasets from the UCI repository, the Promise repository, and mldata.org were exposed to a perturbation experiment with random generation of monotone missing data pattern. In [19] six different types of ANN are proposed as IM: MLP and its variations (the Time-Lagged Feedforward Network (TLFN)), the Generalized Radial-Basis-Function (GRBF) network, the Recurrent Neural Network (RNN), and its variations (the Time Delay Recurrent Neural Network (TDRNN)). Additionally, the Counterpropagation FuzzyNeural Network (CFNN) along with different optimization methods is applied for infilling missing daily total precipitation and extreme temperature series from 15 weather stations. The standard MLP and TLFN appear to provide the most accurate reconstruction of missing precipitation and daily extreme temperatures records with results for the Rcorrelation coefficient between the observed and the reconstructed daily series close to 1. In [20] a novel nonparametric algorithm named Generalized regression neural network Ensemble for Multiple Imputation (GEMI) is proposed. Additionally, a SI version of this approach (GESI) is proposed. The algorithms were tested on 98 synthetic and real-world datasets. All simulation results showtheadvantagesofGEMIascomparedwithconventional algorithms. GEMI has heavy memory storage requirements but outperformed other SI algorithms. In [21] fifteen real and simulated datasets are exposed to a perturbation experiment, based on the random generation of missing values. Several architectures and learning algorithms for the MLP are tested and compared with three classic imputation procedures: mean/mode imputation, regression, and hot-deck [22]. In [23] a methodology based on Gaussian Mixture Model (GMM) and Extreme Learning Machine (ELM) is developed andtestedonsomedatasetsfromtheUCIMachineLearning Repository and the LIACC regression repository. GMM is used to model the data distribution which is adapted to handle missing values, while ELM enables devising a Multiple Imputation strategy for final estimation. The combination ofGMMandELMisshowntobesuperiorinalmostall tested cases over the method based on conditional mean imputation. In[24]aSIapproachrelyingonaMLPandaMIapproach based on the combination of MLP and K-NN is proposed. The models are applied to 18 real and simulated datasets like domains such as biology, medicine, chemistry, electronics, social surveys, census, and business. For datasets with only quantitative variables MIMLP model provided the best results, with IMLP being the best method for datasets with categorical variables. In [25] a two-stage hybrid model for filling the missing valuesusingfuzzyc-meansclusteringandMLPisproposed. It is applied to a Wine dataset with a 1% to 5% of generated missing values and the accuracy of the model is checked using theMeanAbsolutePercentageError(MAPE).TheMAPE obtained for stage 2 (MLP regression to the obtained dataset as a result of applying fuzzy c-means in stage 1) is 4.95% for 1% missing-value records and 8.36% for 5% missing-value records. In the case of air-quality data, few imputation methods have been proposed up to now. In [13], an important set of SI: Listwise, Unconditional mean, Modified Median, Principal Component-based, Expectation-Maximization (EM) (Regularized-EM), and MI methods are applied to three datasets with the most important pollutant variables (NO, NO2,NO 𝑥,CO,O 3, PM10, and PM2.5) and a percentage of missing data among the 3.85% and the 23.52% depending on the year. Missing data of the eight variables are imputed in order to assess the effectiveness of the methods applied. In general, MI tends to yield more scattered values than its counterparts, mainly when the variables have many voids and they correlate poorly to the other variables like CO with 43.5% of missing data in 2006 and they correlate poorly to the other variables. In [11] some methods for the imputation of missing airquality data are compared: in the context of SI (linear, spline, and nearest neighbor interpolations), MI (regression-based imputation, multivariate nearest neighbor, Self-Organizing Maps (SOM), and Multilayer Backpropagation (MLBP) nets) and hybrid methods of the aforementioned. The dataset uses the most common pollutants: NO𝑥,NO 2,O 3, PM10, SO2, and CO concentrations, all on a time-scale of one per hour (hourly averaged), together with four meteorological parameters. The performance of the proposed univariate missing data interpolation was limited, and in general they were able tofillonlyveryshortgapsofcontiguousmissingdata.The general performance of the applied imputation methods was fair good when considering the pollutants (NO𝑥,NO 2,O 3, PM10, SO2, and CO) which are the most important ones in terms of air-quality modelling, but not so good regarding meteorological variables. The results suggested that SOM andMLBParethemethodsofchoiceforair-qualitydata imputation and even better results can be achieved by using the MI. 1.4. Main Contributions. The main contributions of this work are as follows: (i) Deep study of the real-life human health protection task in Spanish region of Castilla y Le´ on. (ii) Multisensor of O3data analysis. 4Complexity (iii) Experimental evaluation of the proposed approach based on multiple-regression techniques together with ANN models. To the best of authors knowledge, this is the first approach of imputation methods of O3based on both MLP and RadialBasis-Function Networks. The rest of this paper is organized as follows. Section 2 presents the techniques and models applied. Section 3 details the real-life case study that is addressed in present work, while Section 4 describes the experiments and results. Finally, Section 5 sets out the main conclusions and future work. 2. Regression Techniques and ANN Models In order to fill missing or corrupted values of O3in high dimensional datasets with air-quality information, two regression techniques and two ANN models have been applied in present study. This set of techniques applied as imputation methods is described in this section. 2.1. Regression Techniques. Linear regression attempts to model the relationship between two variables by fitting a linear equation to observed data. One variable is considered to be an explanatory variable, and the other is considered to be a dependent variable [26]. The general purpose of multiple regressions [27] is to learn more about the relationship between several independent or predictor variables and a dependent or criterion variable. 2.1.1. Multiple Linear Regression. Multiple linear regression (MLR) attempts to model the relationship between two or more explanatory variables and a response variable by fitting a linear equation to observed data [28]. Every value of the independent variable (𝑥) is associated with a value of the dependent variable (𝑦). The population regression line for 𝑝 explanatory variables 𝑥1,𝑥2,...,𝑥𝑝(1) isdefinedtobe 𝑢𝑦=𝛽0+𝛽1𝑥1+𝛽2𝑥2+⋅⋅⋅+𝛽𝑝𝑥𝑝.(2) This line describes how the mean response 𝑢𝑦changes with the explanatory variables. The observed values for y vary about their means 𝑢𝑦andareassumedtohavethesame standard deviation 𝜎.Thefittedvalues𝑏0,𝑏1,...,𝑏𝑝estimate the parameters 𝛽0,𝛽1,...,𝛽𝑝of the population regression line. Since the observed values for yvary about their means u𝑦, the multiple-regression models include a term for this variation. The model is expressed as DATA = FIT + RESIDUAL, where the “FIT” term represents the expression 𝛽0+ 𝛽1𝑥1+𝛽2𝑥2+⋅⋅⋅+𝛽𝑝𝑥𝑝. The “RESIDUAL” term represents the deviations of the observed values 𝑦from their means 𝑢𝑦, which are normally distributed with mean 0 and variance 𝜎. The notation for the model deviations is 𝜀. Formally, the model for multiple linear regression, given nobservations, is [28] 𝑌𝑖=𝛽0+𝛽1𝑥𝑖1 +𝛽2𝑥𝑖2 +⋅⋅⋅+𝛽𝑝𝑥𝑖𝑝 +𝜀𝑖 for 𝑖=1,2,...,𝑛. (3) 2.1.2. Multiple Nonlinear Regression. AMultipleNonlinear Regression (MN-LR) is a form of regression analysis in which observational data are modelled by a function which is a nonlinear combination of the model parameters and depends on one or more independent variables [29]. The data are fitted by a method of successive approximations. The parameters can take the form of an exponential, trigonometric, power, or any other nonlinear function. To determine the nonlinear parameter estimates, an iterative algorithm is typically used. 𝑦=𝑓(𝑋,𝐵)+𝜀, (4) where 𝐵represents nonlinear parameter estimates to be computed, 𝑋is the dependent or criterion variables, and 𝜀represents the error terms. 2.2. Artificial Neural Networks. Artificial Neural Networks (ANN), also known as Artificial Neural Systems (ANS), connectionist systems, adaptive networks, and distributed and parallel processing are simplified models of natural neural systems. The following definition, given by HechtNielsenin1989[30],formalizestheconceptofANN: An ANN is a parallel processing computer system distributed, consisting of a set of elementary processing units equipped with a small local memory and interconnected in a network through connections with associated weights. Each processing unit has one or more input connections and a single output connection that links to many collateral connections as desired. All processing associated with an elementary unit is a local, i.e. depends only on the values that take input signals from the unit and the internal state of the same. 2.2.1. Multilayer Perceptron (MLP). TheMLPconsistsofa system of simple interconnected neurons or nodes. The nodes are connected by weights and output signals which are a function of the sum of the inputs to the node modified by a simple nonlinear transfer, or activation, function. The architecture consists of several layers of neurons; the input layerservestopasstheinputvectortothenetwork.Theterms “input vectors” and “output vectors” refer to the inputs and outputsoftheMLPandcanberepresentedassinglevectors [31]. A MLP may have one or more hidden layers and finally anoutputlayer.MLParefullyconnected,witheachnode connected to every node in the next and previous layer. To perform a comprehensive comparison, the MLP is trained with the following algorithms: (1) Levenberg-Marquardt backpropagation (LM) Complexity 5 (2) Gradient Descent with momentum and adaptive learning rate backpropagation (GDX) [32] (3) Batch Training with weight and bias learning rules (TB) (4) Scaled Conjugate Gradient backpropagation (SCG) (5) Bayesian Regularization backpropagation (BR). 2.2.2. Radial-Basis-Function Networks (RBFN). In a RBFN [33] each unit in the hidden layer of this network has its own centroid, and, for each input vector 𝑥=(𝑥𝑙,𝑥2,...,𝑥𝑛), it computes the distance between 𝑥and its centroid. Its outputoftheunitiscalculatedasanonlinearfunctionofthis distance. Assuming that there are rinput nodes and moutput nodes, the overall response function without considering nonlinearity in an output node has the following form [34]: 𝑀 ∑ 𝑖=1𝑊𝑖∗𝐾(𝑥−𝑧𝑖 𝜎𝑖)=𝑀 ∑ 𝑖=1𝑊𝑖∗𝑔(󵄨󵄨󵄨󵄨󵄨󵄨󵄨󵄨𝑥−𝑧𝑖󵄨󵄨󵄨󵄨󵄨󵄨󵄨󵄨 𝜎𝑖), (5) where 𝑀∈Nis the number of units in the hidden layer, 𝑊𝑖∈R𝑚is the vector of weights linking the 𝑖th hidden-layer unit to the output nodes, xis an input vector, Kis a radially symmetric kernel function of a unit in the hidden layer, z𝑖 and 𝜎𝑖are the centroid and smoothing factor of the 𝑖th kernel node, respectively, and 𝑔:[0,∞)→Ris a function called the activation function, which characterizes the kernel shape. 3. Case Study In present study, data from air-quality stations in Castilla y Le´ on (CyL) are analyzed. CyL is a Spanish region located at the north-center of the Iberian Peninsula. It is composed of nine provinces and it is the most extensive region of Spain withatotalsurfaceof94,226squarekilometersandthesixth with more population: 2,435,797 habitants. Gross Domestic Product (GDP) in CyL represents the 5.3% of country’s GDP [35]. Climate in CyL approaches what is known as the continental ocean, characterized by cold winters and hot summers with short spring and autumn periods. CyL region provides a wide network of stations [36] for the acquisition of air-quality data. These data are public available according to the Open Data Initiative from the Spanish Government [37]. Stations from this network have some interesting characteristics: (1) Stations are classified in types: urban, background, andorientedtothevegetationprotection[36]. (2) These stations collect the fundamental air-quality pollutants, and among them is the O3,whichisthe objective pollutant of this study. Daily averages data [38] of each pollutant are provided in each location. (3) This data presents empty or corrupted data in all of its variablesinsomerowsandinareasonablepercentage to be estimated. Figure 1: Location of the six selected stations in CyL, by Google Maps. Inthepresentstudy,pollutantdatarecordedinsix different stations from the CyL network are analyzed. Daily data averages from years 2000 to 2008 have been selected. For some periods of time within the selected time window, data are not available for all the variables and, thus, the whole example is rejected for the study. Three of the stations are located in the center of the cities and labeled as urban stations; these stations are oriented to the protection of the human health. The other three stations are background stations and are also oriented to the protection of the human health. These stations measure a greater number of pollutants than the other type of stations and are the most important ones in terms of air quality, and many of themarenotcollectedatthestationsforthevegetation protection. This fact is important for the determination of the O3missingvalues,asthisgasisespeciallyharmfulforhuman health. The three urban stations considered in present study are as follows: (1) ´ Avila. “Bus Station” station. Geographical coordinates: 40.65914, −4.68237; 1150 meters above sea level (masl). (2) Aranda de Duero. “Jardines de Don Diego” station. Geographical coordinates: 41.67111, −3.68388; 801 masl. (3) Le´ on.“Avda.SanIgnaciodeLoyola”station.Geographical coordinates: 42.60388, −5.58722; 838 masl. The three background stations are as follows: (1) Burgos. “Fuentes Blancas” station. Geographical coordinates: 42.33611, −3.63611; 929 masl. (2) Segovia. “Acueducto” station. Geographical coordinates: 40.95555, −4.11055; 951 masl. (3) Medina del Campo (Valladolid). “Bus Station” station. Geographical coordinates: 41.31638, −4.90916; 721 masl. Figure1showsthelocationofthesixselectedstationsthat have been studied in the present paper. The pollutants gathered in the above-mentioned stations andanalyzedinthepresentstudyareasfollows: 6Complexity Table 1: Correlation matrix of the six variables in the dataset. O3CO NO NO2PM10 SO2 O31.000 −0.123 −0.161 −0.202 0.072 −0.013 CO −0.123 1.000 0.360 0.412 0.358 0.299 NO −0.161 0.360 1.000 0.540 0.233 0.330 NO2−0.202 0.412 0.540 1.000 0.330 0.257 PM10 0.072 0.358 0.233 0.330 1.000 0.251 SO2−0.013 0.299 0.330 0.257 0.251 1.000 Table 2: Percentage of missing and corrupted data for each one of the analyzed variables. O3NO NO2CO PM10 SO2 Missing 8.104% 8.020% 8.034% 8.554% 9.131% 8.196% Corrupted 1.857% 2.047% 1.815% 2.926% 2.413% 1.801% Total 9.961% 10.067% 9.849% 11.480% 11.544% 9.997% (1) Ozone (O3), 𝜇g/m3, secondary pollutant. See Section 1. (2) Carbon monoxide (CO), mg/m3,primarypollutant. It is an odorless, colorless gas formed by the incomplete combustion of fuels. When people are exposed to CO gas, the CO molecules will displace the oxygen in their bodies and lead to poisoning [39]. (3) Nitric oxide (NO), 𝜇g/m3,primarypollutant.NOis a colorless gas which reacts with ozone undergoing rapid oxidation to NO2, predominant in the atmosphere [39]. (4) Nitrogen dioxide (NO2), 𝜇g/m3,primarypollutant. From the standpoint of health protection, nitrogen dioxide has set exposure limits for long and short duration [39]. (5) Particulate matter (PM10), 𝜇g/m3,primarypollutant. These particles remain stable in the air for long periods of time without falling to the ground and can be moved significant distances by the wind. It is defined by the ISO as follows: “particles which pass through a size-selective inlet with a 50% efficiency cut-off at 10 𝜇m aerodynamic diameter. PM10 corresponds to the ‘thoracic convention’ as defined in ISO 7708:1995, Clause 6” [40]. (6) Sulphur dioxide (SO2), 𝜇g/m3,primarypollutant.It is a gas. It smells like burnt matches. Its smell is also suffocating. SO2is produced by volcanoes and in various industrial processes. In the food industry, it is also used to protect wine from oxygen and bacteria [39]. Primary pollutants are injected into the atmosphere directly. Secondary pollutants are formed in the atmosphere through chemical and photochemical reactions from the primary pollutants [36]. All data from these six variables were normalized for the study. On the other hand, all of them are highly decorrelated. Table 1 shows the correlation matrix of the six pollutants of the case study. It is worth mentioning that O3is the most independent pollutant, as its correlation coefficients with the rest of the variables are close to zero. There are a total of 13,526 samples, as one sample per day (daily average) was collected for the twelve months of every year, between years 2000 and 2008, in the six stations analyzed in this study. Missing or corrupted data appear in all the variables in some rows, which are omitted for the study. Table 2 shows the percentage of missing or corrupted data presentedineachvariableinthewholedataset. Allthesampleswithatleastonemissingorcorrupted value were removed from the dataset. 4. Experiments, Results, and Discussion The main target of this paper is to fill missing O3values in air pollution datasets. To do so, several imputation methods are comprehensively compared as described below. 4.1. Experimental Settings. The imputation methods described in Section 2 are applied to different datasets, all of them with the six variables described in Section 3: (1) The Whole Dataset (WD), comprising the 13,526 samples: results for this datasets are shown in Section 4.2. (2) The Season Dataset (SD): samples in WD are split in four subsets according to the four seasons of the year: spring (3,453 samples), summer (3,349 samples), autumn (3,295 samples), and winter (3,429 samples). Results for this dataset are shown in Section 4.3. (3) The Type station Dataset (TD): samples in WD are split into two subsets according to the type of the station where the data come from; “urban” (6,763 samples) or “background” (6,763 samples). Results for this datasets are shown in Section 4.4. For the three datasets, both statistical and neural imputation methods were applied and the performance is calculated through n-fold Cross-Validation (CV). The main idea behind CV is to split data, normally many times, for estimating the Complexity 7 Table 3: Linear regression and nonlinear regression results for the WD. Method MSE Time (s) Mean STD Mean STD MLR 5.490𝐸−06 2.311𝐸−08 0.089 0.216 MN-LR 5.415E −06 2.437𝐸−08 2.143 0.254 Table 4: Radial-basis function network results for the WD. #ofneurons MSE Time (s) Mean STD Mean STD 10 5.104E −06 2.723𝐸−08 0.050 1.091 30 5.108𝐸−06 1.273𝐸−08 0.050 1.091 50 5.105𝐸−06 2.513𝐸−08 0.047 0.098 risk, error, or performance of each algorithm. Part of data (the training samples) is used for training each algorithm, and the remaining part (the validation samples) is used for validating the algorithm(s). Then, CV selects the algorithm with the smallest estimated risk [41]. CV prevents from overfitting because the training sample is independent of the validation sample. The number of the 𝑘parameters (data partitions) was 10 for all the experiments in the present study. It means that 90% of the data are used for training and 10% for validation. In the case of neural models, the training process is repeated ten times (one for each fold). In the case of MLP, training is also repeated for each training algorithm (see Section 2.2). For all the experiments the Mean and the Standard Deviation (STD) of the Mean Square Error (MSE) for the ten folds arepresentedinTables3–11.TheMeanandtheSTDofthe execution time (in seconds) are also presented in Tables 3–11 for the 10 folds. For MLP and RBFN different network topologies have been applied: combinations of 10, 20, and 30 neurons in the hidden layer. Additionally, in the case of MLP, the model is trained 10 times with the same combination of parameters to reduce the effect of randomness and get more statistically significant results. 4.2. Results from the Whole Dataset. In this section, results in terms of MSE and execution time when applying MLR, MNLR, RBFN, and MLP to the WD are presented. In Tables 3 and 4, it can be observed that the MSE Mean values for the determination of the O3are very similar for the three applied methods (MLR, MN-LR, and RBFN). In the case of RBFN, slightly lower values of MSE are obtained, with the lowest one being obtained with 10 neurons in the hidden layer. Regarding execution times, the MN-LR method turns out to be the slowest and RBFN the quicker. The high values ofSTDfortheruntimeinthecaseofRBFNareduetothefact that it greatly varies from one fold to the others. As it can be seen in Table 5, the LM, SCG, and BR training algorithms present the lowest values of MSE Mean in all cases (10, 30, and 50 neurons) and very close to those shown in Tables 3 and 4. The lowest value of MSE was obtained with theLMlearningalgorithmand50neurons.Thelearning algorithm that attained the worst results (in terms of MSE) is GDX. With respect to execution time, the SCG algorithm attained the best results, while LM and BR are the second best ones, while TB was the slowest of the five algorithms. Obviously, the training algorithms take more time when 50 neurons are defined in the hidden layer, the TB algorithm being the one with greatest effect. 4.3. Results from the Season Dataset. In Tables 6–8 results of applying MLR, MN-LR, RBFN, and MLP to subsets with data from the four seasons of the year (spring, summer, autumn, and winter) are presented. In Tables 6 and 7 the 3 methods present similar values in MSE Mean, and the lowest MSE Mean is achieved by the RBFN with 50 neurons in the hidden layer for the summer season. The MSE Mean values are higher than that observed fortheWD.Theseasonoftheyearwiththelowestvalues of MSE Mean is the summer. One reason may be that there are few variations in pollution conditions during summer time. This is due to the small variation in weather conditions during summer as well as low industrial activity and traffic in urban areas due to vacation time. Furthermore, correlation coefficients in more than 20 pollutants analyzed in [42] are higher for measurements in the summer compared with correlations for measurements over all days combined. The season of the year with the worst results in the calculation of the MSE has been the autumn in the case of the two regression techniques and RBFN, although the differences between the three seasons (spring, summer, and autumn) is not significant. In terms of execution time, it is probed once again that MN-LR is the slowest method, while RBFN is the quickest one, returning very similar results for the four seasons of the year. In Table 8, similarly to Table 5, the training algorithms that achieve the best results in terms of MSE Mean are LM, SCG, and BR. LM achieved the best value of MSE Mean in 10 of the 12 cases shown in Table 8, being exceeded by BR by a minimum value for the winter and spring seasons with a configuration of 10 neurons. GDX records the worst MSE valuesinthe12casesshowninTable8.Again,thebestMSE Mean is obtained for the summer season, reducing the MSE Mean in comparison with those registered by RBFN. The season of the year with the worst results in the calculation of 8Complexity Table 5: Multilayer perceptron results for the WD. # of neurons Training algorithm MSE Time (s) Mean STD Mean STD 10 LM 4.731𝐸−06 5.143𝐸−08 0.070 0.287 GDX 1.129𝐸−04 6.825𝐸−05 0.317 0.287 TB 5.889𝐸−05 4.973𝐸−05 0.642 0.022 SCG 5.216𝐸−06 1.514𝐸−07 0.060 0.001 BR 4.775𝐸−06 1.092𝐸−07 0.074 0.003 30 LM 4.599𝐸−06 1.015𝐸−07 0.102 0.442 GDX 4.045𝐸−04 3.523𝐸−04 0.481 0.442 TB 4.223𝐸−05 2.087𝐸−05 1.420 0.025 SCG 5.162𝐸−06 1.19𝐸−07 0.063 0.001 BR 4.727𝐸−06 5.667𝐸−08 0.102 0.005 50 LM 4.512E −06 8.952𝐸−08 0.160 1.080 GDX 1.541𝐸−04 3.722𝐸−04 0.648 1.080 TB 4.812𝐸−05 3.032𝐸−05 2.156 0.051 SCG 5.014𝐸−06 8.322𝐸−08 0.068 0.001 BR 4.731𝐸−06 1.099𝐸−07 0.161 0.010 Table 6: Linear regression and nonlinear regression results for the Season Dataset. Subset Method MSE Time (s) Mean STD Mean STD Spring MLR 1.895𝐸−05 1.406𝐸−07 0.085 0.208 MN-LR 1.895𝐸−05 1.242𝐸−07 0.169 0.298 Summer MLR 2.101𝐸−05 1.447𝐸−07 0.085 0.215 MN-LR 1.343E −05 1.365𝐸−07 0.665 0.321 Autumn MLR 2.106𝐸−05 2.079𝐸−07 0.085 0.208 MN-LR 2.101𝐸−05 1.447𝐸−07 0.677 0.259 Winter MLR 1.895𝐸−05 1.406𝐸−07 0.088 0.214 MN-LR 1.895𝐸−05 1.242𝐸−07 0.168 0.274 Table 7: Radial-basis function network results for the Season Dataset. Subset # of neurons MSE Time (s) Mean STD Mean STD Spring 10 1.845𝐸−05 1.687𝐸−07 0.046 0.096 30 1.847𝐸−05 1.06𝐸−07 0.050 0.098 50 1.846𝐸−05 1.297𝐸−07 0.047 0.096 Summer 10 1.308𝐸−05 1.549𝐸−07 0.045 0.098 30 1.310𝐸−05 1.389𝐸−07 0.045 0.096 50 1.306E −05 1.344𝐸−07 0.047 0.096 Autumn 10 1.986𝐸−05 1.176𝐸−07 0.046 0.096 30 1.987𝐸−05 2.323𝐸−07 0.046 0.097 50 1.987𝐸−05 2.199𝐸−07 0.045 0.095 Winter 10 1.845𝐸−05 1.687𝐸−07 0.046 0.100 30 1.847𝐸−05 1.060𝐸−07 0.045 0.093 50 1.846𝐸−05 1.297𝐸−07 0.046 0.096 Complexity 9 Table 8: Multilayer perceptron results for the Season Dataset. Subset # of neurons Training algorithm MSE Time (s) Mean STD Mean STD Spring 10 LM 1.696𝐸−05 2.550𝐸−07 0.063 0.086 GDX 3.298𝐸−04 8.39𝐸−04 0.183 0.086 TB 5.870𝐸−05 2.590𝐸−05 0.383 0.049 SCG 1.959𝐸−05 3.439𝐸−07 0.056 0.001 BR 1.695𝐸−05 1.988𝐸−07 0.076 0.015 30 LM 1.531𝐸−05 4.107𝐸−07 0.068 0.217 GDX 6.652𝐸−04 1.051𝐸−03 0.210 0.217 TB 1.069𝐸−04 3.909𝐸−05 0.473 0.019 SCG 1.870𝐸−05 2.857𝐸−07 0.055 0.005 BR 1.562𝐸−05 3.915𝐸−07 0.071 0.225 50 LM 1.473𝐸−05 4.138𝐸−07 0.080 0.181 GDX 6.286𝐸−04 1.200𝐸−03 0.253 0.181 TB 1.564𝐸−04 1.402𝐸−04 0.770 0.075 SCG 1.809𝐸−05 3.768𝐸−07 0.056 0.001 BR 1.580𝐸−05 3.697𝐸−07 0.089 0.490 Summer 10 LM 1.009𝐸−05 1.203𝐸−07 0.059 0.0663 GDX 4.718𝐸−04 4.791𝐸−04 0.166 0.0663 TB 7.65𝐸−05 6.231𝐸−05 0.355 0.0396 SCG 1.217𝐸−05 2.749𝐸−07 0.053 0.0007 BR 1.010𝐸−05 1.436𝐸−07 0.064 0.0059 30 LM 9.171𝐸−06 2.713𝐸−07 0.065 0.440 GDX 8.070𝐸−04 1.446𝐸−03 0.202 0.440 TB 9.117𝐸−05 5.051𝐸−05 0.500 0.047 SCG 1.118𝐸−05 4.118𝐸−07 0.056 0.001 BR 9.822𝐸−06 3.107𝐸−07 0.073 0.004 50 LM 8.673E −06 3.284𝐸−07 0.083 0.225 GDX 4.572𝐸−04 9.864𝐸−04 0.253 0.225 TB 1.269𝐸−04 6.38𝐸−05 0.743 0.019 SCG 1.089𝐸−05 1.561𝐸−07 0.057 0.001 BR 9.851𝐸−06 2.165𝐸−07 0.085 0.003 Autumn 10 LM 1.622𝐸−05 2.146𝐸−07 0.061 0.101 GDX 1.598𝐸−04 3.096𝐸−04 0.168 0.101 TB 7.248𝐸−05 3.589𝐸−05 0.351 0.058 SCG 1.904𝐸−05 2.617𝐸−07 0.055 0.001 BR 1.628𝐸−05 7.680𝐸−07 0.071 0.009 30 LM 1.495𝐸−05 3.564𝐸−07 0.067 0.196 GDX 1.045𝐸−03 1.520𝐸−03 0.204 0.196 TB 1.09𝐸−04 6.737𝐸−05 0.506 0.048 SCG 1.808𝐸−05 2.756𝐸−07 0.054 0.001 BR 1.522𝐸−05 3.456𝐸−07 0.069 0.001 50 LM 1.401𝐸−05 1.926𝐸−06 0.079 0.307 GDX 5.676𝐸−04 1.700𝐸−03 0.240 0.307 TB 1.005𝐸−04 6.447𝐸−05 0.734 0.029 SCG 1.758𝐸−05 5.509𝐸−07 0.054 0.001 BR 1.559𝐸−05 6.591𝐸−07 0.083 0.103