How Measuring Human Movement Helps Study Sleep and Wake Patterns

Over the past four decades, numerous heuristic sleep–wake scoring algorithms have been developed, with more recent approaches incorporating machine learning [1]. In this study, we used one representative algorithm designed for epoch-level activity data and another for high-temporal-resolution acceleration data to assess the impact of actigraphic signal processing on sleep–wake classification, aiming to detect sustained, consolidated sleep periods for the sake of comparability.

The Munich Actimetry Sleep Detection Algorithm (MASDA, also known as Roenneberg algorithm) is a device-independent activity-data-based sleep–wake scoring algorithm that was validated on ActTrust devices produced by Condor Instruments and subsequently applied to activity data collected with MotionWatch 8 devices developed by CamNtech Ltd. [28,29]. Briefly, MASDA estimates the trend of activity data using a centered 24 h moving average and classifies an epoch as sleep if its activity value falls below 15–25% of the trend; the resulting raw sleep–wake classification then undergoes two-stage automatic rescoring, yielding long stretches of consolidated sleep and wake [30]. For this analysis, the open-source pyActigraphy (v1.2.1) implementation [31] of MASDA was used with default parameters, and the threshold was set to 25% as per the recommendation [28].

The van Hees algorithm from GGIR (a research-community-driven R package developed for transparent, vendor-agnostic processing and analysis of actigraphic acceleration data) is specially designed for acceleration-data-based sleep–wake scoring. Briefly, it relies on two parallel steps whose combination yields the final classification [32]: the detection of sustained inactivity bouts (SIBs [12]) and sleep period time (SPT [33]) window. SIBs are identified when the change in wrist angle remains within 5° for at least 5 min. SPT-window detection also relies on wrist angle but includes an automated cleaning step to ensure there is a single contiguous SPT per day (noon-to-noon). Sleep is defined based on the overlap between the contiguous SPT and multiple shorter SIBs. In this study, GGIR (v3.2.6) was run with default settings, except that the SPT threshold was set to 0.4 following visual inspection; which is permitted for methodological research according to the creators of the package [34]. Note that since the acceleration data we have examined (sampled at 10 Hz—see Section 2.3.1) included real-time clock-based, millisecond precision time stamps, it was necessary to reduce the resolution of the time stamps to obtain time-equidistant samples spaced exactly 100 milliseconds apart, as required by the algorithm.

The characterization of biological rest–activity rhythms has a long history. In this task, features are extracted from activity data to quantify distinct aspects of the rhythm. Here, the effects of actigraphic signal processing were assessed for the 5 most-used indicators in nonparametric circadian rhythm analysis [9].

L5 and M10 capture the least-active 5 h and most-active 10 h windows within a day with onset times (L5onset and M10onset) and mean activity levels (L5val and M10val) [9]. Thus, while L5 indexes nocturnal activity, M10 indexes diurnal activity: a lower L5val implies more restful sleep, and a higher M10val indicates greater daytime activity. Relative Amplitude (RA) is used to compute the normalized difference between M10val and L5val (see Equation (1)). RA ranges from 0 to 1; larger values indicate greater separation between nighttime and daytime activity and are generally taken to reflect a more robust rest–activity rhythm. We determined L5 and M10 in the traditional way [13,14,15,17,18,19,20,21,22,23,24,35,36]: first, a 24 h mean profile was obtained by averaging the minute-level activity data across days; then, L5 and M10 were derived by sliding a window in one-minute steps across the 24 h mean profile in a circular manner and computing the mean activity within each window.(1)RA=M10val−L5valM10val+L5val

Interdaily Stability (IS) quantifies the day-to-day predictability (i.e., similarity) of a subject's activity pattern and is computed using Equation (2), where xi denotes the mean activity in each hour in the multi-day activity data, xm is the grand mean across the entire time series, and xh is the mean activity of the h-th hour in the 24 h mean profile [9]. IS ranges from 0 to 1; larger values denote more consistent day-to-day activity patterns.(2)IS=124∑h=124xh−xm21N∑i=1Nxi−xm2

Intradaily Variability (IV) measures within-day fragmentation of the rest–activity rhythm [9] and is computed as shown in Equation (3), where xi denotes the mean activity in each hour in the multi-day activity data, and xm is the grand mean across the entire time series. IV ranges from 0 to 2; higher values indicate more frequent rest–activity transitions (e.g., daytime naps and nighttime awakenings), i.e., greater fragmentation.(3)IV=1N−1∑i=2Nxi−xi−121N∑i=1Nxi−xm2

In most cases, actigraphs using MEMS (micro-electromechanical system) accelerometers measure the x-, y-, and z-axis components of limb acceleration (UFX, UFY, and UFZ, collectively termed UFXYZ). During preprocessing, the primary aims are to combine these axial signals into a single series and remove the effect of Earth's gravitational acceleration (i.e., the gravitational component, g). The former can be achieved by computing the Euclidean norm of the per-axis accelerations to obtain the unfiltered magnitude data (UFM). The effect of g can be removed from this magnitude data by subtracting 1 g and either taking the absolute value (UFNM) or truncating negative differences to zero (ENMO). Alternatively, the DC (gravitational) component can be attenuated using digital filters (e.g., a 3rd-order Butterworth band-pass filter with 0.25 and 2.5 Hz cutoffs), applied either to the UFM data (FMpost) or to the per-axis signals (FX, FY, and FZ, collectively denoted as FXYZ) before computing the Euclidean norm (FMpre).

Following preprocessing, the acceleration data is partitioned into fixed-length epochs (in our case, 1 min), and an activity metric is applied to yield one activity value per epoch. The most classical metrics either simply integrate the acceleration (Proportional Integration Method—PIM) or use fixed thresholds to determine the crossings (Zero Crossing Method—ZCM) or the time the signal remains above such thresholds (Time Above Threshold—TAT). Contrary to the name, ZCM, the threshold level of these intersection-based activity metrics should be set above the noise level rather than at zero, though its exact value has not been disclosed by the manufacturers. Previously [2], we determined its optimal level, and we have accordingly set it equal to the standard deviation (note that for UFM, 1 g should be added to account for the presence of the gravitational component) of all the preprocessed acceleration data in this analysis. Other metrics are based on calculating the standard deviation (Mean Amplitude Deviation—MAD) or the variance (Activity Index—AI) of the acceleration data in a given epoch. The metric High-Pass Filtered Euclidean Norm (HFEN) simply averages the acceleration data, but it requires a specific preprocessing scheme: the per-axis signals must be filtered with a high-pass—not a band-pass—filter (4th-order Butterworth high-pass filter with 0.2 Hz cutoff) before magnitude calculation (denoted as HFMpre rather than FMpre).

Preprocessing techniques and activity metrics can be properly combined in 35 different ways, each of which is a different activity determination method that we denote here in accordance with our past work: the activity metric is treated as a function whose argument is the preprocessing technique. Thus, for example, PIM(UFNM) denotes such an activity signal that has been obtained by minutely integrating the normalized magnitude of acceleration. For more details of these procedures, see our prior work [2].

Given the widespread use of devices from Ametris LLC [13,14,15] and CamNtech Ltd. [16,17] in NPCRA studies, the corresponding activity determination methods have also been included in the present analysis. We have used the Python implementation of Ametris's algorithm (agcounts, v0.6.9), which the company published 3 years ago to enhance transparency [37] to address the scientific community's demand for accessible and reproducible methods [3]. Their solution is comparatively complex and does not align one-to-one with the methods in Section 2.2.1, as epoch-level activity values (the so-called Activity Count; hereafter referred to as AC) are computed via multiple stages (resampling, band-pass filtering, rescaling, thresholding, and summing—see [3] for details). CamNtech disclosed the activity determination procedure of their MotionWatch 8 devices in their user manual [38]: each axis is band-pass-filtered prior to magnitude computation (FMpre-like preprocessing, with 3 and 11 Hz cutoffs); per-second peak acceleration is then determined, and peaks exceeding 0.1 g are summed within each epoch to determine the activity value (hereafter referred to as MW). As the acceleration data we examined were recorded at a sampling rate of 10 Hz, we have approximated this band-pass filter using the one specified in our generalized framework.

Beyond these approaches, the PIM(ENMO) activity-determination method in our generalized scheme corresponds directly to the widely used GGIR package (v3.2.6), which applies ENMO by default and derives epoch-level values through simple averaging, although it uses 5 s epochs by default [32]. However, since correlation-based examinations require consistent epoch length across activity determination methods and in circadian rhythm analysis—where the focus is on the overall multi-day activity pattern rather than short bursts—a 1 min epoch length is preponderant and recommended [13,16,17,18,20,23,39] and the use of ENMO with 1 min epochs is not unprecedented [40], we employed it as well. Furthermore, the ZCM(FMpre), TAT(FMpre), and PIM(FMpre) methods correspond to those implemented in Condor Instruments' ActTrust devices [18,19,20,21], with the sole difference being that ActTrust employs a 0.5–2.7 Hz band-pass filter rather than the 0.25–2.5 Hz band-pass used here [41]. Note that the ZCM, TAT, and PIM metrics are also supported by long-standing and widely used devices produced by Ambulatory Monitoring Inc. [21,22,23,24,25]; however, because these devices use piezoelectric rather than MEMS sensors and only fragmented information is available regarding their algorithmic details [42], the correspondence is necessarily indirect.

The publicly available actigraphic dataset we collected and examined in several of our previous studies [2,6,27], as well as in this work, contains raw triaxial acceleration data obtained from 42 healthy, free-living subjects [43]. The 10-day-long recordings were acquired on each subject's non-dominant wrist at a sampling rate of 10 Hz with a ±8 g dynamic range and 16 mg resolution using a special-purpose MEMS-accelerometer-based actigraphic device we developed [2]. To reduce deterministic measurement errors (scale, offset, and orthogonality), the acceleration signals were subjected to post-calibration using a publicly available algorithm [44].

The impact of actigraphic signal processing on sleep–wake scoring and each NPCRA indicator was examined through separate similarity matrices, in which each cellquantifies the discrepancy across the 42 subjects that would arise in the given indicator value or sleep–wake score if the subjects had been effectively measured using an alternative actigraphic device—and thus their locomotor activity would have been quantified through a different, manufacturer-dependent method (e.g., MW vs. AC). M i j ,

For NPCRA indicators with scalar outcomes (L5val, M10val, RA, IS, and IV), the Symmetric Mean Absolute Percentage Error (SMAPE) was used as a similarity measure to facilitate comparability with the correlation-based similarity matrix (Figure 2). The method of computing a single similarity-matrix cell Mi,j is defined in Equation (4), where N is the number of participants; for participant s, Vsi and Vsj are the values of a given indicator (e.g., L5val), derived from actigraphic signals produced by signal-processing pipelines i and j. SMAPE ranges from 0% to 200%, with lower values indicating smaller discrepancies in the given indicator that are attributable to the choice of signal-processing pipeline (e.g., device selection). Note that L5val and M10val are not ratio-type measures (unlike RA, IS, and IV); therefore, the differences in the value ranges of differently computed activity data might introduce order-of-magnitude discrepancies in these indicators [13,45]. Thus, for fair comparison of L5val and M10val, we have normalized the actigraphic data by their respective means to reduce such distortions.(4)Mi,j=1N∑s=1NVsi−VsjVsi+Vsj/2·100

For NPCRA indicators for which the outcome of interest is the onset (L5onset and M10onset), the similarity matrices were constructed from the overlap (ranging from 0% to 100%) of the designated windows, as defined in Equation (5). Here,denotes the number of participants; for participant,andare the designated windows for the given indicator, derived from actigraphic signals produced by signal-processing pipelinesand, andis the segment length (5 h for L5onset and 10 h for M10onset). N s W s i W s j i j L (5) M i j , = · 1 N ∑ s = 1 N W s i W s j ∩ L 100

For sleep–wake scoring algorithms, two aspects have been evaluated: the timing of sleep and Total Sleep Time (TST). In the former case, analogous to onset-based NPCRA indicators, overlap was assessed; however, because the lengths of sleep segments were not uniform, overlap was quantified via the intersection over union (IoU), i.e., the Jaccard index, as in Equation (6). In the formula,denotes the number of participants; for participant,andare the designated sleep intervals from actigraphic signals, produced by signal-processing pipelinesand. N s S s i S s j i j (6) M i j , = ⋅ 1 N ∑ s = 1 N S s i S s j ∩ S s i S s j ∪ 100

Additionally, we employed dendrograms to reveal the structure and hierarchical organization of the similarity matrices. Prior to this, matrices were converted to dissimilarities and rescaled between 0 and 1, with lower values indicating smaller differences. Dendrograms were computed using SciPy (v1.7.3)—with the default Euclidean metric and complete linkage—which is less sensitive to outliers and noise [46].

In this work, we also connect our prior findings about the correlation-coefficient-based relationships between the different activity determination methods (see [2] for details) to present results. Figure 2 displays the correlation-based similarity matrix, which not only aids in the interpretation of our current results (see Section 3) in light of our previous findings but also serves to illustrate the rationale behind our matrix-based comparison. In the correlation matrix, the closer the r value is to 1, the stronger the linear relationship between activity signals derived from the same acceleration data (see Section 2.3.1) but computed using different methods (Section 2.2.1). These relationships are mainly driven by how the acceleration signal was preprocessed, an observation we reinforced by color coding the labels to indicate different preprocessing families: purple (raw per-axis acceleration), light and dark blue (per-axis filtered acceleration), red (magnitude-filtered acceleration), and green (unfiltered magnitude of acceleration).

Figure 3 illustrates how the mean values of L5, M10, and RA across 42 subjects depend on which signal-processing pipeline in our generalized framework was used to determine the acceleration or activity signal. For each NPCRA indicator, 43 mean values are represented with colored dots: 35 for those activity signals that previously formed the rows and columns in the correlation matrix of Figure 2, and 8 for acceleration signals that have undergone various preprocessing procedures. The coloring of the dots in Figure 3 indicates three categories: blue denotes values computed from epoch-level activity data derived with threshold-based metrics (TAT and ZCM), red denotes values from activity data derived with non-threshold metrics, and black denotes values computed directly from higher-temporal-resolution acceleration data preprocessed in different ways. As shown, applying the L5 indicator to threshold-based activity signals (blue) yields systematically lower values than applying it to activity signals derived with any non-threshold metric (red). By contrast, for M10, the pattern reverses: threshold-based metrics yield higher values. Notably, the indicator values computed from acceleration signals (black) closely resemble those computed from activity signals obtained with non-threshold metrics (red). Since RA is derived from L5 and M10 values, it also forms two distinct clusters, similar to those observed for the two aforementioned indicators.

In Figure 4, the counterpart to the previous figure is shown for the mean values of the IV and IS indicators. In this case, however, the color coding of the mean values—contrary to Figure 3—is based on the preprocessing families, since they group according to the preprocessing techniques applied to the acceleration signal and not based on whether the activity metric was threshold-based or not. These clusters that clearly resemble the structure of the correlation-based similarity matrix (see Figure 2) appear in the same order for IS and IV—albeit with opposite directionality.

Figure 3 and Figure 4 clearly demonstrate that the values of all the NPCRA indicators considered can differ substantially depending on how the actigraphic signals are derived within our generalized framework. For example, IV was reported to be around 0.67 for the controls and 0.811 for subjects suffering from bipolar disorder [47], while IS was measured at approximately 0.446 for the controls and 0.526 for schizophrenic subjects [48]; but as shown in Figure 4, even greater discrepancies can occur within the same group of healthy subjects due to the lack of standardization in actigraphy. Moreover, different NPCRA indicators are more strongly influenced by different aspects of the actigraphic methodology: for the values of L5, M10, and RA, the activity metric exerts the greater influence, whereas for IS and IV, the preprocessing technique applied to the acceleration signal has a larger effect. The impact of these aspects has been quantified using similarity matrices and dendrograms, as defined in Section 2.3.2.

Figure 5 presents the SMAPE-based similarity matrices and their corresponding dendrograms for L5 values; the labels follow the same categorical color coding as for Figure 3. It makes what was inferred from the means in Figure 3 even clearer: the L5 values differ substantially (the differences are generally greater than 80%) depending on whether L5 was applied to activity signals derived using threshold-based metrics (blue labels), to activity signals derived using other metrics (red labels), or directly to acceleration signals (black labels). The corresponding rationale will be addressed in Section 4. For the values of M10 and RA (corresponding matrices are provided in Figure S1 and Figure S2 of the Supplementary Material, respectively), the same clustering pattern was observed. However, while the SMAPE-based difference between the two clusters was generally greater than 80% for L5, it is below 10% for M10 and around 20% for RA. Since RA is directly derived from the value of M10 and L5, it means that the observed discrepancies for RA are mainly driven by the differences in the value of L5.

As shown in Figure 6, for the IS indicator (the corresponding matrix for IV can be found in Figure S3 of the Supplementary Material), the differently determined acceleration and activity signal types form hierarchical clusters along the choice of acceleration-signal-preprocessing technique—mirroring the correlation-based similarity matrix of Figure 2 and the separation of mean values observed for IS and IV in Figure 4—but the differences between the clusters are generally smaller compared to what we observed for the value of L5 in Figure 5. An important implication follows: for a given preprocessing technique applied to the acceleration signal, IS and IV computed directly based on acceleration data generally agree closely with those computed from activity signals derived from the same acceleration data using any activity metrics. In contrast, the values of L5, M10, and RA can differ substantially depending on whether they are computed from acceleration signals or activity signals if the latter were derived using threshold-based metrics.

Figure 7 illustrates the extent of overlap between the 5 least active and 10 most active 10 h periods identified by L5 and M10. As shown, both matrices exhibit similar structures; however, the onset of L5 appears less affected by the choice of actigraphic methodology, as the designated periods generally overlap to a much greater extent compared to M10. Additionally, the patterns of these matrices resemble those observed for IS and IV (see Figure 6) or in the correlation analysis (see Figure 2), suggesting that the onset of L5 and M10 is primarily influenced by how the acceleration data were preprocessed rather than the choice of activity metric used for epoch-level data compression, a finding that is the opposite of what we observed for the values of L5 and M10 (see Figure 5). This further reinforces the necessity of considering the entire activity-determination process when comparing results across studies, as L5 and M10 are largely influenced by both main steps of the actigraphic signal-processing pipeline.

Figure 8 presents the spread of L5, M10, and RA indicator values (each on separate panels, (a) to (c)) across the 42 subjects for specific devices. Since AC and MW are also threshold-based activity metrics like ZCM and TAT, the values form clusters within each indicator depending on whether the activity data were calculated using threshold-based activity metric. This clustering is further substantiated by the dendrogram in panel (d), which is based on the RA. These findings are in line with the observations we made regarding the generalized activity determination methods (see Figure 3). For example, the device specificities of RA are mainly due to the differences in the value of L5, as the interquartile ranges of different devices overlap for the value of M10, but not for the value of L5 and RA. It is crucial to highlight that ZCM, TAT, and PIM are the three classical activity metrics supported by both AMI and Condor Instruments devices. However, since selecting between them is up to scientists when preparing for measurement, only this choice can have a substantial impact on the NPCRA outcomes as shown. For example, the RA value was, on average, close to 0.70 for PIM and 0.86 for ZCM (i.e., an average difference of 0.16 for an indicator that ranges between 0 and 1), and we found that their subject-level difference was statistically significant using a paired t-test (p < 0.001). Since RA values were derived for the same motion, these discrepancies are purely a consequence of how the actigraph was set up (or which device was used), and, remarkably, they are greater for 38 out of the 42 subjects compared to what was observed between pediatric subjects and adults by others using the same device [13].

Figure 9 shows the counterpart of the previous figure for the IS and IV indicator values in panels (a) and (b), while panel (c) depicts the dendrogram for IS. The dendrogram clearly shows that the IS derived from activity data based on the activity determination method from GGIR (i.e., PIM(ENMO)) substantially differs from that derived via the rest of the methods. Since PIM(ENMO) is the only activity determination method in this narrower set where digital filtering was not involved in the preprocessing of acceleration data, this behavior aligns with the observation we made regarding the generalized activity determination methods (see Figure 4): the differences device selection can introduce in the IS and IV are primarily driven by how the acceleration data have been preprocessed rather than the activity metric. For both IV and IS, the greatest discrepancies occurred between PIM(ENMO) from the GGIR package and ZCM(FMpre) supported by the AMI and Condor Instruments devices (for example, on average, 0.99 vs. 0.75 for IV, respectively).

The limited set of specific devices' activity determination methods did not allow us to comprehensively describe the similarity structure based on NPCRA indicator onset times as we did for the generalized framework, for which we compared 43 methods (Section 3.1.2) instead of only 6; therefore, we instead examined another timing-related metric in this case: sleep–wake scoring. We assessed the impact of device choice on sleep–wake scoring by classifying sleep and wakefulness from raw acceleration data using the van Hees algorithm and by applying the MASDA algorithm to six different types of activity signals computed from the same acceleration data to emulate specific devices commonly used in the related literature. In Figure 10, one can see, for an exemplary subject, that even inter-device selection between ZCM and PIM activity metrics in AMI and ActTrust devices can result in several hours of discrepancies in sleep–wake scoring.

To systematically assess such discrepancies, we first examined the total amount of sleep accumulated over 10 days, on average, among 42 participants, depending on whether sleep was identified from raw acceleration data using the van Hees algorithm or from activity signals using the MASDA algorithm; this comparison is illustrated in Figure 11a. For MASDA, highly similar total sleep durations were obtained across activity signals derived from threshold-based metrics (TAT(FMpre), ZCM(FMpre), AC, MW) ranging from 79.3 to 80.3 h. In contrast, activity signals based on integration yielded shorter sleep durations: PIM(FMpre) produced roughly 4 h less sleep on average than the threshold-based signals, while PIM(ENMO) showed an exceptionally low TST of only 50.6 h. As shown earlier (see Figure 9), the PIM(ENMO) signal also behaved anomalously for the IV and IS indicators, likely due to the absence of digital filtering during acceleration preprocessing. Sleep–wake detection from the raw acceleration signal yielded a mean total sleep duration comparable to that estimated by MASDA for the PIM(FMpre) signal.

Beyond the similarity in TST estimated via the different methods, it is also crucial to assess how consistently they identified the onset of sleep. To this end, we examined the overlap of sleep segments identified through the different approaches using the IoU metric, as shown in Figure 11b. Given that the MASDA algorithm produced different TSTs when applied to PIM(ENMO), it is evident that it also identified sleep at substantially different times compared to other methods, with overlaps below 62% in all cases. Comparing acceleration- and activity-based sleep–wake scoring highlights the importance of overlap analyses, since TST alone fails to capture the temporal agreement between methods. While sleep–wake scoring based on acceleration data (the van Hees algorithm applied to UFXYZ) yielded an almost identical TST to the MASDA algorithm applied to activity data, the IoU-based overlap of corresponding sleep segments was consistently below 78% in all cases, underscoring the methodological difference. For instance, although TST differed by only 0.67% between UFXYZ and PIM(FMpre), their sleep segment overlap reached only 74.8%; this disagreement is also illustrated in Figure 12 for one representative participant.

The parameters of the MASDA and van Hees sleep–wake scoring algorithms could be fine-tuned for each signal type to reduce their discrepancies, especially for MASDA when applied to PIM(FMpre) activity data. However, even tuning on specific datasets is not a common practice in typical use cases in the related literature, which we also aimed to replicate in order to examine the extent to which manufacturer-specific actigraphy signal processing affects sleep–wake scoring. Even the fact that we did not simply retain the default parameters, but instead selected MASDA settings within the recommended range that yielded acceptable results across the differently derived activity signals, already represents a more careful approach. Note that we used sleep–wake scoring solely to illustrate the discrepancies introduced by the actigraphic methodology of specific devices from a timing-oriented perspective (similar to the onset of L5 and M10 in) without focusing on sleep detection performance; therefore, we did not aim to evaluate the specificity or selectivity of the different sleep–wake scoring algorithms in the absence of a gold standard. Section 3.1.2