Automatic detection of rest/activity periods with pyActigraphy

The pyActigraphy package implements several rest/activity detection algorithms.

Broadly speaking, they can be divided in two categories:

  • Epoch-by-epoch rest/activity scoring (inspired by PSG):

    • Cole-Kripke’s, Sadeh’s, Scripps’ and Oakley’s algorithms

  • Detection of consolidated periods of similar activity patterns

    • Crespo’s, Roenneberg’s algorithms

All algorithms have been implemented to return a binary time serie (0 being rest or activity depending on the definition made in the original article…)

[1]:
import pyActigraphy
[2]:
import plotly.graph_objs as go
[3]:
import os
fpath = os.path.join(os.path.dirname(pyActigraphy.__file__),'tests/data/')
[4]:
raw = pyActigraphy.io.read_raw_rpx(
    fpath+'test_sample_rpx_eng.csv', start_time='2015-07-04 12:00:00', period='6 days', language='ENG_UK'
)

Epoch-by-epoch rest/activity scoring algorithms

All based on the same idea;

  • use a rolling window over the data and convolute them with a “gaussian”-like kernel.

  • if the resulting data is above a predefined threshold, classify as activity. Rest otherwise…

My opinion regarding these algorithms

The weights of these “gaussian”-like kernels as well as the applied thresholds are very much device and population dependant. Using these algorithms would thus require to adapt these parameters to the population under study and the actigraph device… Possible but seldom done.

[5]:
layout = go.Layout(title="Rest/Activity detection",xaxis=dict(title="Date time"), yaxis=dict(title="Counts/period"), showlegend=False)

Cole-Kripke

[6]:
help(raw.CK)
Help on method CK in module pyActigraphy.sleep.scoring_base:

CK(settings='30sec_max_non_overlap', threshold=1.0, rescoring=True) method of pyActigraphy.io.rpx.rpx.RawRPX instance
    Cole&Kripke algorithm for sleep-wake identification.

    Algorithm for automatic sleep scoring based on wrist activity,
    developped by Cole, Kripke et al [1]_.


    Parameters
    ----------
    settings: str, optional
        Data reduction settings for which the optimal parameters have been
        derived. Available settings are:

        * "mean": mean activity per minute
        * "10sec_max_overlap": maximum 10-second overlapping epoch per
          minute
        * "10sec_max_non_overlap": maximum 10-second nonoverlapping epoch
          per minute
        * "30sec_max_non_overlap": maximum 30-second nonoverlapping epoch
          per minute

        Default is "30sec_max_non_overlap".
    threshold: float, optional
        Threshold value for scoring sleep/wake. Default is 1.0.
    rescoring: bool, optional
        If set to True, Webster's rescoring rules are applied [2]_.
        Default is True.

    Returns
    -------
    ck: pandas.core.Series
        Time series containing the `D` scores (0: sleep, 1: wake) for each
        epoch.

    Notes
    -----

    The output variable D of the CK algorithm is defined in [1]_ as:

    .. math::

        D = P*(
            [W_{-4},\dots,W_{0},\dots,W_{+2}]
            \cdot
            [A_{-4},\dots,A_{0},\dots,A_{+2}])

    with:

    * D < 1 == sleep, D >= 1 == wake;
    * P, scale factor;
    * :math:`W_{0},W_{-1},W_{+1},\dots`, weighting factors for the present
      minute, the previous minute, the following minute, etc.;
    * :math:`A_{0},A_{-1},A_{+1},\dots`, activity scores for the present
      minute, the previous minute, the following minute, etc.


    .. warning:: This algorithm yields scores with a period of 1 minute.
                 Therefore, the time series returned by the CK function has
                 a 1-min period. In the original paper, this algorithm is
                 validated for devices configured in "zero-crossing mode".

    Webster's rescoring rules [2]_:

    1. After at least 4 minutes scored as wake, the next 1 minute scored as
       sleep is rescored as wake
    2. After at least 10 minutes scored as wake, the next 3 minutes scored
       as sleep are rescored as wake;
    3. After at least 15 minutes scored as wake, the next 4 minutes scored
       as sleep are rescored as wake;
    4. If a period of 6 minutes or less that is scored as sleep is
       surrounded by at least 15 minutes scored as wake, then rescore to
       wake;
    5. If a period of 10 minutes or less that is scored as sleep is
       surrounded by at least 20 minutes scored as wake, then rescore to
       wake.

    References
    ----------

    .. [1] Cole, R. J., Kripke, D. F., Gruen, W., Mullaney, D. J.,
           & Gillin, J. C. (1992). Automatic Sleep/Wake Identification
           From Wrist Activity. Sleep, 15(5), 461–469.
           http://doi.org/10.1093/sleep/15.5.461
    .. [2] Webster, J. B., Kripke, D. F., Messin, S., Mullaney, D. J., &
           Wyborney, G. (1982). An Activity-Based Sleep Monitor System for
           Ambulatory Use. Sleep, 5(4), 389–399.
           https://doi.org/10.1093/sleep/5.4.389

    Examples
    --------

[7]:
CK = raw.CK()
[8]:
layout.update(yaxis2=dict(title='Classification',overlaying='y',side='right'), showlegend=True);
[9]:
go.Figure(data=[
    go.Scatter(x=raw.data.index.astype(str),y=raw.data, name='Data'),
    go.Scatter(x=CK.index.astype(str),y=CK, yaxis='y2', name='CK')
], layout=layout)

Sadeh’s and Scripps’ algorithms

[10]:
sadeh = raw.Sadeh()
scripps = raw.Scripps()
[11]:
go.Figure(data=[
    go.Scatter(x=raw.data.index.astype(str),y=raw.data, name='Data'),
    go.Scatter(x=sadeh.index.astype(str),y=sadeh, yaxis='y2', name='Sadeh'),
    go.Scatter(x=scripps.index.astype(str),y=scripps, yaxis='y2', name='Scripps')
], layout=layout)

Oakley’s algorithm

This is the sleep/wake scoring algorithm used by the Actiware software from Respironics. The various activity thresholds have thus been tailored for Actiwatch devices. Be careful when applying this algorithm on data acquired with other devices.

[12]:
help(raw.Oakley)
Help on method Oakley in module pyActigraphy.sleep.scoring_base:

Oakley(threshold=40) method of pyActigraphy.io.rpx.rpx.RawRPX instance
    Oakley's algorithm for sleep/wake scoring.

    Algorithm for automatic sleep/wake scoring based on wrist activity,
    developed by Oakley [1]_.


    Parameters
    ----------
    threshold: float or str, optional
        Threshold value for scoring sleep/wake. Can be set to "automatic"
        (cf. Notes).
        Default is 40.

    Returns
    -------

    oakley: pandas.core.Series
        Time series containing scores (1: sleep, 0: wake) for each
        epoch.


    Notes
    -----

    The output variable O of Oakley's algorithm is defined as:

    .. math::

        O = (
            [W_{-2},W_{-1},W_{0},W_{+1},W_{+2}]
            \cdot
            [A_{-2},A_{-1},A_{0},A_{+1},A_{+2}])

    with:

    * O <= threshold == sleep, 0 > threshold == wake;
    * :math:`W_{0},W_{-1},W_{+1},\dots`, weighting factors for the present
      epoch, the previous epoch, the following epoch, etc.;
    * :math:`A_{0},A_{-1},A_{+1},\dots`, activity scores for the present
      epoch, the previous epoch, the following epoch, etc.


    The current implementation of this algorithm follows the description
    provided in the instruction manual of the Actiwatch Communication and
    Sleep Analysis Software (Respironics, Inc.) [2]_ :

    * 15-second sampling frequency:
        .. math::

            W_{-8} &= \ldots = W_{-5} = W_{+5} = \ldots = W_{+8} = 1/25 \\
            W_{-4} &= \ldots = W_{-1} = W_{+1} = \ldots = W_{+4} = 1/5 \\
            W_{0} &= 4

    * 30-second sampling frequency:
        .. math::

            W_{-4} &= W_{-3} = W_{+3} = W_{+4} = 1/25 \\
            W_{-2} &= W_{-1} = W_{+1} = W_{+2} = 1/5 \\
            W_{0} &= 2

    * 60-second sampling frequency:
        .. math::

            W_{-2} &= W_{+2} = 1/25 \\
            W_{-1} &= W_{+1} = 1/5 \\
            W_{0} &= 1

    * 120-second sampling frequency:
        .. math::

            W_{-1} &= W_{+1} = 1/8 \\
            W_{0} &= 1/2


    The *Automatic Wake Threshold Value* calculation is this [2]_:

    1. Sum the activity counts for all epochs of the data set.
    2. Count the number of epochs scored as MOBILE for the data set
       (the definition of MOBILE follows).
    3. Compute the MOBILE TIME (number of epochs scored as MOBILE from
       step 2 multiplied by the Epoch Length) in minutes.
    4. Compute the Auto Wake Threshold = ((sum of activity counts from
       step 1) divided by (MOBILE TIME from step 3)) multiplied by 0.88888.


    *Definition of Mobile* [2]_:

    An epoch is scored as MOBILE if the number of activity counts recorded
    in that epoch is greater than or equal to the epoch length in 15-second
    intervals. For example,there are four 15-second intervals for a
    1-minute epoch length; hence, the activity value in an epoch must be
    greater than, or equal to, four, to be scored as MOBILE.


    References
    ----------

    .. [1] Oakley, N.R. Validation with Polysomnography of the Sleepwatch
           Sleep/Wake Scoring Algorithm Used by the Actiwatch Activity
           Monitoring System; Technical Report; Mini-Mitter: Bend, OR, USA,
           1997
    .. [2] Instruction manual, Actiwatch Communication and Sleep Analysis
           Software
           (https://fccid.io/JIAAWR1/Users-Manual/USERS-MANUAL-1-920937)

Medium threshold

A threshold of 40 corresponds to a medium wake threshold.

[13]:
oakley = raw.Oakley(threshold=40)

Automatic threshold

It is also possible to compute a threshold value automatically, based on activity data.

[14]:
oakley_auto = raw.Oakley(threshold='automatic')
[15]:
go.Figure(data=[
    go.Scatter(x=raw.data.index.astype(str),y=raw.data, name='Data'),
    go.Scatter(x=oakley.index.astype(str),y=sadeh, yaxis='y2', name='Oakley (thr: medium)'),
    go.Scatter(x=oakley_auto.index.astype(str),y=scripps, yaxis='y2', name='Oakley (thr: automatic)')
], layout=layout)

Consolidated activity/rest period detection

Crespo’s algorithm

This is a threshold-based algorithm that used morphological filters to “clean” short periods of activity (rest) surrounded by periods of rest (acitivity).

[16]:
help(raw.Crespo)
Help on method Crespo in module pyActigraphy.sleep.scoring_base:

Crespo(zeta=15, zeta_r=30, zeta_a=2, t=0.33, alpha='8h', beta='1h', estimate_zeta=False, seq_length_max=100, verbose=False) method of pyActigraphy.io.rpx.rpx.RawRPX instance
    Crespo algorithm for activity/rest identification

    Algorithm for automatic identification of activity-rest periods based
    on actigraphy, developped by Crespo et al. [1]_.

    Parameters
    ----------
    zeta: int, optional
        Maximum number of consecutive zeroes considered valid.
        Default is 15.
    zeta_r: int, optional
        Maximum number of consecutive zeroes considered valid (rest).
        Default is 30.
    zeta_a: int, optional
        Maximum number of consecutive zeroes considered valid (active).
        Default is 2.
    t: float, optional
        Percentile for invalid zeroes.
        Default is 0.33.
    alpha: str, optional
        Average hours of sleep per night.
        Default is '8h'.
    beta: str, optional
        Length of the padding sequence used during the processing.
        Default is '1h'.
    estimate_zeta: bool, optional
        If set to True, zeta values are estimated from the distribution of
        ratios of the number of series of consecutive zeroes to
        the number of series randomly chosen from the actigraphy data.
        Default is False.
    seq_length_max: int, optional
        Maximal length of the aforementioned random series.
        Default is 100.
    verbose: bool, optional
        If set to True, print the estimated values of zeta.
        Default is False.

    Returns
    -------
    crespo : pandas.core.Series
        Time series containing the estimated periods of rest (0) and
        activity (1).

    References
    ----------

    .. [1] Crespo, C., Aboy, M., Fernández, J. R., & Mojón, A. (2012).
           Automatic identification of activity–rest periods based on
           actigraphy. Medical & Biological Engineering & Computing, 50(4),
           329–340. http://doi.org/10.1007/s11517-012-0875-y

    Examples
    --------

        >>> import pyActigraphy
        >>> rawAWD = pyActigraphy.io.read_raw_awd(fpath + 'SUBJECT_01.AWD')
        >>> crespo = rawAWD.Crespo()
        >>> crespo
        2018-03-26 14:16:00    1
        2018-03-26 14:17:00    0
        2018-03-26 14:18:00    0
        (...)
        2018-04-06 08:22:00    0
        2018-04-06 08:23:00    0
        2018-04-06 08:24:00    1
        Length: 15489, dtype: int64

[17]:
crespo = raw.Crespo()
/Users/ghammad/Work/pyActigraphy/pyActigraphy/sleep/scoring_base.py:1388: FutureWarning:

The `center` argument on `expanding` will be removed in the future

/Users/ghammad/Work/pyActigraphy/pyActigraphy/sleep/scoring_base.py:1392: FutureWarning:

The `center` argument on `expanding` will be removed in the future

[18]:
crespo_6h = raw.Crespo(alpha='6h')
[19]:
crespo_zeta = raw.Crespo(estimate_zeta=True)
[20]:
go.Figure(data=[
    go.Scatter(x=raw.data.index.astype(str),y=raw.data, name='Data'),
    go.Scatter(x=crespo.index.astype(str),y=crespo, yaxis='y2', name='Crespo'),
    go.Scatter(x=crespo_6h.index.astype(str),y=crespo_6h, yaxis='y2', name='Crespo (6h)'),
    go.Scatter(x=crespo_zeta.index.astype(str),y=crespo_zeta, yaxis='y2', name='Crespo (Automatic)')
], layout=layout)

Like for the other algorithms, the default parameters correspond to those described in the original papers, which, most likely, have been tuned to the population used to validate the algorithm.

Since the output is made of consolidated periods, it is possible to return the offset and onset times of each period:

[21]:
aot = raw.Crespo_AoT()
[22]:
aot
[22]:
(DatetimeIndex(['2015-07-04 13:00:00', '2015-07-05 06:04:30',
                '2015-07-06 05:18:30', '2015-07-07 05:38:00',
                '2015-07-08 05:30:00', '2015-07-09 05:36:30',
                '2015-07-10 05:51:00', '2015-07-10 12:00:00'],
               dtype='datetime64[ns]', freq=None),
 DatetimeIndex(['2015-07-04 12:00:30', '2015-07-04 21:59:00',
                '2015-07-05 22:02:30', '2015-07-06 21:08:00',
                '2015-07-07 21:54:00', '2015-07-08 21:17:30',
                '2015-07-09 21:33:30', '2015-07-10 11:00:30'],
               dtype='datetime64[ns]', freq=None))
[23]:
aot[0]-aot[1]
[23]:
TimedeltaIndex(['0 days 00:59:30', '0 days 08:05:30', '0 days 07:16:00',
                '0 days 08:30:00', '0 days 07:36:00', '0 days 08:19:00',
                '0 days 08:17:30', '0 days 00:59:30'],
               dtype='timedelta64[ns]', freq=None)

Roenneberg’s algorithm

Also threshold-based but uses correlations with test series of various lengths to find the consolidated period that best matches best the data.

[24]:
help(raw.Roenneberg)
Help on method Roenneberg in module pyActigraphy.sleep.scoring_base:

Roenneberg(trend_period='24h', min_trend_period='12h', threshold=0.15, min_seed_period='30Min', max_test_period='12h', r_consec_below='30Min', rsfreq=None) method of pyActigraphy.io.rpx.rpx.RawRPX instance
    Automatic sleep detection.

    Identification of consolidated sleep episodes using the
    algorithm developped by Roenneberg et al. [1]_.

    Parameters
    ----------
    trend_period: str, optional
        Time period of the rolling window used to extract the data trend.
        Default is '24h'.
    min_trend_period: str, optional
        Minimum time period required for the rolling window to produce a
        value. Values default to NaN otherwise.
        Default is '12h'.
    threshold: float, optional
        Fraction of the trend to use as a threshold for sleep/wake
        categorization.
        Default is '0.15'
    min_seed_period: str, optional
        Minimum time period required to identify a potential sleep onset.
        Default is '30Min'.
    max_test_period : str, optional
        Maximal period of the test series.
        Default is '12h'
    r_consec_below : str, optional
        Time range to consider, past the potential correlation peak when
        searching for the maximum correlation peak.
        Default is '30Min'.
    rsfreq: str, optional
        Resampling frequency used to evaluate the sleep periods. The final
        time series with rest/activity scores is returned with a frequency
        equal to one of the input data. If set to None, no resampling is
        performed.
        Default is None.

    Returns
    -------
    rbg : pandas.core.Series
        Time series containing the estimated periods of rest (1) and
        activity (0).

    Notes
    -----

    .. warning:: The performance of this algorithm has been evaluated for
                 actigraphy data aggregated in 10-min bins [2]_.

    References
    ----------

    .. [1] Roenneberg, T., Keller, L. K., Fischer, D., Matera, J. L.,
           Vetter, C., & Winnebeck, E. C. (2015). Human Activity and Rest
           In Situ. In Methods in Enzymology (Vol. 552, pp. 257-283).
           http://doi.org/10.1016/bs.mie.2014.11.028
    .. [2] Loock, A., Khan Sullivan, A., Reis, C., Paiva, T., Ghotbi, N.,
           Pilz, L. K., Biller, A. M., Molenda, C., Vuori‐Brodowski, M. T.,
           Roenneberg, T., & Winnebeck, E. C. (2021). Validation of the
           Munich Actimetry Sleep Detection Algorithm for estimating
           sleep–wake patterns from activity recordings. Journal of Sleep
           Research, April, 1–12. https://doi.org/10.1111/jsr.13371

    Examples
    --------

[25]:
roenneberg = raw.Roenneberg()
roenneberg_thr = raw.Roenneberg(threshold=0.25, min_seed_period='15min')
[26]:
go.Figure(data=[
    go.Scatter(x=raw.data.index.astype(str),y=raw.data, name='Data'),
    go.Scatter(x=roenneberg.index.astype(str),y=roenneberg, yaxis='y2', name='Roenneberg'),
    go.Scatter(x=roenneberg_thr.index.astype(str),y=roenneberg_thr, yaxis='y2', name='Roenneberg (Thr:0.25)')
], layout=layout)
[27]:
aot = raw.Roenneberg_AoT(threshold=0.25, min_seed_period='15min')
[28]:
aot
[28]:
(DatetimeIndex(['2015-07-04 19:34:30', '2015-07-05 06:58:00',
                '2015-07-05 16:02:30', '2015-07-06 06:59:00',
                '2015-07-06 19:32:00', '2015-07-07 07:22:30',
                '2015-07-07 16:18:30', '2015-07-07 22:00:30',
                '2015-07-08 07:06:30', '2015-07-08 21:40:00',
                '2015-07-09 07:14:30', '2015-07-10 07:24:00'],
               dtype='datetime64[ns]', name='Date_Time', freq=None),
 DatetimeIndex(['2015-07-04 19:16:00', '2015-07-04 21:14:00',
                '2015-07-05 14:28:00', '2015-07-05 20:25:00',
                '2015-07-06 18:15:00', '2015-07-06 20:17:30',
                '2015-07-07 15:20:00', '2015-07-07 19:03:30',
                '2015-07-07 22:37:30', '2015-07-08 19:13:30',
                '2015-07-08 22:01:30', '2015-07-09 20:29:30'],
               dtype='datetime64[ns]', name='Date_Time', freq=None))
[29]:
aot[0]-aot[1]
[29]:
TimedeltaIndex(['0 days 00:18:30', '0 days 09:44:00', '0 days 01:34:30',
                '0 days 10:34:00', '0 days 01:17:00', '0 days 11:05:00',
                '0 days 00:58:30', '0 days 02:57:00', '0 days 08:29:00',
                '0 days 02:26:30', '0 days 09:13:00', '0 days 10:54:30'],
               dtype='timedelta64[ns]', name='Date_Time', freq=None)

Et voilà! Easy, isn’t it?