Ajaykumar Manivannan

Elias J Willemse

Wei Chien

Benny Chin

Bige Tunc ¸er Yuren Zhou

Alain Barrat

Member, IEEE Roland Bouffanais
email: roland.bouffanais@uottawa.ca

Bige Tunc

A Framework for the Identification of Human Vertical Displacement Activity Based on Multi-Sensor Data

Keywords: Wearable sensors, multi-sensor identification, human activity recognition, vertical displacement activity

With new cities increasingly expanding vertically, there is a pressing need to shed light on human vertical mobility, which can readily be achieved with existing sensor technology. To date, the methodology to track and identify vertical movement from large-scale unstructured data sets is lacking. Here, we design and develop such a framework to accurately and systematically identify the sparse human vertical displacement activity that is typically buried into the predominantly horizontal mobility. Our framework uses sensor data from barometer, accelerometer and Wi-Fi scanner coupled with an extraction step involving a combination of feature engineering and data segmentation. This methodology is subsequently integrated into a machine-learning-based classifier to automatically distinguish vertical displacement activity from its horizontal counterpart. We confirm the high accuracy of this approach by a thorough validation and testing showing a 98% overall accuracy and a 92% F1-score in classifying vertical displacement activity.We illustrate the potential of the developed framework by applying it to an unstructured large-scale data set associated with over 16,000 participants going about their daily activity in the city-state of Singapore. This gives us access to all the vertical movements of this large population, and we investigate the statistical distribution of vertical activity, both in terms of number of events and size of vertical jumps, and their temporal heterogeneity across the day. The approach developed here could be used in massive human experiments to uncover the hidden patterns of human vertical mobility. This new knowledge would have significant ramifications for the architectural design of vertical cities.

I. INTRODUCTION

T HE urbanization of our planet is rapidly increasing with 55% of the world population now living in cities [START_REF]World's cities in 2018-data booklet[END_REF]. In 2030, this number is projected to increase to 60% [START_REF]World's cities in 2018-data booklet[END_REF]. In the face of this unabated urbanization trend, cities are struggling to accommodate the population influx through urban sprawl alone, primarily because of land scarcity and the induced strain on transportation networks. An alternative to urban sprawlcurrently predominant in the rapidly urbanizing Asia-consists have developed several new frameworks to analyze and identify specific patterns of mobility. Those patterns hidden in troves of high-resolution mobility data have been uncovered thanks to large-scale experiments or from massive commercial databases-e.g., call detailed record (CDR) [START_REF] Zhou | Understanding urban human mobility through crowdsensed data[END_REF]. However, it is important to note that these human mobility studies are limited to horizontal movements, i.e. based on a two-dimensional representation of the urban landscape. In these studies, roads, railway networks, and pedestrian pathways are modeled on a planar surface.

As already mentioned, the rapid change in the topology of cities in the developing world, especially in Asia, is prominently three-dimensional [START_REF] Yeh | Introduction: High-rise living in asian cities[END_REF]. The development of sustainable and livable cities therefore heavily depends on studies of human mobility across all three dimensions. While brand-new skyscrapers are being erected every day globally, we know surprisingly little about vertical human mobility. For instance, some recent studies propose to study vertical displacements in vertically integrated mixed-use developments as a means to identify key spatial connectors that have a direct influence on social interactions [START_REF] Gopalakrishnan | User-driven emergent patterns of space use in vertically integrated urban environments[END_REF]. Moreover, a systematic study of vertical human mobility would benefit urban/infrastructure planning in many ways: e.g., targeted facility allocation in high-rise buildings, optimal placement of vertical nodes based on the typology of the building and the estimated vertical transportation load, effective vertical integration of a building in its neighborhood, etc. [START_REF] Bouffanais | The rise of homo verticalis[END_REF].

Horizontal mobility has been extensively studied using a wide variety of data sources: e.g., census data, travel surveys, CDR, location-based social network services, GPS [START_REF] Barbosa | Human mobility: Models and applications[END_REF], and smart travel/transit cards [START_REF] Hasan | Spatiotemporal patterns of urban human mobility[END_REF]. However, none of these approaches and sensors can be used effectively to track vertical displacements. Interestingly though, the sensor technology required to accurately track such human vertical mobility is readily available. What is missing is a methodology that enables the accurate identification of various types of possible vertical displacements from the output of large-scale human experiments with sufficient statistical significance.

Here, we report a contribution towards that goal by introducing and validating a methodology to accurately and systematically identify the sparse human vertical displacement activity (VDA) [START_REF] Manivannan | On the challenges and potential of using barometric sensors to track human activity[END_REF] that is deeply embedded within the predominantly horizontal displacement activity. This particular methodology is then integrated into a machine-learning-based classifier capable of dealing with large-scale data sets collected in free-living and unstructured urban environments. Classically, barometers have been the primary type of sensor used to track motion in the vertical direction. Indeed, barometric pressurepossibly augmented by other sensors-is commonly used in the field of Human Activity Recognition (HAR) to recognize the particular VDA class, which is of prime interest to us [START_REF] Manivannan | On the challenges and potential of using barometric sensors to track human activity[END_REF]. Specifically, VDA is a particular human activity class that deals with the vertical displacement of individuals in the built environment through commonly available modes of vertical mobility such as stairs, escalators, elevators, or slopes. In this work, the term VDA is intended to solely encompass human movements in vertically built structures. That means that we are discarding changes in elevation associated with any vehicle motion (motor vehicle, train, bicycle, cable car, etc.).

Our methodological advancement is thoroughly tested and validated using a big data set obtained from a large-scale human experiment carried out in Singapore: the so-called National Science Experiment (NSE). The NSE was a cityscale experiment that involved 50,000 students in Singapore between 2015 and 2017. The wearable devices specifically designed for this large-scale experiment were carried by students continuously for 5 days, and contained several sensors including a barometer and an accelerometer. By fully understanding the complex interplay of factors that influence barometric pressure, we develop several preprocessing methods to alleviate the effects of those factors, and with the end goal of achieving the highest possible accuracy in the VDA identification process. Moreover, the VDA extraction process must be robust enough to handle inherent limitations associated with such large-scale human experiments-i.e., low sampling rate, heterogeneity in devices and participant population, missing data, sensor errors, etc. As part of this process, we label manually a large number of training data (81 subjects, 81 devices, for a time period of 24 hours). This step is followed by a validation using a shortterm video-annotated data set (2 subjects, 5 devices, for a timeperiod of 6 hours). We integrate the developed VDA identification methodology into a machine-learning based classifier, and subsequently applied to the large-scale NSE data set to extract some unique features of human vertical mobility associated with the student population participating in the NSE.

The main contributions of this paper can be summarized as follows (see Fig. 1):

• A novel and accurate multi-sensory identification of vertical displacement activity is developed and validated against a sparse data set from a large-scale human experiment involving over 16, 000 individuals going about their daily activity within a densely urban environment. • The accuracy of this VDA identification process is found to be strongly dependent on a number of constraints associated with the sensing of key physical quantities. Specifically, we design a feature extraction step involving a combination of feature engineering and data segmentation. In addition, the properties of the sensors and how they are used in such large-scale experiments create a number of challenges, which are identified and addressed. • A machine-learning based classifier using our novel VDA identification process allows us to carry out the first large-scale analysis of human vertical mobility in a cityscale experiment. Interestingly, our results reveal a highly heterogeneous distribution of vertical activity, both in terms of the number of events and of the size of vertical jumps. These results have far-reaching implications for the architectural design of dense urban environments.

II. RELATED WORKS

In the field of HAR, sensor data requires pre-processing followed by the application of recognition models to classify the activity classes of interest [START_REF] Lara | A survey on human activity recognition using wearable sensors[END_REF], [START_REF] Bulling | A tutorial on human activity recognition using body-worn inertial sensors[END_REF]. Algorithms for identifying VDA from sensor data range from simple threshold-based models to sophisticated Machine Learning Fig. 1: Human activity recognition (HAR) process in the identification of VDA [START_REF] Straczkiewicz | A systematic review of human activity recognition using smartphones[END_REF]. This flowchart describes the overall process developed to achieve our VDA analysis represented in the last rectangular box. The key contributions of this work are in the boxes with a red outline.

(ML) algorithms like deep-learning. In the present work, we employ a ML-based classification model and we review here several prior studies concerning the use of ML to classify VDA events.

One of the pioneering works on identifying VDA was performed in 1998 by Sagawa et al. [START_REF] Sagawa | Classification of human moving patterns using air pressure and acceleration[END_REF], using accelerometer and barometer data to identify VDA. The study used a threshold-based model with a small training sample (83 minutes and 6 subjects). Since then, the combined effects of the digital revolution, the ubiquity of mobile devices, and the advances in sensors, Big Data and ML have paved the way for the possibility of real-time recognition of activity classes of individuals evolving in complex environments.

Accelerometers are the most widely used sensors to track human activities [START_REF] Janidarmian | A comprehensive analysis on wearable acceleration sensors in human activity recognition[END_REF], [START_REF] Twomey | A Comprehensive Study of Activity Recognition Using Accelerometers[END_REF]. In studies using an accelerometer as the stand-alone sensor, the prevalent mode of tracked vertical mobility is stairs climbing [START_REF] Lara | A survey on human activity recognition using wearable sensors[END_REF], [START_REF] Chen | Performance Analysis of Smartphone-Sensor Behavior for Human Activity Recognition[END_REF]- [START_REF] Cornacchia | A survey on activity detection and classification using wearable sensors[END_REF], although the classification accuracy of stairs climbing tends to be lower than other activity classes [START_REF] Yan | Energy-efficient continuous activity recognition on mobile phones: An activity-adaptive approach[END_REF], [START_REF] Bao | Activity Recognition from User-Annotated Acceleration Data BT -UbiComp 2002: Ubiquitous Computing[END_REF]. This led researchers to include additional data sources to improve the accuracy of VDA detection, such as gyroscope, magnetometer and barometer data [START_REF] Liu | Double-windowsbased motion recognition in multi-floor buildings assisted by a built-in barometer[END_REF]- [START_REF] Halabi | Integrating pressure and accelerometer sensing for improved activity recognition on smartphones[END_REF], and also to consider other modes of vertical mobility such as escalator and elevator rides. For instance, Liu et al. [START_REF] Liu | Double-windowsbased motion recognition in multi-floor buildings assisted by a built-in barometer[END_REF] added barometer data to a model that used accelerometer, gyroscope, and magnetometer sensor data, notably improving the classification accuracy from ∼ 80% to ∼ 90% (number of subjects: 10).

Additionally, some studies have acknowledged that accelerometers are effectively less robust than barometers for VDA classification. Muralidharan et al. [START_REF] Muralidharan | Barometric phone sensors: More hype than hope![END_REF] compared the VDA recognition performance using an accelerometer versus a barometer, and showed that their VDA classification performances were similar -with the barometer-based framework performing slightly better at nearly 100% accuracy (number of subjects: 2). However, the accuracy of accelerometer-based framework dropped drastically when the mobile device was used to take calls or play games [START_REF] Muralidharan | Barometric phone sensors: More hype than hope![END_REF]. Similarly, Vanini et al. [START_REF] Vanini | Using barometric pressure data to recognize vertical displacement activities on smartphones[END_REF] showed that the classification performance of VDA was comparable for accelerometer-only and barometer-only study (∼ 99%), but that the barometer was more energy efficient and less dependent from the on-body position than accelerometers (number of subjects: 10).

Our study focuses on the recognition of VDA as a general class of activity, and uses the following sensor data : (1) location data derived from Wi-Fi Access Points (APs), (2) magnitude of 3-axis accelerometer, and (3) barometric pressure data.

The magnitude of the 3-axis accelerometer signal is orientation independent [START_REF] Lee | Human activity recognition from accelerometer data using convolutional neural network[END_REF], [START_REF] Cleland | Collection of a Diverse, Realistic and Annotated Dataset for Wearable Activity Recognition[END_REF], and does not require complex data post-processing [START_REF] Lara | A survey on human activity recognition using wearable sensors[END_REF]. On the other hand, raw barometric sensor data can entail noise introduced by random sensor errors, limited sensor resolution, and high sampling frequency (> 2 Hz). Filtering techniques like moving average filters [START_REF] Bollmeyer | Precise indoor altitude estimation based on differential barometric sensing for wireless medical applications[END_REF]- [START_REF] Sabatini | A sensor fusion method for tracking vertical velocity and height based on inertial and barometric altimeter measurements[END_REF], Finite Impulse Response (FIR) filters [START_REF] Ghimire | Pedestrian motion state classification using pressure sensors[END_REF], and Infinite Impulse Response (IIR) filters [START_REF] Moncada-Torres | Activity classification based on inertial and barometric pressure sensors at different anatomical locations[END_REF], [START_REF] Ghimire | Pedestrian motion state classification using pressure sensors[END_REF]- [START_REF] Voleno | Energy expenditure estimation using triaxial accelerometry and barometric pressure measurement[END_REF] are commonly used to alleviate the noise effects. To increase precision in extracting elevation changes, signal modeling such as sinusoidal fitting model [START_REF] Massé | Improving activity recognition using a wearable barometric pressure sensor in mobility-impaired stroke patients[END_REF] and sigmoidal nonlinear fitting [START_REF] El Achkar | Instrumented shoes for activity classification in the elderly[END_REF] are also used. In the present study, the spectral resolution of the sensor data collection is 0.06 Hz. Such low sampling rate allows the system to side-step noise appearing at high sampling rate that particularly affects precise extraction of elevation changes. Therefore, the barometric pressure sensor data we consider here are not filtered, and other sources of noise due to sensor resolution are used to quantify the uncertainty in the magnitude of the predicted VDA. Barometric pressure data is usually converted to several common feature types such as : (1) statistical [START_REF] Elhoushi | Using portable device sensors to recognize height changing modes of motion[END_REF], [START_REF] Figueira | Body location independent activity monitoring[END_REF], (2) spectral [START_REF] Elhoushi | Using portable device sensors to recognize height changing modes of motion[END_REF], [START_REF] Figueira | Body location independent activity monitoring[END_REF], (3) temporal [START_REF] Elhoushi | Using portable device sensors to recognize height changing modes of motion[END_REF], [START_REF] Figueira | Body location independent activity monitoring[END_REF], and (4) waveletbased features [START_REF] Ejupi | Accuracy of a wavelet-based fall detection approach using an accelerometer and a barometric pressure sensor[END_REF], [START_REF] Leuenberger | Classification of stair ascent and descent in stroke patients[END_REF]. The most commonly used fea-tures are the rate of change of pressure (vertical velocity or slope) [START_REF] Liu | Double-windowsbased motion recognition in multi-floor buildings assisted by a built-in barometer[END_REF], [START_REF] Bollmeyer | Precise indoor altitude estimation based on differential barometric sensing for wireless medical applications[END_REF], [START_REF] Elhoushi | Using portable device sensors to recognize height changing modes of motion[END_REF], [START_REF] Sankaran | Using mobile phone barometer for low-power transportation context detection[END_REF], [START_REF] Monteiro | Using smartphone pressure sensors to measure vertical velocities of elevators, stairways, and drones[END_REF] and differential pressure (dp) [START_REF] Sankaran | Using mobile phone barometer for low-power transportation context detection[END_REF], [START_REF] Wu | Monitoring building door events using barometer sensor in smartphones[END_REF].

Vanini et al. [START_REF] Vanini | Using barometric pressure data to recognize vertical displacement activities on smartphones[END_REF] used barometric pressure data alone to recognize VDA using the rate of change of pressure along with the standard deviation of differential pressure as features and found that Long Short-Term Memory (LSTM) neural network framework produces a 99% accuracy compared to a decision tree approach (96%) and naïve Bayes classifiers (93%). However, their data collection was of short duration (30 minutes for each class) and conducted in limited environments. Muralidharan et al. [START_REF] Muralidharan | Barometric phone sensors: More hype than hope![END_REF] detected floor changes with an accuracy of 99% using the J48 decision tree model. Even though several factors that affect barometric pressure are considered, the data collected were of short duration (few minutes), conducted in limited structured environments, and lacked any entanglement with transportation modes. Liu et al. [START_REF] Liu | Double-windowsbased motion recognition in multi-floor buildings assisted by a built-in barometer[END_REF] classified vertical displacement activities from horizontal displacement activities (HDA) using inertial measurement units (IMU, including magnetometer), and barometer sensor data, with barometric sensor features derived from the standard deviation of pressure and rate of change of pressure. By training various classifiers such as Random Forest, J48 decision trees, Artificial Neural Networks (ANN), SVM and Naïve Bayes, they obtained that Random Forest classifiers produced the highest accuracy of 92%. Also in this study, each activity class was performed only for a few minutes and limited to ambulation [START_REF] Liu | Double-windowsbased motion recognition in multi-floor buildings assisted by a built-in barometer[END_REF].

The review of the literature on HAR [START_REF] Lara | A survey on human activity recognition using wearable sensors[END_REF], [START_REF] Bulling | A tutorial on human activity recognition using body-worn inertial sensors[END_REF], [START_REF] Wang | Deep learning for sensor-based activity recognition: A survey[END_REF]- [START_REF] Attal | Physical human activity recognition using wearable sensors[END_REF] makes it clear that no classifier can be considered as the best one in a universal way, i.e. without considering the context in which it is used. As each data set comes with its own set of distinct characteristics, the classifier working best for a particular data set and activity type might not have the best performance for a distinct problem or different circumstances (i.e., not generalizable) [START_REF] Lara | A survey on human activity recognition using wearable sensors[END_REF]. In this study, we have chosen two ensemble models (XGBoost and Random Forest), a Bayesian model (Naive-Bayes), and an instance-based nearest neighbor model (k-Nearest Neighbor) to evaluate and compare the performances of each model on our data set (see Fig. 1).

It is common in the HAR literature to use short-duration training data collected in segments that contain only one or two activity classes and are performed in semi-natural or laboratory conditions, with limited variability in environments. However, real-life activities of human occur in complex and unstructured environments, with a wide range of possible sequences, spanning heterogeneous activity classes with heterogeneous durations. Our study collocates itself in such a framework, as it uses a long-term (5 days) data set collected in a large-scale student population (∼ 50, 000 students) during their routine weekdays. Hence, it requires a different approach than those reported in the literature.

Indeed, long-term monitoring of human activities requires a thorough understanding of all the factors that affect the sensor data in different static and dynamic environments. In particular, the factors that influence barometric pressure data are climate and weather, air velocity during motion, built environment, altitude, and sensor accuracy [START_REF] Manivannan | On the challenges and potential of using barometric sensors to track human activity[END_REF]. For a detailed review on the use of barometers to track human activity and the many factors that affect barometric pressure, we refer the reader to the recent review paper [START_REF] Manivannan | On the challenges and potential of using barometric sensors to track human activity[END_REF].

III. DESCRIPTION OF DATA

A. National Science Experiment data

The National Science Experiment (NSE) was designed and commissioned by the Singapore National Research Foundation (NRF) and the Singapore University of Technology and Design (SUTD), with other private and government bodies in Singapore [START_REF] Sutd | National science experiment, singapore[END_REF]. The primary objective of this island-wide science experiment carried out by Singapore students-themed "Step Out for Science"-was for students to monitor and evaluate their own carbon footprint, travel mobility patterns, amount of time spent indoors and outdoors, and more.

Almost 50, 000 students from 92 schools distributed nationwide participated in the NSE in 2016 (Table I). Each student carried a wearable device called SENSg (see Fig. 2), which consisted of built-in environmental, motion sensors and communication units. The devices were able to record and transmit the sensed data related to the Activities of Daily Living (ADL) of the students to a cloud server [START_REF] Wilhelm | Wearable environmental sensors and infrastructure for mobile large-scale urban deployment[END_REF]. The data was recorded every 13 ∼ 18 seconds over weekdays, from Monday to Friday, for 8 different weeks in 2016. To reduce battery consumption, the device goes to a sleep mode when the processed IMU signal shows no user movement. The SENSg devices were handed over to the students on Mondays and collected back on Fridays. Hence, the full time scale of daily ADL is only available on Tuesdays, Wednesdays and Thursdays. Table II shows the cleaned NSE 2016 database after removing devices based on two criteria (a) data coverage for less than 6 hours per day (b) percentage of missing location data larger than 50% in a day. These sensor types are widely available in modern smartphones, rendering them ideal for this particular HAR. The SENSg device comprises the IMU sensor MPU9250 from InvenSense, and the barometer sensor BMP280 from Bosch Sensortech. The device also collected and stored up to a maximum of 20 Wi-Fi Access Points (AP) with the highest Receiver Signal Strength Indication (RSSI). We used Skyhook, a mobile location service from Boston, Massachusetts, that has geolocation of billions of Wi-Fi APs around the world, to convert the Wi-Fi APs to location coordinates [START_REF] Skyhook | Skyhook location services[END_REF]. The location coordinates (latitude and longitude) have a typical location accuracy of ±100 meters. The location accuracy is increased by applying regression to the time series of location data (see Appendix B). When the Wi-Fi APs are sparse or absent, location data is considered as missing. Interpolation is then applied to the time series data to predict these missing values (see Appendix A). The full sensor characteristics are shown in Table III [START_REF] Wilhelm | Wearable environmental sensors and infrastructure for mobile large-scale urban deployment[END_REF]. The embedded barometer sensor is capable of detecting up to 1 meter changes in height, i.e. ±12 Pa [START_REF] Wilhelm | Wearable environmental sensors and infrastructure for mobile large-scale urban deployment[END_REF]. The raw values measured by the accelerometer along its three axis were processed on-board the SENSg, and only descriptive statistics of these raw data were recorded:

(1) max(M acc )-the maximum value of the accelerometer's signal magnitude (2) std(M acc )-the standard deviation of accelerometer's magnitude, both sampled at 100 Hz during the one-second data acquisition temporal window, which occurred in its turn with frequency ∼ 0.0625 Hz. The SENSg device along with its working cycle are shown in Fig. 2.

B. Video-annotated data

As the measure of the barometric pressure behavior is influenced by many factors, it is important to have at our disposal a data set with the corresponding ground truth for validation purposes. To this aim, we collected approximately 6 hours of sensor data annotated by means of video recording. This data set was recorded across different modes of horizontal (walking, idle, train, bus, car, and cycle) and vertical (elevator, escalator, and stairs) activity (Table IV). Data was collected by 2 researchers on different days and time using a SENSg device that was hung using a lanyard similar to the one used by students in the NSE. The video was recorded using a Go Pro Hero 6 mounted on the chest. We will use this data set to validate manual labeling methods described in Sec. IV-A.

IV. METHODOLOGY

As the NSE was conceived as a large scale data collection with relatively high temporal resolution and for long durations, sensor data collection was optimized to save battery life and data bandwidth. This resulted in a compromise in sensor resolution, sampling rate, and type of sensor data collected, which inevitably makes the VDA identification more challenging. This section thus details the machine learning framework that we developed to identify and extract VDAs from continuous temporal segments of the NSE data (Fig. 1). Section IV-A explains the manual labeling techniques, and sections IV-B, IV-C, IV-D, and IV-E encompass the machine learning framework. We refer to the Appendices for details related to the pre-processing of the location data (App. A and B), while App. C describes the classification model parameters, and App. D deals with model tuning.

A. Manual labeling

Accurate annotation of sensor data is highly manpower intensive. One solution is to manually label a small subset of the data based on expertise, and use it to validate the classification framework. We were indeed able to perform such a manual labeling in a subset of the NSE database (see Sec. IV-D), using our general understanding of the factors affecting barometric pressure [START_REF] Manivannan | On the challenges and potential of using barometric sensors to track human activity[END_REF] and leveraging the data annotated with the help of the ground truth video (Sec III-B). In this manual labeling, we labeled each data point as either VDA or HDA, based on the unique characteristic profile of vertical transportation. As noted in [START_REF] Manivannan | On the challenges and potential of using barometric sensors to track human activity[END_REF], pressure changes caused by factors other than vertical transportation are mostly long-term variations (e.g. diurnal pressure cycle) and/or brief and intense transient spikes (e.g. Indoor-to-outdoor transition).

First, and based on our pilot experiments illustrated in [START_REF] Manivannan | On the challenges and potential of using barometric sensors to track human activity[END_REF], we distinguish VDA from HDA in our manual labeling by monitoring the following characteristics: (1) a low horizontal travel velocity, (2) body movements picked up by the accelerometer, and (3) an increase or decrease in barometric pressure. We note that high magnitude pressure changes that can be mistaken for VDA occur during transportation modes [START_REF] Manivannan | On the challenges and potential of using barometric sensors to track human activity[END_REF], but these transportation modes can be identified by tracking the location data trend, which should then reveal a high horizontal velocity.

For validation of the procedure, the classifier will be trained on the manually labeled data set, and then applied to the data set described in Sec.III-B, for which the ground truth is available.

B. Data segmentation

We segment the time series sensor data into regions of significant and minor pressure changes to allow the classifier to focus on the main characteristic of VDA-i.e., the pressurealtitude relation. First, each data point is considered to a have significant pressure change based on a cut-off value (dp i,cut-off ), determined by considering three points: (1) we want to discard small pressure jumps that could correspond to other factors that yield pressure changes of similar magnitude such as slopes, indoor-outdoor transitions, etc. [START_REF] Manivannan | On the challenges and potential of using barometric sensors to track human activity[END_REF], (2) we want to consider pressure changes corresponding to at least 50% of the minimum vertical displacement of a single floor, and (3) the changes need to be consistent with the sensor resolution (±12 Pa).

The change in pressure (dp i ) for each data point is calculated from the difference of the time series data dp i = P i -P i+1 , where P i is the pressure datum at instant t i . The time interval dt i = abs(t i -t i+1 ) associated with dp i should be less than a cut-off value dt i,cut-off to take into account that data can be missing in times of inactivity, and the resulting variation in pressure values for large values of dt i might then be due to the diurnal pressure cycle. The consecutive significant pressure changes in the same direction (positive change or negative pressure drop) are then grouped to form a segment, i.e. S Pi,n+1 = {P i , P i+1 , . . . , P n , P n+1 } for a pressure change sequence of S dpi,n+1 = {P i -P i+1 , P i+1 -P i+2 , . . . , P n -P n+1 }. Other features and sensor data are grouped using the same groups of indices as for the segmented pressure sequences.

In the case of manually labeled time series data set, we first label each data point as either VDA or HDA (Sec. IV-A). Therefore, each segmented data might contain both data points labeled VDA and data points labeled HDA. We thus label each segment using a majority rule, and assigning a VDA label in case of a draw. A perfect data segmentation would allow each VDA segment to indicate a complete VDA event with no false positives or false negatives. The choices of dt i,cut-off and dp i,cut-off ultimately determine the performance of this method. Hence, we compute the F 1 score (a.k.a. Fmeasure of balanced F -score) of capturing a complete VDA event in each segment labeled as VDA in the manually labeled training data for a range of dt i,cut-off ([30, 50, 90, 120] sec) and dp i,cut-off ([20, 23, 25, 27, 30] Pa) values (Table V). The best F 1 score is obtained for values dp i,cut-off and dt i,cut-off equal to 25 Pa and 120 seconds, respectively. We thus perform the data segmentation with these parameter values, and the final classification described below will be performed on these segmented data.

C. Feature engineering

Model explainability is a growing focus in Machine Learning. To improve explainability, it is natural to start from features based on domain-specific knowledge. We use here our data exploration and our understanding of the sensor data and the target event to be recognized to design several domainspecific features. Specifically, we compute the following features from accelerometer, barometer and Wi-Fi localization data.

Rate of pressure change dp/dt: It accounts for the pace of the VDA. This distinguishes the elevation change based activities from phenomena that unravel over slow temporal scales such as sensor drift, diurnal pressure cycle, etc.

Modified zero-crossing rate zcr: The zero-crossing rate zcr is a temporal feature that counts the number of sign changes during a particular time window for a given signal. Here, we modified this feature to count changes in sign only if the corresponding magnitude difference in pressure is ≥ 20 Pa. This conditioned zcr can indeed identify the pressure spikes due to factors such as weather and climate, built environment, air velocity during motion or sensor accuracy, and distinguish them from one-directional pressure changes that occur during changes in elevation.

Horizontal travel velocity dx/dt: as the location data are recorded in terms of latitude and longitude, we use the haversine formula (see Eq [START_REF]World's cities in 2018-data booklet[END_REF]) to calculate the great circle distance x between two locations. It is based on the assumption that the Earth is approximately spherical, a valid assumption for small distances such as the ones measured in the NSE data. The horizontal travel velocity dx/dt plays a key role in differentiating significant pressure changes of VDA from transportation based activities [START_REF] Manivannan | On the challenges and potential of using barometric sensors to track human activity[END_REF]. Specifically, the great-circle distance x between location coordinates (ϕ 1 , λ 1 ) and (ϕ 2 , λ 2 ) with ϕ the latitude, λ the longitude, and r the radius of earth, is given by:

x = 2r arcsin sin 2 ϕ 2 -ϕ 1 2 + cos ϕ 1 cos ϕ 2 sin 2 λ 2 -λ 1 2 . (1)
Statistical features of immediate neighborhood in time series data N i : By definition, a VDA event is always preceded and followed by a HDA event. However, during vertical mobility modes like elevator (and sometimes during escalator rides), a person is potentially standing with no significant body movement. Likewise, during transportation modes such as car, bus or train travel, the vehicle stops intermittently, leading to regions of low horizontal travel velocity. Hence, calculating statistical features over an immediate neighborhood of each data point can entail a sequence of events during and around an activity of interest. More precisely, for each data point i, we compute on the time window N = {i -2, i -1, i, i + 1, i + 2} (of approximate width ∼ 80 seconds) and for both the horizontal travel velocity dx/dt and the accelerometer data (max(M acc ) and std(M acc )) the ten following statistics: Minimum, Maximum, Average, Median, Mode, RMS (Root Mean Square), MAD (Median Average Deviation), Standard deviation, Variance, and IQR (Inter Quartile Range).

Statistical features of segmented data S i : Many of the statistical features considered for the immediate neighborhood of time series data are not suitable for segmented data, as the length of the sequences is typically very small for most VDA events-2 ∼ 4 s or even less (1 ∼ 2 s) for transportation modes during which the sign of slope dp/dt changes very often. Hence, we compute only the mean, median, and mode of each data segment. This is done both for the original sensor data and for the statistical features of immediate neighborhood in the time series of the modified zero-crossing rate zcr, of the horizontal travel velocity dx/dt, and of the accelerometer data (max(M acc ) and std(M acc )).

In total, we compute 95 features of the segmented data (Sec IV-B). We refer to this set of features as Feature set-I.

D. Training-Validation-Test data

The NSE 2016 data set is very diverse in terms of the number of students, unique devices and demography of individual participants. Hence, each train-validation-test data set should reflect this diversity. To ensure this, we select the data collected during one day by 81 students from 81 different schools, with an appropriate balance of school types (primary, secondary and pre-university) and weekdays (limited to Tuesday, Wednesday, and Thursday) to form a representative sample (see Table VI. We manually label this data sample according to the steps described in Section IV. This sample data set is then randomly divided into training (80%) and test set (20%), with all the data collected by any single device assigned either to training or testing as a whole (data collected by a single student cannot be split between training or testing). The classifier model is trained on the training set with the model's hyperparameters being tuned by a 5-fold cross-validation. Once the best model parameters are identified, it is then tested against the test set for the final performance evaluation.

In addition, to validate the manual labeling, we use the video-annotated data described in Sec.III-B. We apply the trained classifier model to these data labeled using the corresponding video recording to establish the ground truth.

E. Classification models

We have selected four commonly used classifiers in HAR [START_REF] Lara | A survey on human activity recognition using wearable sensors[END_REF]: (1) Extreme Gradient Boosting (XGBoost or XGB), (2) Random Forest (RF), (3) Naive-Bayes model (NB), and (4) k-Nearest Neighbors (kNN). Both XGBoost and Random Forest are decision-tree-based ensemble learning algorithms. The XGBoost algorithm is based on the boosting method that adds weak learners sequentially to reduce the loss function of the model, while the Random Forest model is based on the bagging method that adds weak learners in parallel, and uses majority voting model to make final predictions. On the other hand, the Naive-Bayes model is a probabilistic learning algorithm based on Bayes' theorem that assumes strong independence between the features. As for the k-Nearest neighbor model, it is non-parametric and uses distance-based measures to find the k-nearest samples, and it uses a majority voting model to assign a class.

The hyperparameters of these classifiers are tuned by means of a grid search using a 5-fold cross-validation on feature set-I. A more detailed description of the model parameters and model tuning can be found in App. C and D respectively.

V. RESULTS AND DISCUSSION

The central objective of our study is to obtain the magnitude of vertical displacements of individuals during their daily activities as an indication of vertical movement in a largescale study of human mobility. By leveraging the relationship between barometric pressure and altitude, our specific aim is to develop pre-processing methods that alleviate the adverse influence of the other factors that affect barometric pressure in order to accurately extract the instances of vertical motion.

We globally consider the following procedure. First, a small representative sample of the cleaned NSE 2016 data (Sec IV-D) is manually labeled as described in Sec IV-A. We then segment the whole cleaned NSE 2016 data set from Table II along with the manually labeled data and video annotated data according to Sec IV-B. The classifier model is then trained on the segmented manually labeled data to achieve a binary classification between Horizontal Displacement Activity (HDA) or Vertical Displacement Activity (VDA). The trained classifier is applied to the video annotated data to validate manual labeling, and finally applied on the cleaned NSE 2016 data set.

A. Performance comparison of classification models

The four classifiers described in sec. IV-E are trained using feature set-I with the respective hyperparameters tuned from the 5-fold cross-validation with grid search (see Sec. IV-E). Their respective classification performances on the test data set are reported in Table . VII. The XGBoost model provides the highest overall accuracy (98 %) and F 1 -score for classifying VDA (93 %), closely followed by the Random Forest (overall

B. Validation of manual labeling

To ensure the validity of our manual labeling strategies, the trained XGBoost classifier model is applied to the data set with ground truth annotated from a video recording. The results of the classifier performance is shown in Table VIII. Although the classifier has 100% precision (i.e. it does not capture false positives), it has a recall of 80% (i.e. it only captures 80% of the actual VDA events), thereby yielding an F1-score of 89%. Fig. 3: Distribution of the rate of pressure changes in predicted data in video-annotated data set when using feature set-I.

The recall performance is rather low due to a high proportion of instances such that dp/dt < 1.9 Pa/sec in the video-annotated data set (16.9% of all points marked as VDA) compared to the test data set (5.5%) for example. This issue is also responsible for some false negatives in the predicted data, where 77% of the predicted false negatives in video-annotated data set have dp/dt < 1.9 Pa/sec (see Fig. 3). This is due to the fact that many of the data points with small pressure jumps dp are not labeled as VDA in the manually labeled data set. Indeed, many instances with similar magnitudes for dp cannot be ruled out with high confidence given the known factors influencing barometric pressure [START_REF] Manivannan | On the challenges and potential of using barometric sensors to track human activity[END_REF].

C. Impact of sensor type and feature importance

Three sensors are employed for this study: Barometer, triaxis accelerometer, and Wi-Fi scanner (location data). To understand the impact of these sensor types on the classification of VDA, we conduct an ablation study. Since barometric pressure data is used in the pre-processing step (data segmentation) and is vital to the extraction of vertical displacements, features derived from barometer data are not removed in this ablation study. Table IX reports the results of the classification performance with features from the following sensor(s): (a) Barometer only, (b) Barometer and tri-axis accelerometer, and (c) Barometer and Wi-Fi scanner (location data). The F 1 -score for classifying VDA is 57% (overall accuracy of 89%) with barometer features only, and it is significantly improved by the addition of the tri-axis accelerometer features (overall accuracy of 96% and F 1 -score for classifying VDA at 83%). An even better performance is achieved with the combination of location and barometer features with 97% overall accuracy and 90% F 1 -score for classifying VDA. The addition of the accelerometer data to barometer and location data thus only improves the F 1 -score by 3%. As noted in Sec. II, previous studies [START_REF] Muralidharan | Barometric phone sensors: More hype than hope![END_REF], [START_REF] Vanini | Using barometric pressure data to recognize vertical displacement activities on smartphones[END_REF] have obtained similar accuracy with accelerometer only or barometer only features, but their training data was collected in limited environment (e.g. ambulation only or transportation only activities, short time period, etc.). However, for a data set such as the one considered here, which was collected in a long term free-living environment, it thus seems that barometer only features are not enough to classify VDA with high accuracy; however, adding the data of even only one additional sensor data allows to reach very good performances even in this case. The XGBoost model (see App. C for a description of the model and its hyperparameters) has an embedded feature ranking method that quantifies the importance of each feature to build the predictive model. The most relevant parameter to quantify relative feature importance is the total gain that measures the improvement in accuracy brought on by the feature for each tree in the model. Based on the total gain, the feature importance of all 95 features is calculated from the tuned model. To understand the impact of using reduced feature sets, we re-tune the hyperparameters and re-train the classifier model using the top 10 (out of 95) and the top 5 (out of 95) of the important features, denoted as the feature set-II (features 1-10 in Table X) and feature set-III (features 1-5 in Table X) respectively, using the same procedure as described in App. D. The reduced feature set size of 10 and 5 are arbitrarily chosen, but it is supported by the fact that the total gain drops drastically after the most important feature (S dp/dt ) as shown in found to be slightly lower at 93% for feature set-I and 92% for feature set-II and III. As the number of features are reduced from 95 to 10, and ultimately down to 5, there is a very slight increase in precision with a complementary decrease in recall.

The confusion matrix in Table XIb, shows that there is no significant drop in performance when reducing the feature set size. To strike a balance between performance and number of features, we use the XGBoost model results with feature set-II in this section's further analysis.

D. VDA classification performance

Table XIb shows the confusion matrix-i.e. predicted class distributions and corresponding Type-I (False Positives) and Type-II errors (False Negatives). When closely inspecting the feature space in the test data set, one finds that the Type-II error generally occurs when the rate of pressure changes S dp/dt is small, typically below 1.9 Pa/sec for 65% of the false negatives. This observation may have two possible origins. Either, this is due to the small magnitude of pressure jumps dp (associated to low recording frequency that splits a single VDA event across multiple time intervals). Alternatively, this could also be due to larger time intervals dt (due to irregular recording frequency of 0.076 ∼ 0.016 Hz). In both cases, the low and irregular sampling rate is the key limiting factor. In some cases, a vertical mobility event is following or preceding an horizontal transportation mode-e.g., above-theground or underground train travel-, which may lead to misclassification due to the reliance of Wi-Fi localization for the calculation of the travel velocity. This is the case for instance when the APs are sparse.

E. VDA recognition limitations

The magnitude of the altitude change in a given VDA event is derived from the barometric pressure change of the segmented data S dp using Eq [START_REF] Technavio | Elevator and escalator market by product, enduser, and geography -forecast and analysis 2020-2024[END_REF]. The accuracy of this magnitude is limited by the barometric sensor resolution, which is ±12 Pa for our SenSg device (see Table III). Some of the VDA events with lower values of S dp/dt are not appropriately classified, as stated in Sec. V-D. For a typical sampling frequency of 0.062 Hz, this corresponds to an altitude change of 2.5 meters (30 Pa). Hence, it is reasonable to assume that vertical displacements smaller than this value are not properly captured, and thus the vertical moves can only be accurate with a vertical resolution of ∼ 2.5 meters.

The classification performance when recognizing VDAs in this study is thus limited by the type of sensor data, sensor resolution, and sampling frequency. Location data from GPS with an accuracy of ±10 meters would outperform the often inaccurate Wi-Fi localization, which is only accurate within ±200 meters. Similarly, a higher sampling rate of barometric pressure to the tune of 1 Hz would be ideal compared to the lower sampling frequency of ∼ 0.062 Hz in this data set. Moreover, the use of accurately annotated data set-e.g., by means of video recording-can markedly improve the training performance.

VDA can be further sub-classified into different modes of vertical transportation like elevator, escalator, or stairs. This would however require large amounts of video-annotated training data set with high temporal and spatial resolution. This is because the rate of change of vertical displacement-a key metric in distinguishing between these modes [START_REF] Manivannan | On the challenges and potential of using barometric sensors to track human activity[END_REF]-needs to be sampled at a fairly high resolution, especially to be able to differentiate between stairs climbing and escalators riding. Furthermore, the classification between stairs/escalators and elevators will only be accurate for floor jumps larger than 2 ∼ 3 if the sampling frequency is low, such as that found in our data set. For these reasons, we do not consider VDA subclassification in this study. However, with the use of sensors with a higher sampling rate, such sub-classification of VDA should be attainable as shown in previous works [START_REF] Liu | Double-windowsbased motion recognition in multi-floor buildings assisted by a built-in barometer[END_REF], [START_REF] Muralidharan | Barometric phone sensors: More hype than hope![END_REF], [START_REF] Vanini | Using barometric pressure data to recognize vertical displacement activities on smartphones[END_REF], [START_REF] Lang | Classifying elevators and escalators in 3d pedestrian indoor navigation using foot-mounted sensors[END_REF].

F. Prevalence of vertical mobility in Singapore's student population

We have applied the trained XGBoost model with featurre set-II to the curated NSE data described in Table II. The data was collected for N = 16, 486 students from 89 schools aged 6-19 whose residences and schools were spread throughout the island-state of Singapore. For each subject, the data contains at least 6 hours of coverage during a day. The results presented in this section are aggregated over the duration of a single day for each participant-the selected day corresponds to the one with the largest among of data points collected during their week of carrying the SENSg device as part of the NSE program. The predicted VDAs are postprocessed to remove vertical displacements (less than 9% of total) that are accompanied by significant pressure fluctuations (dp > 20 Pa) that may have been caused by sensor errors (see [START_REF] Manivannan | On the challenges and potential of using barometric sensors to track human activity[END_REF]).

The total number of predicted VDAs in the entire data set reaches 182, 841 events after postprocessing. Some statistics of VDA for each subject over a day are calculated and shown in Fig. 4, namely the number of VDAs, the cummulated vertical displacement and the total time spent in the VDA mode. On a daily average, a subject was found to be engaged in 10 events of vertical mobility, traveled vertically ∼ 83 meters, and spent a total of 4 minutes per day in this mode. As an element of comparison, Americans are found to spend on average ∼ 65 minutes per day eating [START_REF] Hamrick | How much time do americans spend preparing and eating food?[END_REF]. The distributions are however rather heterogeneous, with individuals who traveled as much as 140 meters vertically in a single VDA event, moved up to 60 times in vertical direction and spent up to 24 minutes per day in VDA. The highest vertical displacement in our predicted data (140 meters) is nearly half the size of the highest building in Singapore-the Guoco tower stands at 290 meters.

Fig. 4 also shows the breakdown of the distributions of VDA for primary (PRI), secondary (SEC) and pre-university (PREU) students. Interestingly, the vertical mobility footprint of PREU students is markedly higher than that of PRI or SEC students, even if the descriptive statistics of single VDA events are very similar: the median values are 5.2 m for PRI students, 5.1 m for SEC students, 5.8 m for PREU students; the mean values are respectively 9.3 m (PRI), 8.3 m (SEC), 8.1 m (PREU) and the maximum values are 116 m (PRI), 120 m (SEC), 140 m (PREU). The fact that pre-university students are more active in terms of cumulated VDA can in particular be explained by the fact that they are more likely to travel through public transportation like trains that are either above or below the ground road level [START_REF] Monnot | How bad is selfish routing in practice?[END_REF].

The timeline of VDA reveals a rich structure. In highly vertical cities, people experience significant waiting times during their vertical mobility due to their use of shared public transportation such as elevators or slow pedestrian movements in escalators and stairs. In Singapore, these waiting times can be similar in magnitude to the waiting times between trains or buses (2 ∼ 5 minutes). Actually, businesses operating on appointment only often remind their customers to include the waiting time associated with elevator rides when planning their arrival. While these time scales can still be considered small, the ongoing trend of vertical integration of multi-purpose N is the number of subjects and x denotes the median value (dashed lines). Top row: number of VDA events per subject; Middle row: total cummulated vertical displacement (in meters) per subject; bottom row: total time spent in VDA (in seconds) per subject; left column: total population; right column: subjects grouped by school types. buildings will underscore the importance of understanding congestion in terms of vertical transportation and mobility. Understanding congestion per land use type and time is important to inform better transportation planning and this point will apply as well to vertical transportation. Fig. 5 shows the hourly statistical distribution of VDA according to the three group types under investigation (PRI, SEC and PREU). The group-based total number of subjects with VDA (i.e., active subjects, Fig. 5(a)) and total count of VDA (Fig. 5(b)) exhibit similar trends; both indicate that different daily rhythms of individuals occur between the three groups, i.e,. the peaks in activity for PRI are at 7 a.m., 10 a.m., 1 p.m., and 4 p.m., SEC at 6 a.m., 10 a.m., 2 p.m., and 6 p.m., and PREU students at 8 a.m., 12 p.m., and 5 p.m. In addition, a large number of activities of PREU students is observed during evening and night times, which can readily be explained by the fact that some schools like the Institute of Technical Education (ITE), operate primarily in the evening and by the active after-school life of the late-teen population.

Due to the low number of active users between 12 a.m. and 5 a.m., the average vertical displacement (Fig. 5(c)) and average time spent in VDA (Fig. 5(d)) shows large uncertainty within that particular time period, and become stable again after 5 a.m. In terms of per-subject average vertical displacement and time spent, the primary and secondary students show similar patterns, while pre-university students yield markedly different patterns. For instance, the primary and secondary students' average vertical displacement is high during 5 a.m.-6 a.m. time window, drops to 9 meters between 8 a.m.-12 p.m., then peaks up again to ≥ 20 meters between 4 p.m.-9 p.m., and decreases after that. On the other hand, the pre-university students average vertical displacement is about 20 meters between 6 a.m.-7 a.m., slightly decreases to 14 meters around 11 a.m., then slowly increases until 11 p.m. These differences in the morning could be explained by the style of classes between the primary/secondary and pre-university students: primary and secondary classes are usually fixed and students do not need to change classes between courses, whereas preuniversity students will need to move to other classes for different courses.

VI. CONCLUSION

Over the past two decades, advancements in sensor technology and complexity science have enabled the dissection, with unprecedented accuracy, of the fine details of human mobility in urban areas from large-scale data. Mining these "Big Data" revealed the burstiness and relative predictability of human spatial movements [START_REF] Barabási | Bursts: the hidden patterns behind everything we do, from your e-mail to bloody crusades[END_REF], [START_REF] González | Understanding individual human mobility patterns[END_REF], as well as highly dense spatial areas in cities [START_REF] Willemse | Identifying highly dense areas from raw location data[END_REF]. However, these studies, albeit illuminating, were mainly carried out on large groups of people living and moving about in mostly flat cities, and were therefore limited to horizontal movements in two dimensions. The present study explored the application of identifying and monitoring VDA to highlight the prevalence of vertical mobility at city-scale, and its potential to add to the understanding of human mobility in general.

Previous studies on VDA identification focused on the feasibility of using barometric pressure sensing, hence the associated experiments were designed for and carried out in highly controlled conditions, e.g., structured environments, short-term activities, simplified barometric pressure behavior, limited activity classes, small-population samples, and smallscale areas. In the real world, human vertical movements are typically unstructured due to the nature of human living patterns as well as the constraints imposed by the built environment. Using our understanding of the barometric pressure data behavior and considering a more general scope corresponding to real-world conditions, we have developed an accurate VDA extraction methodology, integrated into a machine learning strategy to precisely identify daily life VDA in "Big Data" sets obtained from a city-scale experiment.

This "Big Data" set consisted of 16,581 students between the age of 6 and 19, from 89 schools spread throughout the city-state of Singapore. The relatively long-term tracking (5 days) and high recording frequency (every ∼ 16 seconds) has provided high spatiotemporal resolution data sufficient to identify VDA in diverse conditions. The classification performance of the model is similar or closely aligned with the results reported in the literature (Sec. II)-with 98% overall accuracy and 92% F 1 -score in classifying VDA-while significantly generalizing its applicability to almost all possible real-world conditions.

The recognition of VDA is instrumental to the effectiveness of a range of critical applications in the fields of indoor positioning and navigation [START_REF] Muralidharan | Barometric phone sensors: More hype than hope![END_REF], [START_REF] Pipelidis | A novel approach for dynamic vertical indoor mapping through crowd-sourced smartphone sensor data[END_REF], [START_REF] Haque | A sensor fusion-based framework for floor localization[END_REF], [START_REF] Shen | Barfi: Barometer-aided wi-fi floor localization using crowdsourcing[END_REF], estimation of energy expenditure, and health monitoring [START_REF] Bianchi | Barometric pressure and triaxial accelerometry-based falls event detection[END_REF], [START_REF] Leuenberger | Classification of stair ascent and descent in stroke patients[END_REF]. The present study markedly expands the scope of applications by rendering possible the performance of large-scale human experiments with the ultimate goal of shedding light on the patterns of human mobility in dense urban areas, and its impact on urban planning.

Our study is ultimately limited by the participants age group and size, the sensor characteristics (primarily resolution and sampling rate), and the training data collection methods (video annotation versus manual labeling). Future studies can most likely achieve more significant results by improvements in one or more of these categories to obtain a more accurate and general detection of VDA.

Lastly, the analysis of the aggregated VDA data by population groups/segments has revealed rich details about the anticipated differences in the profiles of daily activity between late-teens and younger children. The uncovered patterns of vertical human mobility enabled us to accurately quantify, for the first time, the distance traveled and time spent in vertical transportation in a densely built urban environment. It has also revealed some unique patterns of activity related to vertical transportation that is present in many aspects of our human lives. There is no doubt that a more systematic analysis of the non-aggregated data would provide substantial new details and unique information about the dynamics of vertical mobility across several dimensions, including gender, age, socioeconomic status, etc.

APPENDIX

A. Localization & Interpolation

The location data is using an API from a third party company called Skyhook based on the available Wi-Fi APs. It requires a minimum number of APs to triangulate a given location. As the density of Wi-Fi APs varies from place to place, accurate localization may not always be achievable with insufficient data. This is especially pronounced in transportation modes-such as underground subway rides, by-pass roads, tunnels, bridges, etc.-where Wi-Fi APs are sparse or even non-existent, even in highly dense cities. In our dataset, a per-person average of 8% data points are missing location information due to failed localization.

To complete the Wi-Fi localization data, an interpolation method is used to predict the missing values. To identify and validate the best interpolation method, successfully localized time series are selected and a fraction of its data is removed randomly. The Root Mean Squared Error (RMSE) value between original and randomly removed data is calculated for three distinct interpolation methods: (1) linear (2) cubic spline, and (3) piece-wise cubic spline. For a fraction (10 ∼ 50%) of the data removed, the linear interpolation consistently outperforms other methods under consideration with the lowest RMSE values. For fraction> 20% of the data removed, RMSE measured in degree of latitude/longitude is ∼ 0.0001 for linear interpolation while RMSE (degree)∼ 0.005 for other methods.

B. Regression

The Wi-Fi localization data has a low spatial resolution of ±400 meters. The travel velocity estimation based on this data reveals high local errors. To reduce local fluctuations, a regression model is applied to find the best fit. To analyze the model performance, RMSE is used to estimate the residuals and the R-squared value quantifies the proportion of variance explained by the model. The regression model performs poorly when applied to the entire time series sequence. This is due to the different regions of variability in data. The time series location data is hence segmented into a series of local (temporary) variable and global (long-term) variable sequences. The regression model is then applied with different window sizes or knot placements to each of these sequences, higher number of knots for a global variable and lower number of knots for a local variable.

A local variable sequence is defined as S i,n = {i, i + 1, . . . , n} where the location of index -n is at distance ≤ distance cut-off from index -i. A global variable sequence is created by combining several local variable sequences of length < 10 (∼ 2.5 minutes) to reflect long term changes in location.

Two regression models: (1) piece-wise polynomial and (2) natural piece-wise cubic regression spline are considered as well as three smoothing models: (3) Savitzky-Golay smoothing, (4) LOESS model, (5) Exponential smoothing model. These five options are compared against different window sizes or knot placements (local variable = [START_REF] Gopalakrishnan | User-driven emergent patterns of space use in vertically integrated urban environments[END_REF][START_REF] Twomey | A Comprehensive Study of Activity Recognition Using Accelerometers[END_REF][START_REF] Bollmeyer | Precise indoor altitude estimation based on differential barometric sensing for wireless medical applications[END_REF], global variable= [START_REF] Gopalakrishnan | User-driven emergent patterns of space use in vertically integrated urban environments[END_REF][START_REF] Twomey | A Comprehensive Study of Activity Recognition Using Accelerometers[END_REF][START_REF] Bollmeyer | Precise indoor altitude estimation based on differential barometric sensing for wireless medical applications[END_REF]) and distance cut off = [0.1, 0.3, 0.5] km). The natural piece-wise Cubic regression spline and Savitzky-Golay show the lowest RMSE and highest R-Squared values for all knot placements/window sizes and distance cut off. Either of these methods will suitably reduce the local errors in location data. A low RMSE and high R-squared value can also point to an over-fitted data, hence, model selection should be supported by additional considerations. The optimal knot placements should be sparse for local variable data to reduce local fluctuations and denser for global variable data where location varies long-term. Here, Natural Piecewise cubic regression spline method is selected with knot placements at (local variable, global variable)= (15, 1) of the segmented data using distance cut off = 0.3 km. Essentially, the knot placements are designed to smooth the local variable data while leaving the global variable sequences intact.

C. Classification model parameters 1) XGBoost: Gradient boosting allows an ensemble of weak learners to build models that depend on the gradient descent algorithm to optimize an objective function. The objective function measures the model fitness of the training data and consists of training loss and regularization [START_REF]Xgboost documentation[END_REF]. Training loss quantifies the predictive capability of the model and the regularization parameters helps reduce over-fitting by controlling the complexity of the model. The balance between these two terms is commonly known as the biasvariance trade-off [START_REF]Xgboost documentation[END_REF]. XGBoost classifier model encompass several parameters that are categorized as general, booster, and learning task parameters [START_REF] Jain | Xgboost parameter tuning guide to parameter tuning in xgboost[END_REF]. The tree booster parameters allows the specification of learning rate (eta), minimum loss reduction to make a split (gamma), maximum depth of a tree (max depth), minimum sum of weights of all observations in a child (min child weight), fraction of observations to be sampled randomly (subsample), fraction of columns to be sampled randomly in each tree (colsample bytree), control of class imbalance (scale pos weight), L2 regularization term on weights (reg alpha), and L1 regularization term on weights (reg lambda) [START_REF]Xgboost documentation[END_REF], [START_REF] Jain | Xgboost parameter tuning guide to parameter tuning in xgboost[END_REF]. Learning task parameters gives the option to choose an objective function (objective) for optimization against a specified evaluation metric (eval metric). The number of trees or estimators (n estimators) required depends on the learning rate eta. The optimization of the tree booster parameters based on the evaluation metric allows the objective function to look for values that avoids over-fitting. Model over-fitting can be further controlled by reducing the learning rate and increasing the number of estimators [START_REF] Jain | Xgboost parameter tuning guide to parameter tuning in xgboost[END_REF].

2) Random Forest: The Random Forest classifier is an ensemble learning algorithm that tries to create a range of uncorrelated trees by randomly selecting features for each tree and randomly (with replacement) assigning training data to each tree. Both of these properties of Random Forest helps control over-fitting [START_REF] Yiu | Understanding random forest[END_REF], [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF]. The number of trees are given by the parameter n estimators and the maximum number of features considered for splitting a node is given by max f eatures. The maximum number of levels in each tree is controlled by the parameter max depth. The samples given to the tree are controlled by the following parameters: bootstrap, min sample split, min weight f raction leaf , and min sample leaf . The parameter bootstrap can be set to T rue to allow random sampling with replacement. The minimum number of data points assigned to a node before splitting is given by min sample split and the parameter min sample leaf is used to control the minimum number of data points allowed in a leaf node. The sample weight can be adjusted using the parameter min weight f raction leaf . Finally, the class imbalance in a data set can be controlled by setting the parameter class weight to balanced.

3) Naive-Bayes: Based on Bayes' theorem, the Naive-Bayes model is a probabilistic learning algorithm. This study uses the Gaussian Naive-Bayes algorithm that assumes Gaussian distribution for each class [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF], [START_REF] Sharma | Gaussian naive bayes with hyperparameter tuning[END_REF]. Other Naive-Bayes models such as Multinomial and Bernoulli are defined for discrete data values and hence not suitable for our problem [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF]. Gaussian Naive-Bayes model has only two parameters: var smoothing and priors. The prior probability of the classes can be assigned through priors. Since there are no prior probability available for our classes, this parameter is unspecified in this work so that the priors can be learnt from the data. The parameter var smoothing is used to adjust the weight given to data points far from the mean distribution. This is the only parameter tuned for the Naive-Bayes model here.

4) k-Nearest Neighbors: k-Nearest Neighbors classifier is a class of Nearest Neighbors algorithm that identifies k-nearest training data points based on their similarity [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF]. The number of neighbors k is given by the parameter n neighbors. The neighbors can be weighed uniformly or differently using the parameter weight. The distance metric used can be controlled by the parameter p, which denotes the power of the Minkowski distance (p = 1 denotes Euclidean and p = 2 denotes Manhattan). The speed of finding the nearest neighbors depends on the parameter algorithm that can be set to brute-force or more faster methods such as tree-based search algorithms. The tree-based search can be controlled by the parameter leaf size that can be adjusted for faster construction and queries.

D. Classifier model tuning

The four classifier models follow an overall procedural structure for model tuning. The optimal model hyperparameters are selected by grid search and a 5-fold cross-validation using the metric -Area Under the Receiver Operating Characteristic Curve (ROC AUC). The grid search is initialized with a set of hyperparameter values, which are fine tuned until no further change in the results is noticed. The 5-fold cross-validation is used to reduce the bias and makes use of the training set efficiently to understand the model's predictive power on new data during the tuning procedure. Due to the specific nature of the classifier models, some differences in the procedure exist. For example, all the hyperparameters for the Random Forest, Naive-Bayes, and k-Nearest Neighbors are tuned using a single grid search while the hyperparameters for XGBoost are tuned in a series of grid searches. The XGBoost model is thus tuned by the following steps:

XGBoost model tuning steps: The optimal parameter values are calculated in a series of steps, where each step progresses by estimating the parameter under consideration based on the parameters calculated from preceding step. The steps taken to tune the model are as follows: (1) Set a high learning rate and find the optimal number of estimators, (2) For the given learning rate and number of estimators, find optimal max depth and min child weight, (3) Find optimal value of gamma (4) Re-calibrate the optimal number of estimators, (5) find optimal values for the parameters subsample and colsample bytree, (6) find the optimal values for the regularization parameters reg alpha and reg lambda, (7) Recalibrate the optimal number of estimators, (8) Reduce the learning rate and find the optimal number of estimators [START_REF] Jain | Xgboost parameter tuning guide to parameter tuning in xgboost[END_REF].

The cross-entropy loss for binary classification is set as the objective function. For an initial learning rate of 0.3, the above steps are followed from (1) to [START_REF] Barbosa | Human mobility: Models and applications[END_REF]. As the final step [START_REF] Hasan | Spatiotemporal patterns of urban human mobility[END_REF], the learning rate is reduced considerably (from 0.3 to 0.2 for feature set-I, to 0.1 for both feature set-II, and feature set-III) to control over-fitting while keeping the F 1 score of VDA classification higher. The optimal number of estimators is then found to be 43, 125, and 72 for feature set-I, II, and III respectively from 5-fold cross validation for the given learning rate. The class imbalance is acknowledged by setting the parameter scale pos weight to 1 for faster convergence.

Fig. 2 :

 2 Fig. 2: NSE SENSg device details: (a) outside look, (b) internal structure [49], (c) working cycle. (Picture Courtesy: [50])

Fig. 4 :

 4 Fig. 4: Aggregated distribution of VDA for a population of N = 16, 486 individuals. The aggregation period is set to 1 day.N is the number of subjects and x denotes the median value (dashed lines). Top row: number of VDA events per subject; Middle row: total cummulated vertical displacement (in meters) per subject; bottom row: total time spent in VDA (in seconds) per subject; left column: total population; right column: subjects grouped by school types.

Fig. 5 :

 5 Fig. 5: Temporal statistics of VDA according to group types: (a) number of subjects with VDA (active subjects) by hours, (b) total VDA count by hour, (c) average vertical displacement per active subject by hour, and (d) hourly average time spent in VDA per active subject by hour. The errorbars in (c) and (d) indicate the confidence intervals estimated using a bootstrap approach (n = 1000).

TABLE I :

 I NSE 2016 -Device delivery information

	Dates	No.	No.	Primary	Secondary	Pre-
		schools	devices	(6-12	(12-16	University
				y.o.)	y.o.)	(16-19
						y.o.)
	Apr 11-15	7	3,002	1,590	1,192	220
	Apr 18-22	7	3,290	2,300	990	-
	May 16-20	10	2,673	820	1,624	229
	May 23-27	2	840	-	220	620
	July 11-15	32	10,751	3,834	6,677	240
	July 18-22	23	7,934	5,233	2,181	520
	July 25-29	8	15,029	220	890	13,919
	Aug 17-19	3	6,000	-	-	6,000
	Total	92	49,519	13,997	13,774	21,748

TABLE II :

 II Curated NSE 2016 database

	Dates	No.	No. IDs	Primary	Secondary	Pre-
		schools		(6-12	(12-16	University
				y.o.)	y.o.)	(16-19
						y.o.)
	Apr 11-15	7	1,797	1,049	608	140
	Apr 18-22	7	1,793	1,265	528	-
	May 16-20	10	700	239	386	75
	May 23-27	-	-	-	-	-
	July 11-15	33	4,387	2,048	2,114	73
	July 18-22	21	3,086	2,292	666	128
	July 25-29	8	3,345	-	93	3,057
	Aug 17-19	3	1,380	-	-	1,380
	Total	89	16,581	6,939	4,395	4,853

TABLE III :

 III Sensor Characteristics of SENSg device

	Sensor	Model	Range		Accuracy Units	Poll fre-
							quency
	Wi-Fi	SN8205	-		-	-	0.062 Hz
	Barometer	BMP280 300	to	±12 Pa Pa	0.062 Hz
			1, 100 hPa
	Accelerometer MPU9250 ±2 g		±80 mg mg/LSB 100 Hz
							(for every
							0.062 Hz)

TABLE IV :

 IV Video-annotated data

	Mode	Total
		time
		(hours)
	Car	0.4
	Bus	1.5
	Train	0.7
	Cycle	0.8
	Walking	1.5
	Idle	0.3
	VDA	0.4
	Total	5.6

TABLE V :

 V Selection of data segmentation parameters based on F 1 score.

	dt i,cut-off			F 1 score
				(%)
	120	97.68	98.09	98.15	96.98	96.98
	90	97.65	98.06	98.12	97.94	96.95
	60	97.65	97.99	98.06	97.87	96.88
	30	97.14	97.39	97.45	97.23	96.22
	dp i,cut-off 20	23	25	27	30

TABLE VI :

 VI Properties of the Training-Test split data sets

	data set	Ground	Size	in	Number of	Number of	Number of	Number of	Number of	Pri.	Sec.	Pre-
		truth	duration		time-series	segmented	devices	individuals	schools			U.
					data points	data
	Training	Manual	64 (days)	207,380	7,132	64	64	64	21	20	23
	(80%)
	Test	Manual	17 (days)	55,092	2,166	17	17	17	8	2	7
	(20%)
	Video	Video	6 (hours)	1,262	305	5	2	-	-	-	-
	annotated	recording
	data

accuracy of 97 % and F 1 -score for classifying VDA at 89 %). Both the Naive-Bayes model and k-Nearest Neighbor model perform poorly with an F 1 -score for classifying VDA standing at 52 % and 17 % respectively. The ensemble learning methods clearly show a superior performance than Naive-Bayes and kNN models, similar to the results reported by Liu et al.

[START_REF] Liu | Double-windowsbased motion recognition in multi-floor buildings assisted by a built-in barometer[END_REF]

in classifying VDA from HDA. Therefore, we select XGBoost as the classifier of choice for this analysis.

TABLE VII :

 VII Comparison of classifier performance on the test data set with feature set-I

	Model	Overall	Class	Precision	Recall	F 1 -	Support
		accuracy		(%)	(%)	score
		(%)				(%)
	XGB	98	HDA	99	99	99	1899
			VDA	93	92	93	277
	RF	97	HDA	98	99	99	1899
			VDA	96	83	89	277
	NB	80	HDA	97	80	88	1899
			VDA	38	83	52	277
	kNN	87	HDA	88	99	93	1899
			VDA	52	10	17	277

TABLE VIII :

 VIII Classification results of video-annotated data set with feature set-I

		(a) Performance metrics
	Class	Precision (%)	Recall (%)	F1-score (%)	Support
	HDA	95	100	97	246
	VDA	100	80	89	59
	weighted average	96	96	96	395
		(b) Confusion matrix
	Class	HDA	VDA	Support
	HDA	246	0	246
	VDA	12	47	59

TABLE IX :

 IX Impact of sensor type on the classifier performance with selected features of feature set-I. Baro -Barometer, Acc -Accelerometer, and Loc -Location data derived from Wi-Fi scanner.

	Sensor	Overall	Class	Precision	Recall	F 1 -	Support
		accuracy		(%)	(%)	score
		(%)				(%)
	Baro.	89	HDA	94	94	94	1899
			VDA	58	56	57	277
	Baro +	96	HDA	98	97	98	1899
	Acc
			VDA	83	84	83	277
	Baro +	97	HDA	98	99	99	1899
	Loc
			VDA	92	88	90	277
	Baro +	98	HDA	99	99	99	1899
	Loc +
	Acc
			VDA	93	92	93	277

Table X .

 X The classification performance of the XGBoost model with feature sets I, II, and III are shown in Table XI. The overall classification accuracy of the model for all feature sets on the test set stands at 98%, and the F 1 -score for classifying VDA is

TABLE X :

 X Feature importance obtained from the tuned XG-Boost classifier model. Top 10 results are shown. S stands for segmented data.

	No Feature	Total gain
	1	S dp/dt	4460
	2	S min(N RMS(std(Macc)) )	2614
	3	S min(N median(dx/dt) )	988
	4	S min( zcr)	421
	5	S dx/dt	265
	6	S min(N max(Macc) )	184
	7	S min(N max(dx/dt) )	184
	8	S min(N avg(dx/dt) )	141
	9	S max(N max(dx/dt) )	124
	10	S min(N min(dx/dt) )	114

TABLE XI :

 XI Classification performance of the XGBoost model

	Feature	Overall	Class	Precision	Recall	F 1 -	Support
	set	accuracy		(%)		(%)	score
		(%)					(%)
	Feature	98	HDA	99		99	99	1899
	set-I
			VDA	93		92	93	277
	Feature	98	HDA	99		99	99	1899
	set-II
			VDA	93		90	92	277
	Feature	98	HDA	98		99	99	1899
	set-III
			VDA	94		90	92	277
			(b) Confusion matrix
		Feature	Class	HDA	VDA	Support
		set
		Feature	HDA	1870	19		1899
		set-I
			VDA	22	255	277
		Feature	HDA	1870	19		1899
		set-II
			VDA	27	250	277
		Feature	HDA	1874	15		1899
		set-III
			VDA	29	248	277

(a) Performance metrics

The NSE 2016 was supported by the Singapore National Research Foundation (NRF), Ministry of Education (MOE),and the SUTD-MIT International Design Centre.

E. Pressure-altitude relation

The fundamental equation for fluids at rest dictates the relationship between atmospheric pressure and altitude [START_REF] Young | A Brief Introduction to Fluid Mechanics[END_REF]. The change in altitude is indeed given by

where γ = ρg is the specific weight of air with density ρ = 1.225 kg/m 3 and g is the acceleration due to gravity g = 9.81 m/s 2 at standard sea-level conditions. The fluid is assumed to be incompressible and in isothermal condition. This pressure-altitude equation is valid for data collected for elevation less than 10 km from sea-level [START_REF] Young | A Brief Introduction to Fluid Mechanics[END_REF], which is always the case in our study.

ACKNOWLEDGMENT

The NSE 2016 was supported by the Singapore National Research Foundation (NRF), Ministry of Education (MOE), and the SUTD-MIT International Design Centre. On-time production of the sensors was only possible due to strong support from Delta Electronics DRC, IABG, and Taoyuan Factory. Logistics for the distribution of thousands of sensors was facilitated by the Singapore Science Center, and the infrastructure to securely store the data was built and maintained by the A*STAR Institute for High Performance Computing.

We also would like to thank Francisco Benita, Francesco Scandola, and Garvit Bansal for their assistance with SENSg devices and guidance on data collection. A.M. was supported by an MOE-SUTD PhD fellowship and currently supported by University of Ottawa PhD Admission scholarship.