Gait analysis

This tutorial showcases the high-level functions composing the gait pipeline. Before following along, make sure all data preparation steps have been followed in the data preparation tutorial.

To run the complete gait pipeline, a prerequisite is to have both accelerometer and gyroscope data, although a small part of the pipeline requires only accelerometer data. Roughly, the pipeline can be split into seven segments:

  1. Data preprocessing

  2. Gait feature extraction

  3. Gait detection

  4. Arm activity feature extraction

  5. Filtering gait

  6. Arm swing quantification

  7. Aggregation

Using only accelerometer data, the first three steps can be completed.

[!WARNING] The gait pipeline has been developed on data of the Gait Up Physilog 4, and is currently being validated on the Verily Study Watch. Different sensors and positions on the wrist may affect outcomes.

Throughout the tutorial, a small segment of data from a participant of the Personalized Parkinson Project is used to demonstrate the functionalities.

Load data

Load the prepared data into memory. For example, the following functions can be used depending on the file extension of the data:

We use the interally developed TSDF (documentation) to load and store data [1].

from pathlib import Path
from paradigma.util import load_tsdf_dataframe

# Set the path to the data file location
path_to_data = Path('../../tests/data')
path_to_prepared_data =  path_to_data / '1.prepared_data' / 'imu'

# Load the data from the file
df_imu, _, _ = load_tsdf_dataframe(path_to_prepared_data, prefix='IMU')

df_imu
time accelerometer_x accelerometer_y accelerometer_z gyroscope_x gyroscope_y gyroscope_z
0 0.00000 0.550718 0.574163 -0.273684 -115.670732 32.012195 -26.097561
1 0.01004 0.535885 0.623445 -0.254545 -110.609757 34.634146 -24.695122
2 0.02008 0.504306 0.651675 -0.251675 -103.231708 36.768293 -22.926829
3 0.03012 0.488517 0.686603 -0.265550 -96.280488 38.719512 -21.158537
4 0.04016 0.494258 0.725359 -0.278469 -92.560976 41.280488 -20.304878
... ... ... ... ... ... ... ...
72942 730.74468 0.234928 -0.516268 -0.802871 0.975610 -2.256098 2.256098
72943 730.75472 0.245455 -0.514354 -0.806699 0.304878 -1.707317 1.768293
72944 730.76476 0.243541 -0.511005 -0.807177 0.304878 -1.585366 1.890244
72945 730.77480 0.240191 -0.514354 -0.808134 0.000000 -1.280488 1.585366
72946 730.78484 0.243541 -0.511005 -0.808134 -0.060976 -1.036585 1.219512

72947 rows × 7 columns

Step 1: Preprocess data

The single function preprocess_imu_data in the cell below runs all necessary preprocessing steps. It requires the loaded dataframe, a configuration object config specifying parameters used for preprocessing, and a selection of sensors. For the sensors, options include 'accelerometer', 'gyroscope', or 'both'.

The function preprocess_imu_data processes the data as follows:

  1. Resample the data to ensure uniformly distributed sampling rate

  2. Apply filtering to separate the gravity component from the accelerometer

from paradigma.config import IMUConfig
from paradigma.preprocessing import preprocess_imu_data

config = IMUConfig()

df_preprocessed = preprocess_imu_data(
    df=df_imu, 
    config=config,
    sensor='both',
    watch_side='left',
)

print(f"The dataset of {df_preprocessed.shape[0] / config.sampling_frequency} seconds is automatically resampled to {config.sampling_frequency} Hz.")
df_preprocessed.head()

The dataset of 730.79 seconds is automatically resampled to 100 Hz.
time accelerometer_x accelerometer_y accelerometer_z gyroscope_x gyroscope_y gyroscope_z accelerometer_x_grav accelerometer_y_grav accelerometer_z_grav
0 0.00 0.053078 0.010040 -0.273154 -115.670732 32.012195 -26.097561 0.497639 0.564123 -0.000530
1 0.01 0.038337 0.058802 -0.256899 -110.636301 34.624710 -24.701537 0.497666 0.564510 0.002305
2 0.02 0.006824 0.086559 -0.256739 -103.292766 36.753000 -22.942002 0.497698 0.564887 0.005122
3 0.03 -0.009156 0.120855 -0.273280 -96.349062 38.692931 -21.175227 0.497733 0.565254 0.007919
4 0.04 -0.003770 0.159316 -0.289007 -92.585735 41.237328 -20.311531 0.497772 0.565610 0.010696

The resulting dataframe shown above contains uniformly distributed timestamps with corresponding accelerometer and gyroscope values. Note the for accelerometer values, the following notation is used:

  • accelerometer_x: the accelerometer signal after filtering out the gravitational component

  • accelerometer_x_grav: the gravitational component of the accelerometer signal

The accelerometer data is retained and used to compute gravity-related features for the classification tasks, because the gravity is informative of the position of the arm.

Step 2: Extract gait features

With the data uniformly resampled and the gravitional component separated from the accelerometer signal, features can be extracted from the time series data. This step does not require gyroscope data. To extract the features, the pipeline executes the following steps:

  • Use overlapping windows to group timestamps

  • Extract temporal features

  • Use Fast Fourier Transform the transform the windowed data into the spectral domain

  • Extract spectral features

  • Combine both temporal and spectral features into a final dataframe

These steps are encapsulated in extract_gait_features (documentation can be found here).

from paradigma.config import GaitConfig
from paradigma.pipelines.gait_pipeline import extract_gait_features

config = GaitConfig(step='gait')

df_gait = extract_gait_features(
    df=df_preprocessed, 
    config=config
)

print(f"A total of {df_gait.shape[1]-1} features have been extracted from {df_gait.shape[0]} {config.window_length_s}-second windows with {config.window_length_s-config.window_step_length_s} seconds overlap.")
df_gait.head()

A total of 34 features have been extracted from 725 6-second windows with 5 seconds overlap.
time accelerometer_x_grav_mean accelerometer_y_grav_mean accelerometer_z_grav_mean accelerometer_x_grav_std accelerometer_y_grav_std accelerometer_z_grav_std accelerometer_std_norm accelerometer_x_power_below_gait accelerometer_y_power_below_gait ... accelerometer_mfcc_3 accelerometer_mfcc_4 accelerometer_mfcc_5 accelerometer_mfcc_6 accelerometer_mfcc_7 accelerometer_mfcc_8 accelerometer_mfcc_9 accelerometer_mfcc_10 accelerometer_mfcc_11 accelerometer_mfcc_12
0 0.0 0.527557 0.485994 -0.314004 0.045376 0.061073 0.312533 0.187336 0.002437 0.038174 ... 0.622230 0.700691 0.150885 0.458676 0.033595 0.243243 0.178762 0.066051 0.090567 0.128075
1 1.0 0.532392 0.466982 -0.426883 0.045834 0.052777 0.261492 0.196270 0.001498 0.012917 ... 0.288084 0.581373 0.086737 0.447655 -0.129908 0.292478 0.028656 0.123552 0.049143 0.081229
2 2.0 0.545189 0.433756 -0.539777 0.057979 0.080084 0.145126 0.200461 0.001840 0.001876 ... 0.277590 0.481786 0.069445 0.331342 -0.197669 0.244441 -0.151760 0.091923 -0.100824 0.003800
3 3.0 0.556586 0.397208 -0.613691 0.070765 0.122509 0.054681 0.104950 0.003253 0.001985 ... 0.521323 0.321901 0.237257 0.078296 -0.074240 0.283644 -0.239699 0.028845 -0.050743 0.036805
4 4.0 0.571852 0.359068 -0.639196 0.079765 0.144845 0.042924 0.095547 0.002794 0.002084 ... 0.538269 0.111283 0.293136 -0.069686 -0.059406 0.356973 -0.266953 0.050041 0.058963 0.082503

5 rows × 35 columns

Each row in this dataframe corresponds to a single window, with the window length and overlap set in the config object. Note that the time column has a 1-second interval instead of the 10-millisecond interval before, as it now represents the starting time of the window.

Step 3: Gait detection

For classification, ParaDigMa uses so-called Classifier Packages which contain a classifier, classification threshold, and a feature scaler as attributes. The classifier is a random forest trained on a dataset of people with PD performing a wide range of activities in free-living conditions: The Parkinson@Home Validation Study. The classification threshold was set to limit the amount of false-positive predictions in the original study, i.e., to limit non-gait to be predicted as gait. The classification threshold can be changed by setting clf_package.threshold to a different float value. The feature scaler was similarly fitted on the original dataset, ensuring the features are within expected confined spaces to make reliable predictions.

from importlib.resources import files
from paradigma.classification import ClassifierPackage
from paradigma.pipelines.gait_pipeline import detect_gait

# Set the path to the classifier package
classifier_package_filename = 'gait_detection_clf_package.pkl'
full_path_to_classifier_package = files('paradigma') / 'assets' / classifier_package_filename

# Load the classifier package
clf_package = ClassifierPackage.load(full_path_to_classifier_package)

# Detecting gait returns the probability of gait for each window, which is concatenated to
# the original dataframe
df_gait['pred_gait_proba'] = detect_gait(
    df=df_gait,
    clf_package=clf_package
)

n_windows = df_gait.shape[0]
n_predictions_gait = df_gait.loc[df_gait['pred_gait_proba'] >= clf_package.threshold].shape[0]
perc_predictions_gait = round(100 * n_predictions_gait / n_windows, 1)
n_predictions_non_gait = df_gait.loc[df_gait['pred_gait_proba'] < clf_package.threshold].shape[0]
perc_predictions_non_gait = round(100 * n_predictions_non_gait / n_windows, 1)

print(f"Out of {n_windows} windows, {n_predictions_gait} ({perc_predictions_gait}%) were predicted as gait, and {n_predictions_non_gait} ({perc_predictions_non_gait}%) as non-gait.")

# Only the time and the predicted gait probability are shown, but the dataframe also contains
# the extracted features
df_gait[['time', 'pred_gait_proba']].head()

Out of 725 windows, 53 (7.3%) were predicted as gait, and 672 (92.7%) as non-gait.
time pred_gait_proba
0 0.0 0.093240
1 1.0 0.093032
2 2.0 0.107607
3 3.0 0.132656
4 4.0 0.142432

Once again, the time column indicates the start time of the window. Therefore, it can be observed that probabilities are predicted of overlapping windows, and not of individual timestamps. The function merge_timestamps_with_predictions can be used to retrieve predicted probabilities per timestamp by aggregating the predicted probabilities of overlapping windows. This function is included in the next step.

Step 4: Arm activity feature extraction

The extraction of arm swing features is similar to the extraction of gait features, but we use a different window length and step length (config.window_length_s, config.window_step_length_s) to distinguish between gait segments with and without other arm activities. Therefore, the following steps are conducted sequentially by extract_arm_activity_features:

  • Start with the preprocessed data of step 1

  • Merge the gait predictions into the preprocessed data

  • Discard predicted non-gait activities

  • Create windows of the time series data and extract features

from paradigma.pipelines.gait_pipeline import extract_arm_activity_features

config = GaitConfig(step='arm_activity')

df_arm_activity = extract_arm_activity_features(
    df_timestamps=df_preprocessed, 
    df_predictions=df_gait,
    config=config,
    threshold=clf_package.threshold
)

print(f"A total of {df_arm_activity.shape[1]-1} features have been extracted from {df_arm_activity.shape[0]} {config.window_length_s}-second windows with {config.window_length_s-config.window_step_length_s} seconds overlap.")
df_arm_activity.head()

A total of 61 features have been extracted from 50 3-second windows with 2.25 seconds overlap.
time accelerometer_x_grav_mean accelerometer_y_grav_mean accelerometer_z_grav_mean accelerometer_x_grav_std accelerometer_y_grav_std accelerometer_z_grav_std accelerometer_std_norm accelerometer_x_power_below_gait accelerometer_y_power_below_gait ... gyroscope_mfcc_3 gyroscope_mfcc_4 gyroscope_mfcc_5 gyroscope_mfcc_6 gyroscope_mfcc_7 gyroscope_mfcc_8 gyroscope_mfcc_9 gyroscope_mfcc_10 gyroscope_mfcc_11 gyroscope_mfcc_12
0 484.00 0.977134 0.105431 0.045578 0.015722 0.026817 0.015546 0.085454 0.000284 0.000632 ... 0.042879 0.054217 0.275817 0.463253 -0.061382 0.144563 0.243984 -0.174514 -0.254318 0.118443
1 484.75 0.981950 0.102169 0.041451 0.006373 0.019288 0.018058 0.084179 0.000445 0.000663 ... -0.139456 -0.195451 0.357785 0.270232 -0.250594 -0.216808 -0.085842 0.197042 0.004748 0.044400
2 485.50 0.978526 0.115141 0.030499 0.006440 0.022474 0.015048 0.089055 0.000328 0.000408 ... -0.068866 -0.180071 0.474648 0.314642 -0.180227 -0.348940 -0.157215 0.166344 0.019357 -0.024083
3 486.25 0.974692 0.128219 0.025009 0.003139 0.014035 0.007913 0.092204 0.000594 0.000150 ... -0.120759 -0.079569 0.482182 0.455254 -0.317632 -0.206404 -0.065104 0.022460 0.065246 0.070507
4 487.00 0.972715 0.131813 0.029218 0.000881 0.008492 0.013129 0.078673 0.000717 0.000138 ... -0.277626 -0.092605 0.103249 0.414886 -0.316136 -0.200999 -0.090632 0.077541 0.111759 0.061460

5 rows × 62 columns

The features extracted are similar to the features extracted for gait detection, but the gyroscope has been added to extract additional MFCCs of this sensor. The gyroscope (measuring angular velocity) is relevant to distinguish between arm activities. Also note that the time column no longer starts at 0, since the first timestamps were predicted as non-gait and therefore discarded.

Step 5: Filtering gait

This classification task is similar to gait detection, although it uses a different classification object. The trained classifier is a logistic regression, similarly trained on the dataset of the Parkinson@Home Validation Study. Filtering gait is the process of detecting and removing gait segments containing other arm activities. This is an important process since individuals entertain a wide array of arm activities during gait: having hands in pockets, holding a dog leash, or carrying a plate to the kitchen. We trained a classifier to detect these other arm activities during gait, enabling accurate estimations of the arm swing.

from paradigma.classification import ClassifierPackage
from paradigma.pipelines.gait_pipeline import filter_gait

# Set the path to the classifier package
classifier_package_filename = 'gait_filtering_clf_package.pkl'
full_path_to_classifier_package = files('paradigma') / 'assets' / classifier_package_filename

# Load the classifier package
clf_package = ClassifierPackage.load(full_path_to_classifier_package)

# Detecting no_other_arm_activity returns the probability of no_other_arm_activity for each window, which is concatenated to
# the original dataframe
df_arm_activity['pred_no_other_arm_activity_proba'] = filter_gait(
    df=df_arm_activity,
    clf_package=clf_package
)

n_windows = df_arm_activity.shape[0]
n_predictions_no_other_arm_activity = df_arm_activity.loc[df_arm_activity['pred_no_other_arm_activity_proba']>=clf_package.threshold].shape[0]
perc_predictions_no_other_arm_activity = round(100 * n_predictions_no_other_arm_activity / n_windows, 1)
n_predictions_other_arm_activity = df_arm_activity.loc[df_arm_activity['pred_no_other_arm_activity_proba']<clf_package.threshold].shape[0]
perc_predictions_other_arm_activity = round(100 * n_predictions_other_arm_activity / n_windows, 1)

print(f"Out of {n_windows} windows, {n_predictions_no_other_arm_activity} ({perc_predictions_no_other_arm_activity}%) were predicted as no_other_arm_activity, and {n_predictions_other_arm_activity} ({perc_predictions_other_arm_activity}%) as other_arm_activity.")

# Only the time and predicted probabilities are shown, but the dataframe also contains
# the extracted features
df_arm_activity[['time', 'pred_no_other_arm_activity_proba']].head()

Out of 50 windows, 0 (0.0%) were predicted as no_other_arm_activity, and 50 (100.0%) as other_arm_activity.
time pred_no_other_arm_activity_proba
0 484.00 0.000002
1 484.75 0.000004
2 485.50 0.000020
3 486.25 0.000013
4 487.00 0.000011

Step 6: Arm swing quantification

The next step is to extract arm swing estimates from the predicted gait segments without other arm activities. Arm swing estimates can be calculated for both filtered and unfiltered gait, with the latter being predicted gait including all arm activities. Specifically, the range of motion ('range_of_motion') and peak angular velocity ('peak_velocity') are extracted.

This step creates gait segments based on consecutively predicted gait windows. A new gait segment is created if the gap between consecutive gait predictions exceeds config.max_segment_gap_s. Furthermore, a segment is considered valid if it is of at minimum length config.min_segment_length_s. Note that no gait without other arm activities was detected in the prior step, and therefore we shall modify the classification threshold (clf_package.threshold) to ensure we have some quantifications.

from paradigma.pipelines.gait_pipeline import quantify_arm_swing

temporary_clf_threshold = 0.00001
dfs_to_quantify = ['unfiltered', 'filtered']

print(f"The original classification threshold of {clf_package.threshold} is for this tutorial temporarily set to {temporary_clf_threshold}.\n")

clf_package.threshold = temporary_clf_threshold

quantified_arm_swing, segment_meta = quantify_arm_swing(
    df_timestamps=df_preprocessed,
    df_predictions=df_arm_activity,
    classification_threshold=clf_package.threshold,
    window_length_s=config.window_length_s,
    max_segment_gap_s=config.max_segment_gap_s,
    min_segment_length_s=config.min_segment_length_s,
    fs=config.sampling_frequency,
    dfs_to_quantify=dfs_to_quantify
)

print(f"Gait segments are created of minimum {config.min_segment_length_s} seconds and maximum {config.max_segment_gap_s} seconds gap between segments.\n")
for df_name in dfs_to_quantify:
    print(f"A total of {quantified_arm_swing[df_name]['segment_nr'].nunique()} {df_name} gait segments have been quantified.")

# The arm swing quantification is returned as a dictionary with two dataframes: 'unfiltered' and 'filtered'.
# The 'unfiltered' dataframe contains the quantification of all gait segments, while the 'filtered' dataframe
# contains only the segments that are classified as gait and no_other_arm_activity.
dataset_used = 'filtered'

print(f"\nThe first rows of the {dataset_used} dataset")
quantified_arm_swing[dataset_used].head()

The original classification threshold of 0.4691878519863539 is for this tutorial temporarily set to 1e-05.

Gait segments are created of minimum 1.5 seconds and maximum 1.5 seconds gap between segments.

A total of 5 unfiltered gait segments have been quantified.
A total of 3 filtered gait segments have been quantified.

The first rows of the filtered dataset
segment_nr range_of_motion peak_velocity
0 2 7.511588 66.135460
1 2 19.391911 78.645638
2 2 18.928598 95.738234
3 2 14.612863 57.997025
4 2 14.536349 98.900658 
import pprint

print(f"Furthermore, per gait segment the following metadata is stored:\n")
pprint.pprint(segment_meta)

Furthermore, per gait segment the following metadata is stored:

{'filtered': {2: {'segment_category': 'long', 'time_s': 2.25},
              3: {'segment_category': 'long', 'time_s': 0.75},
              4: {'segment_category': 'moderately_long', 'time_s': 6.75},
              5: {'segment_category': 'moderately_long', 'time_s': 6.0},
              6: {'segment_category': 'moderately_long', 'time_s': 0.75}},
 'unfiltered': {1: {'segment_category': 'long', 'time_s': 12.75},
                2: {'segment_category': 'moderately_long', 'time_s': 9.0},
                3: {'segment_category': 'long', 'time_s': 12.0},
                4: {'segment_category': 'moderately_long', 'time_s': 7.5},
                5: {'segment_category': 'moderately_long', 'time_s': 7.5}}}

The segment categories are defined as follows:

  • short: < 5 seconds

  • moderately_long: 5-10 seconds

  • long: 10-20 seconds

  • very_long: > 20 seconds

As noted before, the segments (and categories) are determined based on predicted gait (unfiltered gait). Therefore, for the arm swing of filtered gait, a segment may be smaller as parts of the segment were predicted to have other arm activities, yet the category remained the same.

Step 7: Aggregation

Finally, the arm swing estimates can be aggregated across segments.

from paradigma.pipelines.gait_pipeline import aggregate_arm_swing_params

dataset_used = 'filtered'

arm_swing_aggregations = aggregate_arm_swing_params(
    df_arm_swing_params=quantified_arm_swing[dataset_used],
    segment_meta=segment_meta[dataset_used],
    aggregates=['median', '95p']
)

pprint.pprint(arm_swing_aggregations, sort_dicts=False)

{'moderately_long': {'time_s': 13.5,
                     'median_range_of_motion': 14.394662839301768,
                     '95p_range_of_motion': 22.658748546739268,
                     'median_peak_velocity': 59.475147961606005,
                     '95p_peak_velocity': 105.31356363194875},
 'long': {'time_s': 3.0,
          'median_range_of_motion': 14.612862824791089,
          '95p_range_of_motion': 19.29924826375146,
          'median_peak_velocity': 78.64563778043262,
          '95p_peak_velocity': 98.26817330011067},
 'all_segment_categories': {'time_s': 16.5,
                            'median_range_of_motion': 14.57460571480515,
                            '95p_range_of_motion': 22.180208543070503,
                            'median_peak_velocity': 65.54184788176072,
                            '95p_peak_velocity': 105.21165557470036}}

The output of the aggregation step contains the aggregated arm swing parameters per gait segment category. Additionally, the total time in seconds time_s is added to inform based on how much data the aggregations were created.