2MWT

Introduction

This document describes the backend for Orikami's s2MWT biomarker. The 2 Minute Walking Test is a measure of walking capacity. Orikami's s2MWT implements this test using timestamped smartphone location data, recorded outside to obtain adequate signal accuracy. 120 seconds of location data, collected after the subject is presented an auditory start signal 🔊, is used to calculate the distance walked during those seconds. Smartphone location update frequency is generally ~0.5-2/s (February 2024).

Interface

The s2MWT backend algorithm is implemented in Python and called from the Experiment Service. It does not implement its own pulsar client, but rather communicates with the Experiment Service using a json-formatted input file and a json-formatted result printed to stdout.

Input

Expected format for json-formatted input files:

    [
        {'experiment_id': '2020-03-01-190506_S23M03', 'imuType': 'start_signal', 'data': {'timestamp': 1602243648412.589}},
        {'experiment_id': '2020-03-01-190506_S23M03', 'imuType': 'location', 'data': {'timestamp': 1602243648412.589, 'latitude': 52.525659, 'longitude': 6.0477620000000005, 'altitude': 40.102844, 'accuracy': 3.9}},
        {'experiment_id': '2020-03-01-190506_S23M03', 'imuType': 'location', 'data': {'timestamp': 1602243648428.5193, 'latitude': 52.525661, 'longitude': 6.047763, altitude': 39.027893, 'accuracy': 3.9}},
        ...
    ]

If start_signal input is omitted, calculation starts at the earliest location timestamp available.

Output

The algorithm's result is a dictionary. It is printed to stdout and contains the following key/value pairs:

• experiment_id: copied from the input data set
• timestamp: UTC at the start of calculation
• error_message: user-intelligible error description if no result could be produced. Null if no error
• debug_info: more information on the error, geared towards developers/maintainers
• img_base64: base64 encoded .SVG image. Note: may be up to ~300kB size!
• calculation_time: amount of time [seconds] spent in calculation
• version: versioning of the background code as a whole. Should be bumped whenever generic code or the code implementing one or more models changes

img_base64 is a base64 encoded .SVG image containing a plot of the raw input locations ('unfiltered') and various filtered paths, including path linearity. It also shows a bar graph comparing the distances predicted by the various models. Models that output a distance standard deviation estimate have a grey bar displayed around the end of their bar.

Then, for each 'model' (ws_1_0, ws_1_2, and possibly wheel or analogwheel):

• distance_{model}: floating point number describing the total path distance in meters
• distance_std_{model}: floating point number describing the stddev in the distance estimation (None if not produced by the model)
• r_{model}: floating point number describing track linearity after path filtering, using Pearson's r for the model's output set of horizontal coordinates
• version_{model}: versioning of the model's implementation

Use

The backend collects all relevant sensor data during the s2MWT and presents it as a batch in a json-formatted input file:

    python main.py testdata\2019-03-21-143059_fd_x5.json

Or, in order to suppress plot generation in order to preserve calculation time and/or storage space:

    python main.py testdata\2019-03-21-143059_fd_x5.json --no-plot

The result will be something like:

    {"experiment_id": "2020-03-01-190506_S23M03", "timestamp": 1687511977.9684317, "error_message": null, "debug_info": null, "distance_ws_1_0": 260.02170894138544, "r_ws_1_0": 0.37458377954337024, "distance_ws_1_2": 263.68637999486066, "distance_std_ws_1_2": 7.839612315841891, "r_ws_1_2": 0.373534151841474, "distance_ws_1_3": 265.9595661150232, "distance_std_ws_1_3": 8.821669238780448, "r_ws_1_3": 0.3730079652740305, "distance_unfiltered": 265.66713431937706, "r_unfiltered": 0.38141362574935034, "distance_wheel": 272.0, "distance_analogwheel": 271.48, "img_base64": "kwgNDg5NyAwIApMIDQyODEgMC...", "calculation_time": 0.25095486640930176}

Testing

1. (json) test files from 2MWT study data

Test json files may be created from the 2MWT study data as follows:

    python src\test\study_to_json_files.py study_folder destination_folder

So, e.g.

    python src\test\study_to_json_files.py C:\Users\Orikami\Documents\studies\2019_03 C:\Users\Orikami\Documents\GitHub\diapro-2MWT-1\algorithm\testdata

It is possible to circumvent various 'checks' on data validity. This is for debugging purposes only:

    python main.py testdata\2019-03-21-143059_fd_x5.json --no-checks

The debug_info field will contain a message saying "WARNING: --no-checks" to alert the user to to use of this debugging feature.

2. Bland-Altman plots

The output of one or more s2MWT predictor models may be compared to a reference model using Bland-Altman plots. For this, a folder containing json input files is used, possibly created using study_to_json_files.py described above. By default, analogwheel is the reference model:

    python -m src.test.blandaltmann_json_files testdata/2020

It is also possible to use another model than analogwheel as the reference model:

    python -m src.test.blandaltmann_json_files testdata/2020 -reference mymodel

3. Pearson's r

Concurrent validity of the ws_* models (i.e., the Kalman Filter models) may be established by calculating Pearson's correlation coefficients using filter output and corresponding gold standard 'analog' measurement wheel data on study data. This repository contains a zipped data set from which erroneous data was removed. Pearson's r may be calculated for this data as follows:

    python -m src.test.pearson_r_on_dataset .\src\test\pearson_r_data.zip

The output is a dictionary, mapping each model to its corresponding Pearson's coefficient. It may also contain an errors member containing a comma-separated list of experiment ids for which distance calculation failed. The dictionary is printed to stdout. Examples:

    {'ws_1_0': 0.8879036703298288, 'ws_1_2': 0.8891692123403546, 'ws_1_3': 0.936927088050846}
    {'ws_1_0': 0.88, 'ws_1_2': 0.88, 'ws_1_3': 0.93, 'errors': '2020-03-01-151121_S07M01,2020-03-01-151355_S07M02'}

4. pytest

Testers for KRN-589 (criterion validity) and KRN-446 (basic requirements on test results) are available as pytest testers, that can be run from the krane-2mwt-algorithm root folder:

    python -m pytest

Implementation

Overview

The various model implementations all follow the data flow detailed below. Multiple models may be provided with the same data set in order to compare their characteristics. Each model will yield a distance (expressed in meters), the filtered path (see below) and possibly an estimate of the confidence in that distance (expressed in meters as well, smaller is better).
If requested, a plot is created summarizing the combined results of all models.

Data flow

The overall data flow through the 2mwt backend algorithm is as follows:

1. Verify preconditions:
   • check data availability (>= 120s after optional start signal)
   • check for sufficient (horizontal) accuracy. Every single position measurement is accompanied by the smartphone's estimate of its standard deviation in latitude/longitude and is expressed in meters. 1m would be highly accurate, 5m would be a regular location measurement, 15m would be fairly inaccurate. If the median horizontal accuracy for all position data in a 2MWT would amount to >25m, this 2MWT data set will be flagged inaccurate and no distance will be calculated.
2. Exclude data points that are considered invalid:
   • multiple locations with the same timestamp;
   • location data that are near in time but far away in position (implying physically impossible velocity).
3. Select the first 120 seconds of location data after the start signal, or the first 120 seconds if no start signal is provided. (Note that timestamps should be expressed in seconds, e.g. unix timestamps.) A set of locations ordered by their timestamps forms a path.
4. Transform WGS84 latitude/longitude path to a Euclidean coordinate system with the same unit as that used for location measurement accuracy, i.e. meters.
5. Use a Kalman Filter to find an estimate of the actual path followed by the test subject's smartphone by blending the measured path with a path model (e.g., a straight line with constant velocity). Inaccurate location data (accuracy of multiple meters or more) will force the filter to stick close to the model, whereas accurate location data will more closely follow the data rather than the model.
   Some Kalman Filter models will estimate confidence in the resulting path, in the form of a standard deviation of a distance estimate expressed in meters.
6. Transform the filtered path back from Euclidean coordinates to WGS84.
7. Calculate path length using the Haversine formula across consecutive locations.

Models

The various models differ in steps 4, 5 and 6 above.

Coordinate systems

2 Euclidean coordinate systems are being used, so two sets of coordinate system transforms (steps 4 and 6) are currently in use:

• linearized latitude/longitude
  The origin is translated by subtracting the minimum values for latitude and longitude in the path. These coordinates are then multiplied by a scaling factor, currently 100000 is used. This roughly transforms to meters. Below, $x$, $y$, $latitude$ and $longitude$ are vectors (time series): $\left(\!\begin{array}{c}x \\ y\end{array}\!\right) = 100000\cdot \left(\!\begin{array}{c}latitude - \min (latitude) \\ longitude - \min (longitude)\end{array}\!\right)$
• Universal Transverse Mercator
  UTM is a 2D map projection system that divides the world into 60 zones. Within each zone, north and east are measured in meters, so this coordinate system much more accurately relates to location accuracy specified in meters than the system above.

Kalman Filters

Step 5 is the path filtering. For this step, Kalman Filtering is used. See rlabbe for a readable, thorough introduction into this topic. Its author wrote the FilterPy Python library that we use in our implementation. There are 4 main features of the Kalman Filters that we implemented:

• linear Kalman Filter We use the classical linear Kalman Filter, stemming from our linear path/constant (2D) velocity path model. There exist variants that are suited for nonlinear models, like the Extended Kalman Filter (EKF), Unscented Kalman Filter (UKF) and Cubature Kalman Filter (CKF).
• distance estimate within the filter state
  The filter state contains x and y locations and velocities ([x, vx, y, vy]) to estimate the path that was actually walked. A fifth state variable may be added that tracks accumulated distance, estimated from (x, y) or from (vx, vy). Since each state variable carries a variance on the diagonal of the covariance matrix, this way an estimated 'standard deviation' may be extracted from the distance estimation process. Note that this value will scale with the units used for location, so it basically requires using an accurate projection like UTM described above. Note that accumulated distance introduces nonlinearity in the model, (over)stretching the use of a linear Kalman Filter (see above).
• process variance $Q$
  Process variance, denoted as $Q$, describes the amount of confidence in the model. A small $Q$ expresses high confidence, so only measurements with high confidence (low measurement variance $R$) will be able to pull the filter away from the model's predictions. Note that $Q$ actually is a matrix, and in principle, the x and y elements do not need to be identical.
  Our model assumes walking in a straight line with a fixed velocity. A low $Q$ will result in lagging at the start of the experiment (difficulty in tracking walking velocity that quickly changes from 0 to around 1-2 m/s) and rounding of sharp turns/edges. Finding an appropriate $Q$ is more art than science, but highly relevant, since it profoundly influences filter behaviour.
• decimation or interpolation Location data will generally be sampled with an interval between 0.5 and 1s, occasionally even longer. An audible start signal to the user marks the beginning of 120s of measurement interval. Both start and end of the interval will generally not coincide with a location sample. We offer two different solutions: either the samples in between start and end of the measurement interval are used for Kalman filtering, or all available data is filtered and the distance travelled is interpolated at start and end.

Model versions

1.0

This model uses linearized latitude/longitude for x and y. There is no distance error estimate. It uses the decimation approach in distance calculation.
Process variance is estimated from measurement accuracy as a scalar.

$$q=2\frac{\min {(r,r_{max})}}{r_{max}}$$

Here measurement accuracy $r$, expressed as a standard deviation with unit meters, is maximized to a constant parameter $r_{max}$. So $q$ would range from 0 for $r=0$ to 2 for $r\geq r_{max}$. Matrix $Q$ is formed from $q$ for both x and y state elements.
Paradoxically, this model will tend to deviate less from the straight line assumption for more accurate measurements.

1.0.1

This model differs from version 1.0 in that it uses interpolation rather than decimation.

1.2

This model uses UTM projection for x and y, so x and y accurately describe the 2D location in meters. It does generate a distance error estimate. It uses the decimation approach in distance calculation.
Process variance is calculated like described above under model 1.0.

1.2.1

This model differs from version 1.2 in that it uses interpolation rather than decimation.

1.3.1

This model is similar to model 1.2.1 above, with two exceptions:

• process variance $Q$
  The estimated process variance no longer depends on measurement accuracy, but on estimated velocity. The concept here is that given the basic model assumption of walking in a straight line with a fixed velocity, the error between model and reality may grow larger faster if the test subject is travelling at a higher velocity, because a change of direction will have a larger impact. $q$ requires a minimum value though, otherwise the filter will never stroll from its starting location with velocity 0. The 1.3 implementation uses $q=\max {(q_{min}, v\cdot dt)}$ This has $q$ go up from parameter $q_{min}$ for $v=0$ to the estimated distance walked to the next location. Note that while this is conceptually sound, it should be noted that $q$ represents a standard deviation, so this particular implementation may overstate the model error. It should also be noted that for an even more realistic model, the x and y components could be evaluated separately based on the velocity vector.
• heuristic local velocity filtering
  All models described above filter consecutive location data for unrealistic 'jumps' in location that could signal velocity unattainable by a walking human. All models described above have these location samples removed from the data set. Model 1.3 removes this heuristic and uses the kalman filter with improved process variance to deal with noisy data. 'Never discard data' - Roger R. Labbe Jr.

Future improvements

This section describes a number of suggestions for improving 2MWT calculation from position data.

Architectural improvements

• Pulsar client A 'standalone' 2MWT Background Service could be built that includes a pulsar client and simply listens for position data measurement. It could use each and every position measurement to update a location model for every user it receives position data from. It would then generate a distance estimate after at least 120s of position data is available after a timestamped start signal. This would simplify the Experiment Service.
• Dynamic duration Duration of a walking test matters for patients with compromised physical fitness, which is why 2MWT is used in some patient groups and 6MWT in other groups. Duration could be a simple parameter, e.g. provided with the start signal, in order to implement both with a single piece of background service. Note that distinct durations may need to be validated separately because of the influence of physical fitness. Duration may also influence the error in the estimate.

Model improvements

• Velocity from position data Current smartphone OS's provide an estimate of a velocity vector in addition to position, including an error estimate of that velocity. This could easily be added to the update step in the Kalman Filter to better accomodate change of direction (corners) and change of velocity (starting or turning).
• Distance estimate as part of the state Models 1.2 and above contain a form of distance estimate as part of the state of the filter, but it is compromised by the linear nature of the current Kalman Filter. A better nonlinear estimate could simplify the calculation process - distance, rather than the estimated path, would be the outcome of the filter - and provide a reliable estimate of accuracy (variance) of the result.
• Include (linear) acceleration The current state model uses 2D position and 2D velocity. It could be enhanced by including (linear) acceleration, which would be observed by the smartphone's IMU and added to the update step in the Kalman Filter. There is one major drawback though: this data will depend on walking pattern, possibly requiring validation of the 2MWT/6MWT separately for groups of people with distinct walking patterns. (It can be argued that location data is agnostic to disease). Note that adding an observation to the filter's update step also requires an error estimate of the observation, which may be hard to obtain. There is plenty of literature describing filters that integrate raw or calibrated IMU data and GNSS data into 'fused' position measurements.
• Non-linear filters To accomodate users that do not walk along a straight line with constant velocity, a nonlinear path model could be used. The distance estimate and IMU data additions mentioned above are examples of the introduction of nonlinearity. A nonlinear model requires another class of Kalman Filter. For a small state the Unscented Kalman Filter (UKF) could be used, whereas the elaborate filters fusing IMU and GNSS data described in the literature may employ a Cubature Kalman Filter (CKF) or even more complicated filter types.

Diagrams

Below a generic overview of the 2MWT data collection system within the Krane framework.

The process of performing a 2MWT is initiated by the user and involves dataflow across various subsystems, as detailed below. Use hyperlinks for more details.

The 2mwt Service uses an elaborate communication mechanism to guarantee scalability and stability and prevent problems in case of aborted tests.

2MWT backend processing takes place in a Python environment, called from the 2mwt Service:

Links

1. 2 Minute Walking Test
2. Kalman and Bayesian Filters in Python
3. Universal Transverse Mercator