Cosinor analysis with pyActigraphy

This notebook illustrates how to perform a single-component Cosinor analysis on actigraphy data.

The idea of a Cosinor analysis is to estimate some key parameters of the actigraphy count series by fitting these data with a (co)sine curve with a period of 24h:

\[Y(t) = M + A*cos(\frac{2\pi}{T}*t + \phi)\]

where M is the MESOR (Midline Statistic Of Rhythm), A is the amplitude of the oscillations, T is the period and \(\phi\) is the acrophase.

Then, the fit procedure provides estimates of these parameters which can then help to characterize the 24h rest-activity rhythm of an individual.

However, despite its popularity in circadian rhythm analyses, the Cosinor analysis exhibits some short-comings when dealing with non-stationary signals such as actigraphy signals. Intuitively, such model with a fixed period, amplitude and offset cannot accomodate signals with:

  • rest-activity cycles that do not follow a perfect 24h period;

  • an overall daily activity that changes on a day-by-day basis;

  • an overall activity trend that evolves with time.

For more informations on the topic:

Leise, T. L., & Harrington, M. E. (2011). Wavelet-Based Time Series Analysis of Circadian Rhythms. Journal of Biological Rhythms, 26(5), 454–463. https://doi.org/10.1177/0748730411416330

For a nice comparison of several techniques, including Cosinor and SSA (available in pyActigraphy):

Fossion, R., Rivera, A. L., Toledo-Roy, J. C., Ellis, J., & Angelova, M. (2017). Multiscale adaptive analysis of circadian rhythms and intradaily variability: Application to actigraphy time series in acute insomnia subjects. PLOS ONE, 12(7), e0181762. https://doi.org/10.1371/journal.pone.0181762

Import packages and input data

First, load your favourite python packages:

[1]:

import pyActigraphy

[2]:

import plotly.graph_objects as go

Let us now define the path to the example data files contained in the pyActigraphy package:

[3]:

import os
fpath = os.path.join(os.path.dirname(pyActigraphy.__file__),'tests/data/')

Cosinor analysis

Now, simply import the Cosinor module:

[4]:

from pyActigraphy.analysis import Cosinor

… and define a Cosinor object:

[5]:

cosinor = Cosinor()

Initial fit values

By default, the initial values of the cosine fit functions are the following:

[6]:

cosinor.fit_initial_params.pretty_print()

Name          Value      Min      Max   Stderr     Vary     Expr Brute_Step
Acrophase     3.142        0    6.283     None     True     None     None
Amplitude        50        0      inf     None     True     None     None
Mesor            50        0      inf     None     True     None     None
Period         1440        0      inf     None     True     None     None

These parameters will be estimated from the data. However, in order to ease the fit procedure, it makes sense to start from reasonable values. For example, the initial value for the period is set to 1440, the number of minutes per day. It makes only sense for data with a 1-min sampling frequency.

However, it is trivial to change these intial values.

For example, let’s change the initial value of the period to 2880, in case you would like to analyse data with a 30-sec sampling frequency:

[7]:

cosinor.fit_initial_params['Period'].value = 2880

[8]:

cosinor.fit_initial_params.pretty_print()

Name          Value      Min      Max   Stderr     Vary     Expr Brute_Step
Acrophase     3.142        0    6.283     None     True     None     None
Amplitude        50        0      inf     None     True     None     None
Mesor            50        0      inf     None     True     None     None
Period         2880        0      inf     None     True     None     None

Free and fixed fit parameters

By default, none of the fit parameters are fixed and will be estimated from the data. However, in some cases, it might be convenient to fix the initial of a parameter and not let it vary during the fit procedure.

It might be the case of the period for instance. Let us fix it to 1440.

[9]:

cosinor.fit_initial_params['Period'].value = 1440

[10]:

cosinor.fit_initial_params['Period'].vary = False

[11]:

cosinor.fit_initial_params.pretty_print()

Name          Value      Min      Max   Stderr     Vary     Expr Brute_Step
Acrophase     3.142        0    6.283     None     True     None     None
Amplitude        50        0      inf     None     True     None     None
Mesor            50        0      inf     None     True     None     None
Period         1440        0      inf     None    False     None     None

It is also possible to apply lower and upper bounds to these parameters.

For more informations about these parameters, please see: https://lmfit.github.io/lmfit-py/parameters.html

Analysis of a single subject

Read data from an example file…

[12]:

raw = pyActigraphy.io.read_raw_awd(fpath+'example_01.AWD', start_time='1918-01-24 08:30:00', period="7D")

… and plot it for visual inspection:

[13]:

fig = go.Figure(go.Scatter(x=raw.data.index.astype(str),y=raw.data))

[14]:

fig.show()

To perform a cosinor analysis with the pyActigraphy package, it is as simple as:

[15]:

results = cosinor.fit(raw.data, verbose=True) # Set verbose to True to print the fit output

[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 28
    # data points      = 10081
    # variables        = 3
    chi-square         = 6.8103e+08
    reduced chi-square = 67576.1845
    Akaike info crit   = 112113.910
    Bayesian info crit = 112135.566
[[Variables]]
    Amplitude:  156.998940 +/- 3.66168149 (2.33%) (init = 50)
    Acrophase:  4.86541767 +/- 0.02332300 (0.48%) (init = 3.141593)
    Period:     1440 (fixed)
    Mesor:      174.191952 +/- 2.58907736 (1.49%) (init = 50)

Fit results

To access the best fit parameter values:

[16]:

results.params['Mesor'].value

[16]:

174.19195195975644

It is also possible to transform them to a dictionary:

[17]:

results.params.valuesdict()

[17]:

OrderedDict([('Amplitude', 156.99894000312472),
             ('Acrophase', 4.865417670453099),
             ('Period', 1440),
             ('Mesor', 174.19195195975644)])

The results object contains also informations about the goodness-of-fit:

[18]:

results.aic # Akaike information criterium

[18]:

112113.91042765231

[19]:

results.redchi # Reduced Chi^2

[19]:

67576.1844639616

More informations about the list of the informations accessible via the results object: https://lmfit.github.io/lmfit-py/fitting.html#minimizerresult-the-optimization-result

Fit result visualization

If you are interested in visualizing the cosinor fit that best approximates the data, the Cosinor class implements a convinient function:

[20]:

help(cosinor.best_fit)

Help on method best_fit in module pyActigraphy.analysis.cosinor:

best_fit(ts, params) method of pyActigraphy.analysis.cosinor.Cosinor instance
    Best fit function of the data.

    Parameters
    ----------
    ts : pandas.Series
        Originally fitted time series.
    params: instance of Parameters [1]_
        Best fit parameters.

    Returns
    -------
    bestfit_data : pandas.Series
        Time series of the best fit data.

    References
    ----------

    .. [1] Non-Linear Least-Squares Minimization and Curve-Fitting for
           Python.
           https://lmfit.github.io/lmfit-py/index.html


[21]:

best_fit = cosinor.best_fit(raw.data, results.params)

[22]:

fig = go.Figure(
    data=[
        go.Scatter(x=raw.data.index.astype(str),y=raw.data,name='Raw data'),
        go.Scatter(x=best_fit.index.astype(str),y=best_fit,name='Best fit')
    ]
)

[23]:

fig.show()

Batch analysis

While it is sometimes necessary to perform a cosinor analysis of a single subject to get used with the procedure and get “a feeling” about the data, it is often not pratical to follow the same procedure to analyse multiple subjects.

Fortunately, pyActigraphy implements several convenient functions to ease a batch analysis.

Let us first read a batch of files, using the example files included the pyActigraphy package for demonstration purposes:

[24]:

readers = pyActigraphy.io.read_raw(fpath+'example_*.AWD', reader_type='AWD')

[25]:

len(readers.readers)

[25]:

6

6 files have been found and read by the readers object.

Read and apply a SST log (See https://ghammad.github.io/pyActigraphy/pyActigraphy-SSt-log.html for more informations):

[26]:

readers.read_sst_log(fpath+'example_sstlog.csv')

[27]:

readers.apply_sst(verbose=True)

Could not find an entry in SST log file for example_04
Could not find an entry in SST log file for example_05
Could not find an entry in SST log file for example_02
Found an entry in SST log file for example_03
Found an entry in SST log file for example_01
Could not find an entry in SST log file for example_01_mask

Now, let perform a cosinor analysis on each of these files, in parallel, using 3 cpu’s in order to speed up the computations.

[28]:

results_batch = cosinor.fit_reader(readers, n_jobs=3, prefer='threads') # prefer='threads': add that parameter if running on Mac OS.

For convenience, the output results are formatted as a Pandas.DataFrame (https://pandas.pydata.org/pandas-docs/stable/getting_started/dsintro.html#dataframe) in order to ease further analysis.

[29]:

results_batch

[29]:
Amplitude Acrophase Period Mesor BIC RedChiSq
example_04 7.345742e+01 0.927315 1440 81.733231 321025.485636 28448.100307
example_05 1.352039e+02 5.472072 1440 120.821941 225115.967014 31930.970729
example_02 5.711733e-08 3.032968 1440 183.388579 207392.482663 77803.026821
example_03 3.381232e+02 6.240036 1440 346.219095 171305.166394 146324.160952
example_01 1.577049e+02 4.750332 1440 171.007640 143499.394878 64194.285808
example_01_mask 1.266343e-07 2.986621 1440 127.565387 202689.738565 60700.671139

To write these results to an output file, it is as simple as: results_batch.to_csv('myfile.csv').

More info on to_csv: https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.to_csv.html?highlight=to_csv#pandas-dataframe-to-csv

Et voilà! Easy, isn’t it?