On the Unification of Common Actigraphic Data Scoring Algorithms

Abstract

Actigraphy is a well-known, inexpensive method to investigate human movement patterns. Sleep and circadian rhythm studies are among the most popular applications of actigraphy. In this study, we investigate seven common sleep-wake scoring algorithms designed for actigraphic data, namely Cole-Kripke algorithm, two versions of Sadeh algorithm, Sazonov algorithm, Webster algorithm, UCSD algorithm and Scripps Clinic algorithm. We propose a unified mathematical framework describing five of them. One of the observed novelties is that five of these algorithms are in fact equivalent to low-pass FIR filters with very similar characteristics. We also provide explanations about the role of some factors defining these algorithms, as none were given by their Authors who followed empirical procedures. Proposed framework provides a robust mathematical description of discussed algorithms, which for the first time allows one to fully understand their operation and basics.

1. Introduction

Actigraphy is a well-recognized method for investigation of human activity, with broad applications spanning consumer technology [1], supplementary use in diagnostics, as well as in research (e.g., sleep and circadian studies [2] or ergonomics [3]). Actigraph is typically a wrist-worn device, with a built-in accelerometer, battery and memory storage. Depending on an actigraph’s architecure, acceleration data can be digitized and stored in different modes, including zero-crossing mode (ZCM), time above threshold (TAT), digital integration (DI) [2,4] or raw acceleration [5].

As an actigraph can collect a longitudinal information about one’s movement profile, the application of actigraphy as a complementary or supplementary tool for diagnosis has been investigated for various sleep or circadian rhythm-related disorders. A recent review and meta-analysis by the American Academy of Sleep Medicine (AASM) Task Force [6] indicates the possible usefulness of actigraphy in, e.g., insomnia, circadian dysrhythmia or insufficient sleep syndrome.

For sleep-related actigraphic studies, one of the main objectives of the data analysis process is to distinguish between stages of sleep and wake during bedtime, with a general assumption of different movement intensity in each. Over the decades of research in the field of computational approaches for sleep-wake stages discrimination, various algorithms have been developed and investigated with respect to their accuracy, sensitivity or specificity in distinguishing the stages and bidirectional transitions between them. Part of such approaches relies mainly on an epoch-by-epoch movement intensity analysis, however, currently new methods based on machine learning techniques are being recognized as well [7,8,9].

In the work of Tilmanne et al. [10], several actigraphy sleep scoring algorithms were mentioned: Webster algorithm [11] as the earliest automatic approach; a group of Sadeh algorithms [12,13,14] and Cole-Kripke algorithm, which are the most recognizable; and Sazonov algorithm [15]. Another two algorithms, which seem to be less popular, yet can be found in research works, are: University of California San Diego (UCSD) algorithm [16] (e.g., in [17]) and Scripps Clinic [18] algorithms (e.g., in [19]). In a recent work of Palotti et al. [19], Authors compared the performance of some of the abovementioned sleep-wake scoring algorithms (including algorithms of: Webster, Cole-Kripke, Sadeh [13], Sazonov and Scripps Clinic) with more recent, machine learning-oriented approaches.

Instead of performance-oriented comparison, in our previous work [20] (in press) we performed a preliminary investigation of common mathematical features of two algorithms for actigraphic sleep-wake scoring—Sadeh [12] and Cole-Kripke [21]. In this article, the initial research is extended by four more algorithms: Scripps Clinic [18], Sazonov [15], UCSD [16] and Webster [11]. Additionally, we analyse their properties within the unified framework to show similarities and to unveil their common underlying major mode of operation—low-pass filtering.

2. Materials and Methods

2.1. Unified Framework Proposal

In this section, we present a unified view, linking together popular actigraphy algorithms in a mathematically coherent framework, composed of three steps:

1.

Collapsing data into epochs.

2.

Linear convolution with empirically chosen coefficients:

$$Y\left[n\right]=\sum _{i=-N}^{M}W\left(i\right)X[n-i]$$

(1)

where $W\left(i\right)$ is a vector of coefficients, X is the input signal, and $N,M\in \mathbb{N}$.

3.

Rescoring (optional), e.g. as in [21]: “After at least 4 min scored as wake, the next 1 min scored as sleep is rescored wake […]”.

Novelty of the proposed approach consists of:

1.

Properly treating the first step as downsampling, which is a well-known signal processing procedure. This observation allows, e.g., to identify aliasing introduced by the procedures applied thus far, and design correct procedures.

2.

Observing that the coefficients, used in different algorithms for the convolution step—$W\left(i\right)$ in Equation (1)—in all cases actually result in a low-pass finite impulse response (FIR) filtering with very similar cutoff frequencies. This observation allows for an informed design of these filters for new algorithms using signal processing knowledge and tools, and an efficient analysis of the current approaches.

2.2. Defining Popular Actigraphy Algorithms under Unified Framework

In the following sections, we show how the discussed algorithms can be split into the three abovementioned steps, and how their defining equations are really a form of convolution of the signal and some empirically fitted kernel—effectively an FIR filter.

2.2.1. Cole-Kripke Algorithm Family

Cole-Kripke [21], Webster [11], UCSD [16] and Scripps Clinic [18] algorithms start with signal downsampling, and then apply very similar equations on the resampled signal. This similarity allows to name these algorithms as Cole-Kripke family (descriptions and original equations are available in Appendix A.1 and Appendix A.2).

Cole-Kripke family equations are of identical form, differing in the number of elements in the sum and, obviously, values of the coefficients. In the following, we show that the similarities exceed significantly the mere mathematical form of the equations.

Step (2.) of the Cole-Kripke algorithm family can be defined as:

$$D\left[n\right]=\sum _{i=-N}^{M}W\left(i\right)X[n-i]$$

(2)

where $D\left[n\right]$ is the output value for epoch n, $W\left(i\right)$ is a vector of coefficients, X represents the input signal and $N,M\in \mathbb{N}$. In the following we show that algorithms from the Cole-Kripke family are equivalent to FIR filters.

A general FIR filter can be defined as:

$$y\left[n\right]=\sum _{i=0}^{K}b\left(i\right)x[n-i]$$

(3)

where y is the output vector, x is the input vector, b is the vector of coefficients and K is the order of the filter. After changing indices in the sum in Equation (3), so that $K-N=M$, we obtain the following expression:

$$y\left[n\right]=\sum _{i=-N}^{M}b(i+N)x[n-i-N]$$

(4)

After moving the whole output signal by N samples and changing the vector of coefficients so that $W\left(i\right)\equiv b(i+N)$, we obtain the following equation, which is identical to step (2.) of the Cole-Kripke algorithm family (Equation (2)):

$$D\left[n\right]=y[n+N]=\sum _{i=-N}^{M}W\left(i\right)x[n-i]$$

(5)

In other words, since moving the signal by few samples is irrelevant for further analysis, the actual action of step (2.) in the Cole-Kripke algorithm is equivalent to filtering the input signal with an FIR filter of order $N+M$. Rescoring rules (step (3.)) of a form “after at least X minutes scored as A, the next Y minutes scored as A are rescored B”, where A and B mean sleep or wake, may be applied to the filtered signal afterward.

2.2.2. Sazonov Algorithm

Details of the Sazonov algorithm are described in Appendix A.4. Its output is given in terms of probability. At first glance, it appears to differ from the other discussed algorithms. Probability domain is $[0,1]$ and Sazonov algorithm is basing on it, however, from logit function definition $\mathrm{logit}\left(p\left[n\right]\right)=ln\frac{p\left[n\right]}{1-p\left[n\right]}\equiv \alpha \in [-\infty ,\infty ]$. In this case threshold of $p=0.5$ is equivalent to $\alpha =0$. This allows to easily eliminate the constant d from Equations (A6) and (A7) by adequately changing the threshold $\alpha$, and also makes the similarity to Cole-Kripke algorithm more visible. Therefore, Equation (A7) will be referred to as follows:

$$\alpha \left[n\right]=\sum _{i=0}^{8}A\left(i\right)max[n-i]$$

(6)

where A is the vector of coefficients, $max\left[n\right]$ is a maximum signal value in epoch n, $\alpha$ is the output of the algorithm. Equation (A6) should be understood in a similar fashion. This formulation includes steps (1.) and (2.) from Section 2.1, where $max\left[n\right]$ implements “collapsing into epochs”, which are further convolved with $A\left(i\right)$:

1.

$X\left[n\right]=max\left[n\right]$

2.

$\alpha \left[n\right]={\sum }_{i=0}^{8}A\left(i\right)X[n-i]$

Now we see clearly, that the Sazonov algorithm is a member of the Cole-Kripke algorithm family—except for the application of rescoring rules, which are treated as optional either way. It is also worth noting, that Equation (A6), which was not used in the end, similarly can be perceived as a linear combination of three different Cole-Kripke algorithms with a different resampling method used within each of them.

2.2.3. Sadeh Algorithm

At the beginning, it should be noted that both algorithms defined by Equations (A4) and (A5) contain an FIR component—for the first one it is an FIR filter of order 0, and for the second it is an FIR filter of order 10 combined with moving signal by 5 samples (similarly to the Cole-Kripke algorithm family). Mean from a few epochs can be represented as an FIR filter with all coefficients set to $\frac{1}{j+1}$, where j is the order of a filter. In both equations all components other than FIR filters are nonlinear. This allows to formulate the description of the Sadeh algorithm family as follows:

$$PS\left[n\right]=\sum _{i=-N}^{M}W\left(i\right)X[n-i]+\sum _{i=0}^{K}C\left(i\right){\varphi }_n\left(i\right)+d$$

(7)

where the first component is equivalent to a FIR filter, while the second one is a linear combination of some nonlinear transformations ${\varphi }_n\left(i\right)$. Finally, d is a constant and $N,M,K\in \mathbb{N}$.

As explained further in Section 4, we could treat some of the nonlinear parts of the normalized Sadeh algorithm from Equation (7) as some sort of rescoring rules. However, there are at least two significant obstacles preventing that. Firstly, it is extremely nontrivial to convert input variables like sample variances into rescoring rules akin to Cole-Kripke rescoring rules. Secondly, rescoring rules are applied on the output of the filter post hoc, not during filtering. Additionally, these nonlinear components are obligatory—contrary to the rescoring rules. Although a part of the Sadeh algorithm is an FIR filter and it does fit our unified framework, additional nonlinear terms of Equation (7) prevent us from analyzing the Sadeh algorithm through the lens of the unified Equation (1), as it would be incomplete and could lead to false conclusions.

3. Results

Properties of FIR Filters from the Cole-Kripke Algorithm Family

As presented in the previous section, all the discussed algorithms—except Sadeh’s—are equivalent to FIR filters. This simple observation allows for the first time to address the similarities and differences between these approaches in a unified and well-established framework. Furthermore, most importantly, the properties of FIR filters have been meticulously researched for decades, which greatly simplifies the physical interpretation of their properties.

This fact is both illustrated and summarized in Figure 1, presenting the transmittances (frequency responses) of the FIR filters corresponding step (2.) of different (versions of the) algorithms. Original coefficients, as published by the authors of respective algorithms, are recalculated as corresponding FIR filter coefficients, according to the mathematical transformations of algorithms’ equations which were described in the previous sections. Transmittances of the filters are computed using the Scipy module (v. 1.7.1) in Python [22] in a following way: for a vector of frequencies w response $H\left(e^{iw}\right)$ is returned, which in case of FIR filter of order M with coefficients $W\left(i\right)$ means: $W\left(0\right)+W\left(1\right)e^{-iw}+\dots +W\left(M\right)e^{-iwM}$. $\left|H\right|$ is an estimate of the transmittance.

By inspecting Figure 1 and

Table 1. Cutoff frequencies (at −3dB) of low-pass filters corresponding to step (2.) of analyzed algorithms. Third column presents the corresponding periods in minutes.

Algorithm	Low-Pass Cutoff	Period
Cole-Kripke (max 10 s nonoverlapping)	0.00127 Hz	13 m 08 s
Cole-Kripke (max 10 s overlapping)	0.00119 Hz	14 m 02 s
Cole-Kripke (max 30 s nonoverlapping)	0.00107 Hz	15 m 31 s
Cole-Kripke (mean 60 s)	0.00112 Hz	14 m 50 s
Webster (max 60 s overlapping)	0.00104 Hz	16 m 00 s
Webster (max 10 s overlapping)	0.00078 Hz	21 m 20 s
UCSD	0.00150 Hz	11 m 09 s
Sazonov	0.00153 Hz	10 m 53 s
Scripps Clinic	0.00143 Hz	11 m 38 s
Average:	$0.00121\pm 0.00023$ Hz	14 m 17 s ± 3 m 04 s

we observe that all of these filters are low-pass filters with cutoff frequencies corresponding to a period of approximately 10–20 min. In the following, we logically split the filtering into two formal stages: weakening and the actual filtering. By weakening, we mean the decrease of the signal’s amplitude by a constant, introduced in the process of filtering by multiplying all the filter’s coefficients by the same value.

We discuss only one version of the UCSD algorithm (Appendix A.2) because the only difference between its versions is the value of multiplicative constant P, while all the $W_i$ coefficients are always the same (see Equation (A1)). This means that the only difference between versions of UCSD algorithm is the dB value at which the passband begins, while the shape of the transmittance or cutoff frequency are always exactly the same. In other words, all those versions of the UCSD algorithm are the same filters with different weakening.

The largest differences in transmittances, as observed in Figure 1, are between different algorithms (i.e., Webster vs. Cole-Kripke vs. Scripps Clinic vs. Sazonov vs. UCSD) as opposed to versions of the same algorithms (i.e., different resampling methods used in the same algorithm). This is most likely due to the use of different devices in each study, or, possibly, other different properties of signal, like sampling frequency before resampling.

However, the consequences of the presented result are even more severe, as it seems that not all statements made in [21] regarding limitations of their algorithms are true. An important, yet partially misleading statement made by Cole et al. [21] was that for each different method of collapsing data into epochs, new algorithms—i.e., new sets of coefficients—should be derived. This seems to be only partially true as inspecting Figure 1 one can see, that transmittances of filters derived for following three resampling methods: maximum 10 s of each minute considering overlapping epochs, maximum 10 s of each minute considering non-overlapping epochs and maximum 30 s of each minute considering non-overlapping epochs (see [21]) are virtually identical. Therefore, there is no need for using different filters for these different resampling methods, as opposed to Cole et al. statement. While there is a difference in the weakening stage of the filter with data resampled using mean vs. using max metrics, there is little difference between different methods of downsampling data using max metric with different epoch lengths, overlapping or not. Differences between mean and max metrics used during resampling data as described in the Appendix A.1 are due to the fact that the mean will be less (or equal, but this is rather not a case in real data) than the maximum value of a set. Since no matter what resampling method was used, data is normalized so that $D\left[n\right]>1$ means wake, weakening of the mean method should be lesser than of the max because the maximum value of a set is greater than or equal to the mean.

From the above we may infer that:

1.

Discussed algorithms include a preprocessing stage, which consists of signal downsampling, but downsampling can have nontrivial sample selection. This procedure is not a standard resampling method and may introduce aliasing artifacts.

2.

Sample selection metric (e.g., maximum, average) changes the overall signal weakening constant value, but the overall low-pass filtering cutoff frequency remains the same.

3.

Weakening stage depends on certain hardware factors (like, e.g., sampling frequency).

4.

Modern filter design tools could be used to design a low-pass filter instead of an empirical coefficient fit. Such a potentially unified algorithm design could be easily adaptable to different actigraphy hardware and sampling rates.

4. Discussion

All of the discussed algorithms rely on resampling input data as the first step of analysis, usually via some non-standard procedures. It is worth mentioning that these procedures (originally called “collapsing”) are potential sources of aliasing artifacts, as Nyquist theorem requires low-pass filtering before downsampling to avoid them.

All of the presented algorithms, except Sadeh’s (1989) [12] algorithm, can be interpreted as low-pass filters. Cutoff frequencies are very similar in all cases and correspond to periods from 10 to 25 min. Treating the central step of these algorithms as FIR filters should allow making better use of the signal processing theory, as there are well-established and optimized packages for filter usage and design. Furthermore, constructing similar algorithms in the style analogous to the Cole-Kripke algorithm family should be trivial and quick, while maintaining the assumption that the signal filtered with a low-pass filter should indicate sleep below some threshold and wake above it (or vice versa).

However, similarities are even broader than merely stemming from the presence of a low-pass FIR component. For all the discussed approaches it seems that the main goal of the procedures is to smoothen the input signal. Not only by low-pass filtering but also by other nonlinear operations such as the application of rescoring rules. Their smoothing action is for example due to the fact that short peaks of activity during sleep time will be rescored as sleep.

Sadeh algorithms can be in part described by the unified mathematical framework we proposed, although it can be done only for the linear part of the equations. Similar to the other discussed algorithms, Sadeh (1994) [13] contains a low-pass filter—in a form of a moving average.

Most nonlinear terms of Equation (A4) (Sadeh from 1989) include standard deviation of signal epochs. Standard deviation can be used to detect edges in a signal [23]. These terms are used to detect signal transients, and reduce their influence on the overall output (as mentioned in Appendix A.3), effectively smoothing the output signal. Sadeh from 1994 (Equation (A5)) uses a moving average, which is effectively a low quality low-pass FIR filter (Figure 2). Additionally, one of the nonlinear terms is also standard deviation of signal epochs, and some nontrivial epoch thresholding. All these considerations support the claim that Sadeh algorithms, like all other discussed approaches, aim to smoothen the actigraphic signal before applying some fixed threshold to discriminate sleep and wake.

Accuracies of the discussed algorithms, provided both by Authors [12,13,15,18,21] as well as some independent studies [10,24,25], are very close. It can be explained by the fact that the algorithms are in their core FIR filters with very similar transmittances. The unified view proposed in this paper allows to understand and analyze their properties within a coherent and well-established framework, in which not only the mathematical formulae seem to be just specific cases of one procedure, but also their intended modes of operation (i.e., resampling and smoothing) seem to be the same.

Finally, presented observations may possibly lead to an informed construction of algorithms for processing and analysis of actigraphic data, based upon a firm mathematical and signal processing background, with better argumentation and much less effort, as opposed to the purely empirical approach employed so far.