ST-MVL: Filling Missing Values in Geo-Sensory Time Series Data

1 Introduction

Many sensors have been deployed in the physical world to continuously and cooperatively monitor the environment, such as for air quality and meteorology. These sensors generate massive geo-tagged time series data, helping humans to better understand surrounding conditions [Zheng et al., 2014]. However, we usually lose readings of a sensor (or a set of sensors) at unexpected moments, because of geo-sensor errors (e.g. power outages) or communication errors. Missing readings will not only affect real-time monitoring especially for emergency conditions, but also compromise the performance of further data analysis like prediction and inference.

*The paper was done when the first author was an intern in Microsoft Research under the supervision of the second author.

+Yu Zheng is the correspondence author of this paper.

Filling missing readings in a collection of geo-sensory time series data, however, is challenging for two reasons:

1) Readings can be absent at arbitrary sensors and timestamps. In some extreme cases, we may lose readings from a sensor at consecutive timestamps, e.g. 𝑠2 shown in Figure 1 A), or lose readings of all sensors in one (or more) timestamp(s) simultaneously, like 𝑡2 illustrates in Figure 1 A). We call these extreme cases block missing. It is very difficult for existing models to handle the block missing problem, as we may not be able to find stable inputs for a model. For instance, non-negative matrix factorization (NMF) cannot handle cases where the data of a column or a row are completely missing in a matrix.

2) Affected by multiple complex factors, sensor readings change over location and time significantly and non-linearly. First, the readings of sensors with a shorter distance may not always be more similar than those with a farther distance. As illustrated in Figure 1 B), 𝑠1 is closer to 𝑠2 than 𝑠3 in terms of geographical Euclidian Distance. As shown in Figure 1 C), however, the air quality readings of 𝑠2 are more similar to 𝑠3

than 𝑠1. The reason is that 𝑠2 and 𝑠3 are located in two regions with a similar geographical context, such as Point of Interests and traffic patterns, while 𝑠1 is deployed in a forest and with a lake between it and 𝑠2. These cases violate the First Law of Geography, devaluing some interpolation-based models, like Inverse Distance Weighting (IDW) [George and Wong, 2008]. Second, sensor readings fluctuate tremendously over time, sometimes coming with a sudden change. As illustrated in Figure 1 C), the reading of 𝑠2 at the 31st timestamp drops over 200 in two hours. Such sudden change is actually very important to real-time monitoring and further data analysis, but cannot be well handled by existing smoothing or interpolation methods.

To tackle these challenges, we propose a spatio-temporal multi-view-based learning (ST-MVL) method to collectively fill missing values in a set of geo-sensory readings from a spatio-temporal and global-local perspectives simultaneously. Our contributions are four-fold:

• x ST-MVL simultaneously considers 1) the spatial correlation between different time series and 2) the temporal correlation between readings at different timestamps in the same series, to generate a more accurate estimate.
• x ST-MVL integrates the advantages of global views, i.e. empirical models derived from the data over a long period, and those of local views, i.e. data-driven algorithms that are concerned with recent readings, to achieve better accuracy.
• x ST-MVL can handle the block missing problem, combining the four views in a multi-view learning framework.
• x We evaluate our method using Beijing air quality and meteorological data. The results demonstrate the advantages of our method compared with ten baselines.

2 Overview

Figure 2 presents 𝑚 sensors' readings at 𝑛 consecutive timesamps, which are stored in a form of matrix, where a row stands for a sensor and a column denotes a timestamp. An entry 𝑣𝑖𝑗 refers to the reading of 𝑖 𝑡ℎ sensor at 𝑗𝑡ℎ timestamp. In this matrix, 𝑣2𝑗 and 𝑣1,𝑗+1 are missing. We can estimate the reading of 𝑣2𝑗 based on its spatial neighborhoods, such as 𝑠1 and 𝑠3; we call it a spatial view. 𝑣2𝑗 can also be estimated based on readings of adjacent timestamps, like 𝑡𝑗-1 and 𝑡𝑗+1, etc.; we call it a temporal view. We can also use different time length of data for an estimation, enabling local and global views. For example, we consider adjacent readings of 𝑣2𝑗 from 𝑡𝑗-2 to 𝑡𝑗+2 in a local data matrix, which is regarded as a local view. We define the number of columns in a local data matrix as a window size 𝜔. Alternatively, we can take into account readings over a very long time period, e.g. from 𝑡1 to 𝑡𝑛, which is regarded as a global view. Local view captures instantaneous changes, whereas global view represents longterm patterns. If the values of an entire row or a column in a local data matrix are missing, we call it a block missing problem.

Global

Figure 2: Data matrix of sensor readings

As shown in Figure 3, ST-MVL consists of four views: IDW, Simple Exponential Smoothing (SES), User-based Collaborative filtering (UCF) and Item-based Collaborative filtering (ICF). The four views are then aggregated to generate a final estimate for a missing reading.

IDW computes an estimate for a missing reading of a sensor based on the readings of the sensor's spatial neighborhoods. SES estimates the missing reading of a sensor based on readings at other timestamps of the same sensor. As IDW and SES are actually empirical models derived from data over a long period of history, they denote a global spatial view and global temporal view on the missing reading respectively.

On the contrary, UCF estimates a missing reading only based on the local similarity between a sensor's recent readings and those of the sensor's spatial neighbors, where a sensor is regarded as a user. Likewise, ICF estimates a missing reading based on the local similarity between recent readings of different timestamps, where a timestamp denotes an item. As UCF and ICF only consider local similarity from a spatial and temporal perspective, they stand for local spatial and local temporal views respectively.

To leverage advantages of different views, we propose a multi-veiw learning algorithm that finds a linear combination of different views' predictions with minimal square errors.

3 The ST-MLV Method

In this section, we detail ST-MLV, using a running example to fill the reading of 𝑣2𝑗 shown in Figure 2. In the following, the term "missing value" equals "missing reading".

3.1 Global Spatial View - IDW

To model the global spatial view, we employ a statistic model, IDW, to interpolate a missing value based on its spatial neighborhoods. IDW assigns a weight to each available reading of geospatially adjacent sensors according to their distance to the target sensor, and aggregates these weights by Equation (1) to make prediction 𝑣̂𝑔𝑠.

where 𝑑𝑖 is the spatial distance between a candidate sensor 𝑠𝑖 and the target sensor, and α is a positive power parameter that controls the decay rate of a sensor's weight by 𝑑𝑖 -𝛼 . 𝑑𝑖 -𝛼 assigns a bigger weight to closer sensors' readings, for which a bigger 𝛼 denotes a faster decay of weight by distance.

Figure 4 presents the statistics that motivate us to use IDW to model global spatial view, using two datasets: air quality data and meteorological data in Beijing from May 2014 to May 2015. Here, we calculate the ratio between arbitrary two sensors' readings at the same timestamp. The ratio decreases as the distance between two sensors increases in both datasets. This actually follows the First Law of Geography [Tobler, 1970], i.e. Everything is related to everything else, but near things are more related than distant things, which is an empirical spatial correlation in geo-sensory data.

Here, we demonstrate IDW using the running example. Suppose two sensors 𝑠1 and 𝑠3 are near 𝑠2 with spatial distance 6.5km and 8.5km respectively. If setting 𝛼 = 1, weights to two sensors are 1/6.5 and 1/8.5. Then, we get the prediction 𝑣̂𝑔𝑠=130 after weighting average.

3.2 Global Temporal View - SES

To incorporate the global temporal view, we utilize SES to estimate the missing value based on the readings of same sensor at other timestamps. SES is frequently used in the time series domain as an exponential moving average model, formally defined as [Gardner, 2006]:

where 𝑡𝑗 is a time interval between a candidate reading 𝑣𝑗 and a target reading; 𝛽 is a smoothing parameter with a range of (0, 1). In general, 𝛽 ∗ (1 -𝛽)𝑡-1 gives a bigger weight to recent readings than distance ones, and a smaller 𝛽 denotes a slower decay of weight over the time interval. Traditional SES uses only predecessors of the target timestamp as input. Here, we extend it to use both predecessors and successors of a target timestamp. Given a target timestamp, our SES gives a weight 𝛽 ∗ (1 -𝛽)𝑡-1 to each reading of the same sensor, calculating 𝑣̂𝑔𝑡 by normalizing the weight according to Equation (3):

In implementation, we only select readings within a temporal threshold (12 hours), as distant reading is not very useful.

Our SES model is inspired by the observation from time series data. Figure 5 presents the ratio between arbitrary two readings at two different timestamps of the same sensor, using the same air quality data an meteorological data presented in Figure 4. Both curves in Figure 5 decrease as the time interval increases, showing an empirical temporal correlation in time series, i.e. readings of recent timestamps are more relevant than readings of distant timestamps.

Following the running example, suppose we find four timestamps (𝑡𝑗-2, 𝑡𝑗-1, 𝑡𝑗+1, 𝑡𝑗+2) for 𝑡𝑗. If setting β = 0.5, the weights for the four timestamps are (0.25, 0.5, 0.5, 0.25) respectively. Thus, the final result 𝑣̂𝑔𝑡 equals to 185.

3.3 Local Spatial View - UCF

UCF is a data-driven algorithm that has been wildly used in recommender systems. The general idea behind it is similar users make similar ratings for similar items [Su and Khoshgoftaar, 2009]. Following the running example shown in Figure 2, we regard a sensor as a user, constructing a local data matrix for 𝑣2𝑗 with adjacent readings [𝑣∗(𝑗-𝜔-1 2 ) …𝑣∗(𝑗+𝜔-1 2 ) ], where 𝜔 is the windows size. The similarity between two sensors (𝑠𝑖 , 𝑠2) is measured based on readings in the local data matrix, according to Equation (4):

where 𝑁𝑇 is the number of timestamps that two sensors both have readings. We then calculate a weighted average 𝑣̂𝑙𝑠 using the similarity as a weight, according to Equation 5.

UCF deals with special cases where the empirical spatial correlation mentioned in Section 3.1 does not hold, e.g. the case shown in Figure 1 B), capturing the time-dependent (i.e. local) spatial correlation between sensors' readings.

3.4 Local Temporal View - ICF

Regarding a timestamp as a item, ICF [Sarwar et al., 2001] calculates the similarity between two timestamps (𝑡1, 𝑡2) based on recent readings in the local data matrix built in the UCF model, according to Equation (6):

where 𝑁𝑆 is the number of sensors that have readings in both 𝑡1 and 𝑡2 . We use the similarity as a weight to calculate weighted average 𝑣̂𝑙𝑡 with the readings [ 𝑣𝑗1 …𝑣𝑗2 ] whose timestamp is within window size 𝜔 of the local data matrix, according to Equation (7).

ICF, as local temporal view, deals with special cases where the empirical temporal correlation mentioned in Section 3.2 does not hold, using the time-dependent (i.e. local) temporal correlation learned from recent data.

3.5 Multi-View Learning

ST-MVL integrates the predictions of the aforementioned four views to generate a final result, through a multi-view learning algorithm, according to Equation (8):

where 𝑏 is a residual and 𝑤𝑖 ( 𝑖 =1,2,3,4) is a weight assigned to each view. Algorithm 1 presents the procedure of ST-MVL. When a dataset encounters a block missing problem, in which ICF and UCF do not work very well, we leverage IDW and SES to generate an initial value for those missing entries (see Line 3). Then, we predict each missing entry using ICF, UCF, IDW and SES respectively (Line 4-9), combining the four predictions based on a linear-kernelbased multi-view learning framework (see Line 10 and Equation 8). We train the model for each sensor respectively, by minimizing the linear least square error [Lawson and Hanson, 1974] between predictions and ground truth.

Algorithm 1 ST-MVL

Input: Original Data Matrix 𝑀, 𝜔, 𝛼, 𝛽; Output: Final Data Matrix; 1. O ← Get_All_Missing_Values(𝑀); 2. If there are block missing problem 3. 𝑀 ← Initialization(𝑀, 𝛼, 𝛽); 4. Foreach target t in O 5. 𝑣̂𝑙𝑠 ← UCF(𝑀, 𝑡,𝜔); 6. 𝑣̂𝑙𝑡 ← ICF(𝑀, 𝑡, 𝜔); 7. 𝑣̂𝑔𝑠 ← IDW(𝑀, 𝑡, 𝛼); 8. 𝑣̂𝑔𝑡 ← SES(𝑀, 𝑡, 𝛽); 9. 𝑣̂𝑚𝑣𝑙 ← Muti-view_Learning(𝑣̂𝑙𝑠, 𝑣̂𝑙𝑡, 𝑣̂gs, 𝑣̂𝑔𝑡); 10. Add 𝑣̂𝑚𝑣𝑙 into 𝑀;

//using IDW or SES 11. Return 𝑀;

4 Experiments

4.1 Datasets and Ground Truth

We evaluate our model based on two real datasets: air quality and meteorological data in Beijing from 2014/05/01 to 2015/04/30 [Zheng et al., 2015], which has 8,759 timestamps respectively. The air quality data is collected at 36 air quality monitoring stations, each of which generates a reading every hour, as depicted in Figure 6 A). The meteorological data is collected by 16 sensors shown in Figure 6 B) every hour.

From the two datasets, we select 4 properties, consisting of PM2.5, NO2, Humidity and WindSpeed, to fill missing values whose statistics are shown in Table 1. There are two missing situations: block missing and general missing. The former is further comprised of spatial block missing and temporal block missing; the two may have some overlap. The spatial block missing is referred to as records with all sensors' readings simultaneously absent; the temporal block missing is records of the same sensor with data absent in a certain temporal window size (see Figure 2, e.g. there are 3.15% of missing values in PM2.5 property when ω = 11). General missing is the missing values except for the block missing. For example, about 13.25% of sensor readings in the air quality dataset lose PM2.5 property, including 8.15% general missing and 2.15% spatial block missing.

We partition the 1-year data into two parts, using the 3rd, 6th, 9th and 12th months as a test set and the rest for a training set. To train our model, we select entries (from the training set) whose local matrix is not totally absent of data at the row and column where the entry stays. As illustrated in Figure 7, for the test set, we accommodate each month's data into a matrix respectively. We find all missing entries (denoted by red X) in each month's data. Then, we drop out the values of entries (in the next month's matrix e.g. 𝑣32 in March) that correspond to its predecessor's (i.e. Feb.) missing entries. If those entries are not absent in the next month (e.g. 𝑣32=59 in March), their values are used as a ground truth to measure the accuracy of our predictions. Note that we cannot randomly drop out some non-missing entries in the test set to generate ground truth, as real-world missing patterns are not random.

4.2 Baselines & Metrics & Prameters

We compare our method with 10 baselines shown in Table 2. ARMA& SARIMA: ARMA predicts a missing value based on the readings of a few timestamps ago. SARIMA generalizes ARMA, considering the seasonal factor in time series. stKNN: This uses the average readings of its k nearest spatial and temporal neighbors as a prediction (𝑘=6).

Kriging: Kriging interpolates a missing value with available readings of geospatially adjacent sensors.

AKE (Applying k-nearest neighbor estimation) [Pan and Li, 2010]: It uses a regression model to combine the readings of a missing entry's neighbors.

DESM: Combines previous readings of a missing sensor and its neighbors' readings at current time linearly, using statistic based spatio-temporal correlations. [Gruenwald et al., 2010]. IDW+SES: Fill missing entries in a sensor's local matrix by IDW and SES respectively; return the average of predictions. CF: Initialize a local matrix by IDW+SES; apply UCF and ICF to generate a prediction respectively; the return is the average of the two predictions.

NMF: Initialize a local matrix by IDW+SES; fill the local matrix by a NMF.

NMF-MVL: Similar to our method, but it replaces two local views with a NMF.

Evaluation Metrics: We measure our mehtod by Mean Absolute Error (MAE) and Mean Relative Error (MRE),

where 𝑣𝑖 ̂ is a prediction and 𝑣𝑖 is the ground truth; 𝑛 is the number of cases.

Parameter Settings: We test different 𝛼 for IDW, 𝛽 for SES, and 𝜔 for UCF & ICF, finding a best setting for them, e.g. when 𝛼=4, 𝛽=0.85, and 𝜔=7 achieve the best performance in PM2.5 property.

4.3 Results

4.3.1 Results of Combination Methods

Table 3 presents the results of different combination methods of four views, based on PM2.5 property. ST-MVL outperforms all kinds of combinations, bringing a significant improvement beyond the best single view SES and the best combination of two views (IDW+SES). ST-MVL also achieves a slightly better performance in filling in block missing than others. In the meantime, the combination of two views, e.g. IDW+SES, is better than using a single view; same for UCF+ICF. The four views are complementary to each other, containing long-term patterns, knowledge derived from recent contexts, the spatial correlation between locations, and the temporal correlation between timestamps. Additionally, IDW and SES have better results than UCF and ICF respectively, showing the ability of empirical models dealing with general cases.

Note: Spa means Spatial; tem means Temporal.

4.3.2 Overall Results

Table 4 presents the comparison between our method and ten baselines, where ST-MVL outperforms all baselines under all circumstances. Sudden changes denote the readings changing over the last timestamp by a threshold (e.g. 50 in our settings). We use IDW+SES to initialize the missing values for datadriven algorithms, such as CF and NMF, when the latter faces block missing. Individually, NMF outperforms CF, as it considers user and item similarities simultaneously. However, when combined with other views in the multiview learning framework, CF brings a bigger improvement over NMF, i.e. ST-MVL outperforms NMF-MVL.

4.3.3 Impact of Parameters for IDW & SES

Figure 8 A) plots the performance of IDW changing over 𝛼 in PM2.5 and WindSpeed datasets. PM2.5 has a minimum MRE when 𝛼=4, and WindSpeed has a minimum MRE when 𝛼=1. Figure 8 B) shows the performance of SES changing over 𝛽. When 𝛽=0.85 and 0.55, PM2.5 and WindSpeed have a minimum MRE respectively. As PM2.5 fluctuates more strongly than wind speed in both spatial and temporal spaces, the latter should be given a smaller 𝛼 and 𝛽.

4.3.4 Impact of Window Size for CF

Figure 9 shows the impact of window size on CF. A large window size may lose time-dependency, but a very small window size may not capture the similarity between locations or time stamps. That is why MAE decreases first and then increases after 𝜔=7.

4.3.5 Results in Other Datasets

Table 5 presents the results of all properties of the two datasets. ST-MVL outperforms IDW+SES in all properties, achieving an improvement of 0.009 in MRE with the multi-view learning. It verifies the advantages of our method integrating empirical models and data-driven algorithms.

5 Related Work

5.1 Methods Considering a Single View

Spatial models: Inverse distance weighting, linear regression and Kriging are three wildly used methods. [Chen and Liu, 2012] estimates the rainfall distribution using inverse distance weighting, and [Wu and Li, 2013] fills the missing temperature using ordinary Kriging. These methods leverage the spatial correlation derived from the spatial coordinates data or learnt from data.

Temporal models: Methods in this category can be classified into non-feature-based and feature-based. The former, such as SES, ARMA and SARIMA, solely considers a sensor's readings [Ceylan et al., 2013]. The latter, such as Graph model and regression model [Lee et al., 2008; Fung, 2006], exploits a function of features, considering the temporal correlation of readings over time.

Beyond the aforementioned two categories of methods solely considering spatial information or temporal information, our method simultaneously utilizes spatial and temporal correlations, and therefore is more capable of filling in missing values.

5.2 Methods Considering Spatio-Temporal Views

Many statistic-based methods, such as mean imputation and last seen [Little and Rubin, 2014], can be used to impute the missing data. In addition, there are many data-driven methods for filling missing values. [Gruenwald et al., 2010] proposes a DEMS method to handle the missing data by time-dependent spatial and temporal relationships. [Pan and Li, 2010] proposes a nearest neighbors imputation method entitled AKE to describe the spatial correlation and incorporate the temporal information when meeting no enough spatial information. [Wang et at., 2006; Ma et al., 2007] unify user-based and item-based collaborative filtering to fill missing values in the recommendation system. [Johan at el., 2014] takes advantage of matrix factorization and a salient version called empirical orthogonal functions model is proposed by [Antti et al., 2010]. In general, statistic-based methods are derived from history data, characterizing the global information; Data-driven methods focus on the recent contexts, characterizing the local information. Our method combines the statistic models and data-driven models in a multi-view learning framework to generate a more accurate estimate.

6 Conclusion

In this paper, we propose a multi-view learning-based method to fill missing values for geo-sensory time series data, simultanesouly considering spatial, temporal, global and local views. We evaluate our method ST-MVL based on Beijing's air quality and meteorological data. ST-MVL has a mean relative error around 17% for PM2.5 and NO2, and about 8% for humidity, outperforming baseline methods that solely consider a single view by 26% on average. ST-MVL also surpasses those that only combine two out of the four views by 10% in MRE on average. The code and datasets have been released at: http://research.microsoft.com/apps/pubs/?id=264768.

Acknowledgments

This work was supported by the China National Basic Research Program (973 Program, No. 2015CB352400) and the National Natural Science Foundation of China (Nos. U1401258, 61573292).