Source code for

Calculate physical activity (PA) levels with conversion to activity counts.

Functions for converting raw sensor data to physical activity (PA) or
moderate-to-vigorous physical activity (MVPA) counts, similar to those given
by dedicated accelerometers such as Actigraph devices.

For a uniaxial accelerometer, the signal should first be passed into
:func:``, then the counts should categorized
using :func:``.

For a triaxial accelerometer, an additional step is required. Each axis
should first be passed into :func:`` separately,
then the 3 count vectors should be passed into
:func:`sensormotion.signal.vector_magnitude` to calculate vector
magnitude (VM) of the counts. Finally, the single VM vector should be
categorized using :func:``.

import matplotlib.pyplot as plt
import numpy as np
import sensormotion.signal

from scipy.integrate import simps, trapz

[docs]def convert_counts( x, time, time_scale="ms", epoch=60, rectify="full", integrate="simpson", plot=False, fig_size=(10, 5), ): """ Convert acceleration to physical activity (PA) counts. Given an acceleration signal over a **single** axis, integrate the signal for each time window (epoch). The area under the curve for each epoch is the physical activity count for that time period. Linearly interpolated values are used if exact multiples of `epoch` are not found in the time signal. Parameters ---------- x : ndarray Acceleration signal to be converted to PA counts. time : ndarray Time signal associated with `x`. time_scale : {'ms', 's'}, optional The unit that `time` is measured in. Either seconds (s) or milliseconds (ms). epoch : int, optional The duration of each time window in seconds. Counts will be calculated over this period. PA counts are usually measured over 60 second epochs. :func:`` also requires 60 second epochs, however, if you're using your own cut point set and just want raw count values feel free to use any sized epoch. rectify : {'full', 'half'}, optional Type of rectifier to use on the input acceleration signal. This is to ensure that PA counts take into consideration negative acceleration values. Full-wave rectification turns all negative values into positive ones. Half-wave rectification sets all negative values to zero. integrate : {'simpson', 'trapezoid'}, optional Integration method to use for each epoch. plot : bool, optional Toggle to display a plot of PA counts over time. fig_size : tuple, optional If plotting, set the width and height of the resulting figure. Returns ------- counts : ndarray PA count values for each epoch. """ assert len(x) == len(time), "signal and time must be the same length" assert np.all(np.diff(time) > 0), "time signal is not fully ascending" assert ( integrate == "simpson" or integrate == "trapezoid" ), "integrate method must either be simpson or trapezoid" x = np.asarray(x) time = np.asarray(time) # convert time to seconds if time_scale == "ms": time = time / 1000 elif time_scale == "s": time = time # calculate time diff time = time - time[0] assert max(time) > epoch, "length of signal time shorter than epoch size" # rectify signal x = sensormotion.signal.rectify_signal(x, rectify) # interpolate missing times values to get exact epochs boundary_count = int(max(time) / epoch) + 1 boundary_times = [i * epoch for i in range(boundary_count)] missing_times = np.setdiff1d(boundary_times, time) # epoch times to interp x = np.append(x, np.interp(missing_times, time, x)) # interpolate x values time = np.append(time, missing_times) # sort new time and signal arrays together sort_idx = time.argsort() time = time[sort_idx] x = x[sort_idx] # get index of each epoch/boundary value for slicing boundary_idx = np.where(np.isin(time, boundary_times))[0] # integrate each epoch using Simpson's rule counts = np.ones(len(boundary_idx) - 1) # preallocate array for i in range(len(counts)): lower = boundary_idx[i] upper = boundary_idx[i + 1] + 1 # upper bound should be inclusive cur_signal = x[lower:upper] cur_time = time[lower:upper] if integrate == "simpson": counts[i] = simps(cur_signal, cur_time) elif integrate == "trapezoid": counts[i] = trapz(cur_signal, cur_time) # plot counts if plot: f, ax = plt.subplots(1, 1, figsize=fig_size)[1:], counts, width=epoch - 5) plt.xticks( boundary_times[1:], [ "{} - {}".format(boundary_times[i], boundary_times[i + 1]) for i, x in enumerate(boundary_times[1:]) ], ) for tick in ax.get_xticklabels(): tick.set_rotation(45) plt.suptitle("Physical activity counts", size=16) plt.xlabel("Time window (seconds)") plt.ylabel("PA count") return counts
[docs]def cut_points(x, set_name, n_axis, plot=False, fig_size=(10, 5)): """ Categorize physical activity (PA) counts into intensity categories. Use a pre-specified cut-point set to determine the intensity level of each epoch of PA. Cut-point sets are from published research articles, and a summary can be found in the Actigraph FAQ: **Important**: The cut-point sets used here assume each count epoch is 60 seconds long. If you're using counts from :func:`` make sure you had set the value of epoch to 60. Don't use this function if you use different length epochs. Parameters ---------- x : ndarray Vector of counts calculated by :func:`` or vector magnitude from :func:`sensormotion.signal.vector_magnitude`. This can either be from a single axis or a vector magnitude, but set the `n_axis` parameter accordingly. set_name : {'butte_preschoolers', 'freedson_adult', 'freedson_children', 'keadle_women'} The name of the cut-point set to use. These are listed in the Actigraph FAQ, and the corresponding research article is shown below: `butte_preschoolers`: Butte et al. (2013) - Prediction of Energy Expenditure and Physical Activity in Preschoolers `freedson_adult` (uniaxial): Freedson et al. (1998) - Calibration of the Computer Science and Applications, Inc. accelerometer `freedson_adult` (triaxial): Freedson et al. (2011) - Validation and Comparison of ActiGraph Activity Monitors `freedson_children`: Freedson et al. (2005) - Calibration of accelerometer output for children `keadle_women`: Keadle et al. (2014) - Impact of accelerometer data processing decisions on the sample size, wear time and physical activity level of a large cohort study n_axis : {1, 3} Number of axes over which acceleration was recorded (1 = uniaxial, 3 = triaxial accelerometer). This is used to determine the cut-point set values to use as some sets are for counts over a single axis, whereas others are thresholds for vector magnitude (calculated from 3 axes). plot : bool, optional Toggle to display a plot of PA counts with category thresholds marked. fig_size : tuple, optional If plotting, set the width and height of the resulting figure. Returns ------- category : list List of intensity categories for each epoch of PA, as determined by the cut-point set thresholds. time_spent : ndarray Amount of time (counts) spent at each PA intensity level. """ # new cut-point sets should be added to this dictionary sets = { "butte_preschoolers": { 1: { "sedentary": [-np.inf, 239], "light": [240, 2119], "moderate": [2120, 4449], "vigorous": [4450, np.inf], }, 3: { "sedentary": [-np.inf, 819], "light": [820, 3907], "moderate": [3908, 6111], "vigorous": [6112, np.inf], }, }, "freedson_adult": { 1: { "sedentary": [-np.inf, 99], "light": [100, 1951], "moderate": [1952, 5724], "vigorous": [5725, 9498], "very vigorous": [9499, np.inf], }, 3: { "light": [-np.inf, 2690], "moderate": [2691, 6166], "vigorous": [6167, 9642], "very vigorous": [9643, np.inf], }, }, "freedson_children": { 1: { "sedentary": [-np.inf, 149], "light": [150, 499], "moderate": [500, 3999], "vigorous": [4000, 7599], "very vigorous": [7600, np.inf], } }, "keadle_women": { 1: { "sedentary": [-np.inf, 99], "light": [100, 1951], "moderate": [1952, np.inf], }, 3: { "sedentary": [-np.inf, 199], "light": [200, 2689], "moderate": [2690, np.inf], }, }, } try: cur_set = sets[set_name][n_axis] print("Cut-point set: {} (axis count: {})...".format(set_name, n_axis)) for i in cur_set: print("{}: {} to {}".format(i, cur_set[i][0], cur_set[i][1])) except KeyError: print( "Error: cut-point set not found. Make sure the set name and/or " "number of axes are correct" ) raise # categorize counts category = [] for count in x: for intensity in cur_set: if cur_set[intensity][0] <= count <= cur_set[intensity][1]: category.append(intensity) break # count time spent category_unique, category_count = np.unique(category, return_counts=True) time_spent = np.asarray((category_unique, category_count)) # plot counts with intensity categories if plot: boundaries = [(item, cur_set[item][0]) for item in cur_set] boundaries.sort(key=lambda x: x[1]) f, ax = plt.subplots(1, 1, figsize=fig_size), len(x) + 1), x) for line in boundaries[1:]: if line[1] < max(x): plt.axhline(line[1], linewidth=1, linestyle="--", color="k") t = plt.text(0.4, line[1], line[0], backgroundcolor="w") t.set_bbox(dict(facecolor="w", edgecolor="k")) plt.xticks(range(1, len(x) + 1)) plt.suptitle("Physical activity counts and intensity", size=16) plt.xlabel("Epoch (length: 60 seconds)") plt.ylabel("PA count") return category, time_spent