Improving Accuracy of Anomaly Detection in Spatial-Temporal Population Data through SHAP Values of Reconstruction Errors⋆ Ryo Koyama1 , Tomohiro Mimura1 , Shin Ishiguro1 , Keisuke Kiritoshi2 , Takashi Suzuki1 and Akira Yamada1,* 1 NTT DOCOMO, INC, Sanno Park Tower, 2-11-1 Nagatacho, Chiyoda-ku, Tokyo, Japan 2 NTT Communications Corporation, Otemachi Place West Tower, 2-3-1 Otemachi, Chiyoda-ku, Tokyo, Japan Abstract When accidents, disasters, or other large-scale events occur, they significantly disrupt traffic, leading to congestion and reduced mobility. To effectively address these issues, it is crucial to precisely detect the underlying causes of these disruptions through the analysis of human mobility data. A common approach in anomaly detection is to employ dimensionality reduction techniques to compute reconstruction errors. However, the reconstruction errors generated by traditional methods are influenced by the correlations among features, which may obscure the true causes of anomalies. To overcome this limitation, we introduce an approach that calculates the SHAP (SHapley Additive exPlanations) values associated with the reconstruction errors resulting from dimensionality reduction. We conducted experiments using a dataset of human mobility patterns to evaluate the effectiveness of this method. The results demonstrate that our approach provides a more accurate explanation of anomalies compared to conventional methods. Keywords Anomaly Detection, Spatial-Temporal, SHAP, Dimensionality Reduction, Reconstruction Error 1. Introduction accuracy and interpretability of anomaly detection results, while also extending the applicability of Takeishi’s approach Anomaly detection in spatial-temporal population data has to high-dimensional data. gained significant attention in recent years due to its poten- The main contributions of this paper are as follows: tial applications in urban planning, disaster management, and public safety. Accurate identification of unusual pat- • We introduce a novel anomaly detection method terns in human mobility can help authorities respond more that integrates dimensionality reduction with SHAP effectively to disruptive events such as accidents and disas- values to identify anomalies in spatialtemporal pop- ters. However traditional anomaly detection methods often ulation data. face challenges in capturing the complex spatial and tempo- • We evaluate the effectiveness of our approach using ral dependencies in high dimensional population data. a real-world dataset of human mobility patterns in Previous studies have explored various approaches to a major urban area, demonstrating its superiority anomaly detection in spatial temporal data. Ochiai et al.[1] compared to traditional reconstruction error-based proposed a method that utilizes mesh-based population data methods. derived from mobile communication records to detect non- • We extend the applicability of Takeishi’s Shapley designated evacuation centers during disasters. Their ap- value-based approach to high-dimensional spatial- proach relies on significant reconstruction errors in anomaly temporal data, enhancing its utility for real-world scenarios, which are trained only with data representing scenarios. normal conditions. While this method demonstrates poten- tial, it may not effectively capture the underlying causes of The remainder of this paper is organized as follows. Sec- anomalies. tion 2 provides an overview of related work in anomaly On the other hand, Takeishi[2] demonstrated the effec- detection and spatial-temporal data analysis. Section 3 tiveness of using Shapley values to explain the causes of describes our proposed methodology in detail. Section 4 anomalies in dimensionality reduction scenarios. By apply- presents the experimental setup. Section 5 discusses the ing Shapley values to one-dimensional health data, such as experimental results. Finally, Section 6 concludes the paper myocardial infarction and breast cancer records, Takeishi’s and discusses future research directions. method provides a more interpretable understanding of anomaly detection results. However, the applicability of this approach to high-dimensional spatial-temporal data 2. Related Work has not been explored. Building upon the insights from Ochiai et al. and Takeishi, The detection of anomalies in urban population flows has we propose a novel anomaly detection framework that com- been extensively explored using diverse data sources, in- bines the strengths of both approaches. Our method inte- cluding road sensors[3], surveillance cameras[4], and social grates SHAP (SHapley Additive exPlanations) values with media data[5]. While road sensors and surveillance cam- dimensionality reduction to identify and explain anomalies eras prove effective for identifying local abnormalities, their in spatial-temporal population data. By leveraging the ex- broader application for city-wide anomaly detection is ham- planatory power of SHAP values, we aim to improve the pered by high installation and maintenance costs. On the other hand, social media data facilitates multimodal anomaly STRL’24: Third International Workshop on Spatio-Temporal Reasoning detection methods, such as the integration of bike-sharing and Learning, 5 August 2024, Jeju, South Korea ⋆ and taxi usage history[6], and the semantic interpretation Current affiliation: NTT Corporation of location-based anomalies identified through social media $ ryou.koyama.aw@nttdocomo.com (R. Koyama) analytics[7].  0009-0007-4157-1634 (R. Koyama) © 2024 Copyright for this paper by its authors. Use permitted under Creative Commons License Attribu- tion 4.0 International (CC BY 4.0). CEUR ceur-ws.org Workshop ISSN 1613-0073 Proceedings This research leverages population data extracted from variance, PCA ensures that the most important features of communication logs between mobile devices and base sta- the data are preserved, allowing for effective dimensionality tions to enhance anomaly detection capabilities. In contrast reduction and subsequent analysis. to road sensors and surveillance cameras, mobile device data captures a wide array of individual behaviors throughout the 3.2. Reconstruction Error entire city, thus offering a more comprehensive solution for detecting anomalies. For instance, Yabe et al.[8] utilized sta- Using the principal components obtained from PCA, we tistical methods to detect the emergence of non-designated can perform dimensionality reduction and reconstruction. shelters during disasters, although their studies lacked a Consider a test data vector y ∈ R𝑑 . The reduced repre- quantitative assessment of accuracy. Similarly, Ochiai et sentation y𝑟𝑒𝑑 ∈ R𝑝 is obtained using the mapping func- al. used mobile phone-based population data to detect non- tion 𝑟𝑒𝑑𝑢𝑐𝑒 : R𝑑 → R𝑝 , and the reconstructed vector designated evacuation sites during disasters by focusing on ŷ ∈ R𝑑 is computed using the reconstruction function reconstruction errors. However, as Takeishi has pointed 𝑟𝑒𝑐𝑜𝑛𝑠𝑡𝑟𝑢𝑐𝑡 : R𝑝 → R𝑑 as follows: out, these reconstruction errors, significantly influenced by ŷ = 𝑟𝑒𝑐𝑜𝑛𝑠𝑡𝑟𝑢𝑐𝑡(𝑟𝑒𝑑𝑢𝑐𝑒(y)) (1) feature interactions, may not accurately pinpoint anomaly locations. To overcome this limitation, Takeishi introduced The reconstruction error 𝑒y , a measure of fidelity of re- a method that employs Shapley Values to elucidate the ori- construction, is defined as the squared Euclidean distance gins of anomalies within dimensionality reduction models. between the original and reconstructed vectors: 𝑑 ∑︁ 3. Methodologies 𝑒y = ‖y − ŷ‖22 = (𝑦𝑖 − 𝑦ˆ𝑖 )2 (2) 𝑖=1 This section details our proposed methodology, SHAP Val- This error metric helps identify significant deviations from ues of Reconstruction Error, for anomaly detection, incorpo- normal patterns, which are potential indicators of anoma- rating SHAP values derived from reconstruction errors. We lies. begin by describing Principal Component Analysis (PCA) as Additionally, the reconstruction error for each feature 𝑖 the foundation for dimensionality reduction. Subsequently, is calculated as: we outline the traditional method based on reconstruction 𝑒𝑦𝑖 = |𝑦𝑖 − 𝑦ˆ𝑖 | (3) errors. Then, we introduce our enhanced approach that This feature-specific error is employed to identify which integrates SHAP values to improve the accuracy and ef- specific features are exhibiting anomalies. fectiveness of anomaly detection. Finally, we explain the calculation of SHAP values for multi-dimensional objec- tive variables, extending the methodology to more complex 3.3. SHAP Values of Reconstruction Error scenarios. Typically, the SHAP value 𝜑𝑖 (𝑥) for each feature 𝑖 is cal- culated using all features of the instance 𝑥 as explanatory 3.1. Principal Component Analysis variables and the prediction 𝑦ˆ as the objective variable as follows: PCA is a widely-used technique for reducing the dimension- ality of high-dimensional data while preserving the most 𝜑𝑖 (𝑥) = 𝑠ℎ𝑎𝑝_𝑣𝑎𝑙𝑢𝑒(𝑖; 𝑦ˆ, 𝑥) (4) significant features. By projecting the data onto a lower- dimensional space, PCA identifies the principal components This formulation allows for measuring the impact of fea- that capture the maximum variance in the data. ture 𝑖 on the prediction 𝑦ˆ, providing a concrete metric for Given a dataset X ∈ R𝑛×𝑑 with 𝑛 samples and 𝑑 features, understanding the significance of each feature in the model. the PCA process involves the following steps: Similarly, the SHAP value of reconstruction error 𝜓𝑖 (𝑥) is calculated using all features of the instance 𝑥 as explanatory 1. Standardize the Data: Subtract the mean of each variables and the reconstruction error 𝑒𝑦 as the objective feature from the dataset to center the data around variable using the following function: the origin. 2. Compute the Covariance Matrix: Calculate the 𝜓𝑖 (𝑥) = 𝑠ℎ𝑎𝑝_𝑣𝑎𝑙𝑢𝑒(𝑖; 𝑒𝑦 , 𝑥) (5) covariance matrix C = 𝑛1 X𝑇 X. This formulation measures the impact of feature 𝑖 on 3. Perform Eigenvalue Decomposition: Decom- the reconstruction error 𝑒𝑦 , providing a concrete metric pose the covariance matrix into eigenvalues and for evaluating the significance of each feature in anomaly eigenvectors: C = VΛV𝑇 , where Λ is a diagonal detection. matrix containing the eigenvalues, and V is a matrix whose columns are the corresponding eigenvectors. 4. Select Principal Components: Choose the top 3.4. SHAP Values of Reconstruction Error 𝑝 eigenvectors corresponding to the largest eigen- for Multidimensional Obejective values to form the principal components. These Variables components maximize the variance retained in the lower-dimensional space. The SHAP value of Reconstructtion Error 𝜓𝑖 (𝑥) for feature 𝑖, when the objective variable is represented in one dimension, In this study, we set the threshold for the variance to be is calculated using the following equation [9]: retained at 90%. This means that we select the number of principal components 𝑝 such that the cumulative variance 𝜓𝑖 (𝑥) = ∑︁ (𝑑 − |𝑆| − 1)!|𝑆|! [𝑓 (𝑥𝑆∪𝑖 ) − 𝑓 (𝑥𝑆 )] explained by these components is at least 90%. By retaining 𝑆⊆{1,...,𝑑}∖𝑖 𝑑! the principal components that explain the majority of the (6) In this equation, 𝑓 (𝑥𝑆∪𝑖 ) denotes the model’s predicted value when the feature set 𝑆 includes feature 𝑖, and 𝑓 (𝑥𝑆 ) 𝐾 ⎡ represents the predicted value when the set 𝑆 is used with- (𝑘) 1 ∑︁ ∑︁ (𝑑 − |𝑆| − 1)!|𝑆|! 𝜓𝑔 (𝑥) = ⎣ out feature 𝑖. |𝑆| denotes the number of elements in the 𝐾 𝑘=1 𝑆⊆{1,...,𝑑}∖𝑔 𝑑! feature subset 𝑆, and 𝑑 is the total number of features. ]︁ [𝑓 (𝑘) (𝑥𝑆∪𝑔 ) − 𝑓 (𝑘) (𝑥𝑆 )] (as defined in Equation 7) Subsequently, the SHAP Values of Reconstruction Error (𝑘) (9) 𝜓𝑖 (𝑥) for feature 𝑖 impacting the 𝑘-th dimension of the objective variable is calculated as the average difference Rankings for each grid are obtained by comparing these between the model predictions with and without feature scores against all others in the dataset. 𝑖, across all combinations of features, thus extending equa- tion 6 into the multidimensional context of SHAP Values of Hits@k Calculations For both methods, Hits@k is de- Reconstruction Error as follows: fined and calculated separately for each method to assess the efficacy in identifying the true anomalies within the top 𝑘 ranks of predicted anomalies. The total number of test in- ⎡ 𝐾 (𝑘) 1 ∑︁ ∑︁ (𝑑 − |𝑆| − 1)!|𝑆|! 𝜓 𝑖 (𝑥) = ⎣ 𝐾 𝑘=1 𝑆⊆{1,...,𝑑}∖𝑖 𝑑! stances, denoted as 𝑁 , is used to normalize the calculations, ensuring that the results are proportional to the size of the test dataset. The calculations are as follows: ]︁ [𝑓 (𝑘) (𝑥𝑆∪𝑖 ) − 𝑓 (𝑘) (𝑥𝑆 )] (7) 𝑁 This formulation allows for measuring the impact of fea- 1 ∑︁ Hits@𝑘reconst = 1(rank(𝑒𝑦𝑔𝑖 ) ≤ 𝑘) (10) ture 𝑖 across different combinations of features on each 𝑁 𝑖=1 dimension of the recostruction error. By doing so, it pro- 𝑁 vides a concrete metric for elucidating the significance of 1 ∑︁ Hits@𝑘SHAP = 1(rank(𝜓𝑔(𝑘) (𝑥)) ≤ 𝑘) (11) feature 𝑖 in anomaly detection, offering detailed insights 𝑁 𝑖=1 𝑖 into the causes of anomalies in specific dimensions. where 1(·) is the indicator function, and 𝑁 represents the total number of test instances. These metrics facilitate a di- 4. Preliminaries rect comparison of the methods’ effectiveness in accurately detecting anomalies. 4.1. Definition: Spatial-Temporal Population Data 5. Experiments This study utilizes Mobile Spatial Statistics (MSS) [10], representing population counts recorded across a two- dimensional geographic grid. Each record, denoted as (𝑔, 𝑡, 𝑝𝑜𝑝𝑔,𝑡 ), indicates the population count 𝑝𝑜𝑝𝑔,𝑡 at grid 𝑔 and timestamp 𝑡. 4.2. Problem Statement The goal of this study is to assess the accuracy of anomaly detection in spatial-temporal population data. We compare two methodologies: a traditional approach using recon- struction errors, and a novel approach using SHAP Values of Reconstruction Errors. 4.2.1. Anomaly Insertion Methodology During the test phase, artificial anomalies are introduced by altering population figures within selected grids. For each timestamp 𝑡, a grid 𝑔 is chosen randomly, and 𝑝𝑜𝑝𝑔,𝑡 Figure 1: The study area encompassing a total area of 5km x is modified to the maximum or minimum value seen during 5km around Shibuya Station, divided into a grid with each cell the training period, defined as: measuring 500m x 500m. {︃ max(𝑝𝑜𝑝𝑔,𝑡′ : 𝑡′ ∈ 𝑇train ), if max anomaly 𝑝𝑜𝑝𝑛𝑒𝑤 𝑔,𝑡 = 5.1. Dataset min(𝑝𝑜𝑝𝑔,𝑡′ : 𝑡′ ∈ 𝑇train ), if min anomaly This study utilizes Mobile Spatial Statistics data generated 4.2.2. Evaluation Methodology from communication records between NTT DOCOMO’s base stations and mobile devices. This data is divided into The efficacy of each detection method is quantified using mesh units across Japan in accordance with the Regional the Hits@𝑘 metric, which determines if the true anomalous Mesh standards provided by the Ministry of Internal Af- grid is among the top 𝑘 ranks based on anomaly scores. The fairs and Communications Statistics Bureau[11]. Population scores are calculated using the following equations: counts for each grid are estimated at 10-minute intervals, 𝑒𝑦𝑔 = |𝑦𝑔 − 𝑦ˆ𝑔 | (as defined in Equation 3) (8) considering factors such as number of devices accessing each base station, market share rates, residential areas, age, Table 1 Examples of Anomaly Detection Results for Randomly Sampled Grids: This table compares the detection rankings from reconstruction error and SHAP value methods across three timestamps. Ranks close to 1 indicate higher detection accuracy. (𝑘) Timestamp 𝑡 Randomly Sampled Grid 𝑔 𝑝𝑜𝑝𝑔,𝑡 𝑝𝑜𝑝𝑛𝑒𝑤 𝑔,𝑡 rank(𝑒𝑔 ) rank(𝜓𝑔 (𝑥)) 2022/10/31 19:10 5339-3588-4 1.378 -1.257 1 1 2022/10/31 10:50 5339-3574-1 1.429 2.187 9 8 2022/10/31 9:10 5339-3584-1 1.536 1.682 6 2 and gender. To ensure privacy, the data is prepared in ac- Table 2 cordance with guidelines published by NTT DOCOMO[12]. Comparison of Hits@k Metrics for Reconstruction Error and The experimental area, as shown in figure 1 , comprises 100 SHAP Methods Across All Test Instances: This table presents the grids of 500 square meters each, centered around Shibuya performance of both anomaly detection methods under condi- Station. The population data is treated as 100-dimensional tions of maximum and minimum population changes. Hits@k feature vectors and standardized for each dimension. values range from 0 to 1, where values closer to 1 indicate higher accuracy in detecting anomalies. MAX MIN 5.2. Training Phase 𝑘=1 𝑘=3 𝑘=1 𝑘=3 The training data consists of population records with a 10- Hits@𝑘reconst 0.382 0.458 0.340 0.410 minute resolution from October 17 and October 24, 2022, Hits@𝑘SHAP 0.417 0.472 0.375 0.438 totaling 288 instances (6×24×2), were prepared. A PCA model was trained with these data, setting the dimensional- ity reduction to retain 90% of the variance. and minimum population changes, is summarized in Table 2. Hits@k values range from 0 to 1, with higher values 5.3. Testing Phase indicating more effective anomaly detection. For maximum population changes, the SHAP method The test data consists of 144 population records (6×24) from demonstrates superior performance with Hits@k scores of October 31, 2022, matching the same month and day of the 0.417 at 𝑘 = 1 and 0.472 at 𝑘 = 3, exceeding the scores of week as the training data. Anomalies were inserted using the reconstruction error method, which are 0.382 at 𝑘 = 1 the method described in Section 4.2.1. For each instance, one and 0.458 at 𝑘 = 3. Similarly, in scenarios of minimum pop- grid was randomly selected, and its population count was ulation changes, the SHAP method achieves better scores replaced with either the maximum or minimum population of 0.375 at 𝑘 = 1 and 0.438 at 𝑘 = 3, outperforming the observed for that grid. A total of 288 tests were conducted reconstruction error method’s scores of 0.340 at 𝑘 = 1 and to determine if the grid with the altered population could 0.410 at 𝑘 = 3. These findings confirm the effectiveness be accurately identified. of the SHAP method in consistently identifying anomalies under varied conditions, highlighting its superiority over 5.4. Experimental Results the traditional method. This section presents the results of anomaly detection exper- iments conducted using both the traditional reconstruction 6. Conclusion error method and the proposed SHAP value method. The performance of each method is illustrated through selected This paper evaluated the performance of established recon- examples at various timestamps, as detailed in Table 1. struction error techniques and the SHAP value method for The analysis shows varying levels of detection accuracy anomaly detection in spatio-temporal population datasets. for each method, with lower ranking values indicating The study highlighted the SHAP method’s enhanced capa- higher precision in anomaly detection. Specifically, at 19:10 bility for precise anomaly identification, which is crucial on October 31, 2022, both methods accurately detected the for high-accuracy applications such as urban planning and anomaly in grid 5339-3588-4, achieving the lowest possible emergency management. The experimental datasets were rank of 1. This instance demonstrates the effectiveness of synthetically modified to include anomalies, offering a con- both approaches in scenarios where there is a substantial trolled setting to examine and contrast the performance of change in population, from 1.378 to -1.257. these methods. Future research aims to extend the appli- Conversely, at 10:50 on the same day, the anomaly in grid cation of these techniques to real-world data, particularly 5339-3574-1 was detected with lower accuracy, resulting in scenarios impacted by events like accidents, disasters, or in ranks of 9 and 8 for the reconstruction error and SHAP significant public gatherings. methods, respectively. This indicates a reduced effectiveness of both methods in detecting anomalies associated with smaller changes in population, from 1.429 to 2.187. References Moreover, at 9:10, the SHAP method outperformed the [1] K. Ochiai, M. Terada, M. Hanashima, H. Sano, Y. Usuda, reconstruction error method by more accurately identifying Detection of non-designated shelters by extracting the anomaly in grid 5339-3584-1, with a rank of 2 compared population concentrated areas after a disaster (indus- to 6. This demonstrates the SHAP method’s enhanced ability trial paper), in: Proceedings of the 30th International to detect subtle yet significant changes in population, from Conference on Advances in Geographic Information 1.536 to 1.682. Systems, 2022, pp. 1–9. A comprehensive evaluation using the Hits@k metric, [2] N. Takeishi, Shapley values of reconstruction errors which assesses performance under scenarios of maximum of pca for explaining anomaly detection, in: Proceed- ings of the International Conference on Data Mining Workshops, Beijing, China, 2019, pp. 793–798. [3] C. Xu, W. Wang, P. Liu, A genetic programming model for real-time crash prediction on freeways, IEEE Trans- actions on Intelligent Transportation Systems 14 (2013) 574–586. [4] A. Adam, E. Rivlin, I. Shimshoni, D. Reinitz, Robust real-time unusual event detection using multiple fixed- location monitors, IEEE Transactions on Pattern Anal- ysis and Machine Intelligence 30 (2008) 555–560. [5] X. Lin, A. L. Kenneth, P. R. Spence, Exploring ex- treme events on social media: A comparison of user reposting/retweeting behaviors on twitter and weibo, Computers in Human Behavior 65 (2016) 576–581. [6] Y. Zheng, H. Zhang, Y. Yu, Detecting collective anoma- lies from multiple spatio-temporal datasets across dif- ferent domains, in: Proceedings of the 23rd SIGSPA- TIAL international conference on advances in geo- graphic information systems, 2015, pp. 1–10. [7] B. Pan, Y. Zheng, D. Wilkie, C. Shahabi, Crowd sens- ing of traffic anomalies based on human mobility and social media, in: Proceedings of the 21st ACM SIGSPA- TIAL International Conference on Advances in Geo- graphic Information Systems, 2013, pp. 344–353. [8] T. Yabe, K. Tsubouchi, A. Sudo, Y. Sekimoto, A frame- work for evacuation hotspot detection after large scale disasters using location data from smartphones: case study of kumamoto earthquake, in: Proceedings of the 24th SIGSPATIAL International Conference on Ad- vances in Geographic Information Systems, 2016, pp. 1–10. [9] S. M. Lundberg, S.-I. Lee, A unified approach to in- terpreting model predictions, in: Proceedings of the 31st International Conference on Neural Information Processing Systems, 2017, pp. 4768–4777. [10] M. Terada, T. Nagata, M. Kobayashi, Population esti- mation technology for mobile spatial statistics, NTT DOCOMO Technical Journal 14 (2013). [11] Statistics Bureau of Japan, Area mesh statistics, https: //www.stat.go.jp/data/mesh/m_tuite.html, 2024. [12] NTT DOCOMO, Mobile spatial statistics guide- line, https://www.docomo.ne.jp/english/binary/ pdf/service/world/inroaming/inroaming_service/ Mobile_Kukan_Toukei_Guidelines.pdf, 2024.