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?