(Multi-fractal)Detrended fluctuation analysis

Human activity exhibits a temporal organization characterised by scale-invariant (fractal) patterns over time scales ranging from minutes to 24 hours. The (MF)DFA method allows the quantification of this scale-invariance.

The (MF)DFA methods have been originally described in:

• Peng, C.-K., Buldyrev, S. V., Havlin, S., Simons, M., Stanley, H. E., & Goldberger, A. L. (1994). Mosaic organization of DNA nucleotides. Physical Review E, 49(2), 1685–1689. https://doi.org/10.1103/PhysRevE.49.1685

• Kantelhardt, J. W., Zschiegner, S. A., Koscielny-Bunde, E., Havlin, S., Bunde, A., & Stanley, H. E. (2002). Multifractal detrended fluctuation analysis of nonstationary time series. Physica A: Statistical Mechanics and Its Applications, 316(1–4), 87–114. https://doi.org/10.1016/S0378-4371(02)01383-3

Aging and Alzheimer’s disease, both marked by an alteration of the suprachiasmatic nucleus (SCN), the circadian pacemaker, have been associated with a degradation of this scale-invariant organization. More informations in:

• Hu, K., Van Someren, E. J. W., Shea, S. A., & Scheer, F. A. J. L. (2009). Reduction of scale invariance of activity fluctuations with aging and Alzheimer’s disease: Involvement of the circadian pacemaker. Proceedings of the National Academy of Sciences, 106(8), 2490–2494. https://doi.org/10.1073/pnas.0806087106

• Li, P., Yu, L., Lim, A. S. P., Buchman, A. S., Scheer, F. A. J. L., Shea, S. A., … Hu, K. (2018). Fractal regulation and incident Alzheimer’s disease in elderly individuals. Alzheimer’s & Dementia, 14(9), 1114–1125. https://doi.org/10.1016/j.jalz.2018.03.010

• Li, P., Yu, L., Yang, J., Lo, M. T., Hu, C., Buchman, A. S., … Hu, K. (2019). Interaction between the progression of Alzheimer’s disease and fractal degradation. Neurobiology of Aging, 83, 21–30. https://doi.org/10.1016/j.neurobiolaging.2019.08.023

DFA steps

In a nutshell, the DFA method consists in:

1. removing the global mean and integrating the time series of a signal, that is:

   \[X_{t} = \sum_i^N(x_i - \bar{x})\]

   where \(\bar{x}\) denotes the mean value of the time series \(\{x_i\}_{i\in[1:N]}\)

2. dividing the integrated signal into nonoverlapping windows of the same chosen size n;

3. detrending the integrated signal in each window using polynomial functions to obtain residuals, that is:

   \[\widehat{X_t} = X_{t} - Y_{t}\]

   where \(Y_t\) denotes the trend obtained by polynomial fit and \(\widehat{X_t}\) the integrated time series after detrending;

4. calculating the root mean square of residuals in all windows as detrended fluctuation amplitude F(n), that is:

   \[F_n = \sqrt{\frac{1}{N} \sum_{t=1}^N \widehat{X_t}^2}\]

Imports and Input data

First, import your favourite packages:

import pyActigraphy

import numpy as np

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

import plotly.graph_objects as go
from plotly.subplots import make_subplots

# The DFA methods are part of the Fractal module:
from pyActigraphy.analysis import Fractal

help(Fractal.dfa)

Help on method dfa in module pyActigraphy.analysis.fractal:

dfa(ts, n_array, deg=1, overlap=False, log=False) method of builtins.type instance
    Detrended Fluctuation Analysis function

    Compute the q-th order mean squared fluctuations for different segment
    lengths.

    Parameters
    ----------
    ts: pandas.Series
        Input signal.
    n_array: array of int
        Time scales (i.e window sizes). In minutes.
    deg: int, optional
        Degree(s) of the fitting polynomials.
        Default is 1.
    overlap: bool, optional
        If set to True, consecutive windows during segmentation
        overlap by 50%. Default is False.
    log: bool, optional
        If set to True, returned values are log-transformed.
        Default is False.

    Returns
    -------
    q_th_order_msq_fluc: numpy.array
        Array of q-th order mean squared fluctuations.

Read test file:

raw = pyActigraphy.io.read_raw_awd(fpath+'example_01.AWD', start_time='1918-01-24 08:00:00', period='9 days')

raw.duration()

Timedelta('9 days 00:01:00')

In order to illustrate how the method works, let’s review its different steps:

Signal detrending and integration

profile = Fractal.profile(raw.data.values)

# Create figure with secondary y-axis
fig = make_subplots(specs=[[{"secondary_y": True}]])

# Add traces
fig.add_trace(
    go.Scatter(x=raw.data.index.astype(str),y=raw.data.values, name='Data'),
    secondary_y=False,
)
fig.add_trace(
    go.Scatter(x=raw.data.index.astype(str),y=profile, name='Profile'),
    secondary_y=True,
)

# Add figure title
fig.update_layout(
    title_text="Detrended and integrated data profile"
)

# Set x-axis title
fig.update_xaxes(title_text="Date time")
# Set y-axes titles
fig.update_yaxes(title_text="Activity counts", secondary_y=False)
#fig.update_yaxes(title_text="<b>secondary</b> yaxis title", secondary_y=True)

fig.show()

By definition, the profile values at the beginning and the end of the recording is zero, as shown on the plot.

As a side remark, it is interesting to note how the total daily activity fluctuates; the mean activity of the first three days is significantly lower than the activity of the remaining days.

Signal segmentation into non-overlapping windows

The next step of the DFA method consists to divide the signal into consecutive non-overlapping segments of equal length.

help(Fractal.segmentation)

Help on method segmentation in module pyActigraphy.analysis.fractal:

segmentation(x, n, backward=False, overlap=False) method of builtins.type instance
    Segmentation function

    Segment the signal into non-overlapping windows of equal size.

    Parameters
    ----------
    x: numpy.array
        Input array.
    n: int
        Window size.
    backward: bool, optional
        If set to True, start segmentation from the end of the signal.
        Default is False.
    overlap: bool, optional
        If set to True, consecutive windows overlap by 50%.
        Default is False.

    Returns
    -------
    segments: numpy.array
        Non-overlappping windows of size n.

# Example of segmentation with a window size of 1000 elements.
n = 1000
segments = Fractal.segmentation(profile, n)

len(profile)

12961

segments.shape

(12, 1000)

As expected, the profile has been segmented into 12 windows of size 1000. The remaining points are discarded. In order to avoid any bias from these discarded points, the DFA method, implemented in pyActigraphy averages over the results obtained over segments made from the start (foward) and from the end (backward) of the recording.

Local and global q-th order fluctuations

help(Fractal.local_msq_residuals)

Help on method local_msq_residuals in module pyActigraphy.analysis.fractal:

local_msq_residuals(segment, deg) method of builtins.type instance
    Mean squared residuals

    Mean squared residuals of the least squares polynomial fit.

    Parameters
    ----------
    segment: numpy.array
        Input array.
    deg: int
        Degree(s) of the fitting polynomials.

    Returns
    -------
    residuals_msq: numpy.float
        Mean squared residuals.

In the next step, the local fluctuations (i.e mean squared residuals of a least-square fit using a polynomial of order 1) are computed for each segment:

local_fluctuations = [Fractal.local_msq_residuals(segment,deg=1) for segment in segments]

Finally, the local fluctuations are averaged over all segments:

help(Fractal.q_th_order_mean_square)

Help on method q_th_order_mean_square in module pyActigraphy.analysis.fractal:

q_th_order_mean_square(F, q) method of builtins.type instance
    Qth-order mean square function

    Compute the q-th order mean squares.

    Parameters
    ----------
    F: numpy.array
        Array of fluctuations.
    q: scalar
        Order.

    Returns
    -------
    qth_msq: numpy.float
        Q-th order mean square.

The DFA method corresponds to \(q=2\):

Fractal.q_th_order_mean_square(local_fluctuations,q=2)

17613.53594297675

This value correspond the global fluctuation observed at a time scale of \(n=1000\). The previous steps are repeated for various values of \(n\) in order to estimate the function \(F_q(n)\).

In case of scale-invariance, this function should scale as:

\[F_q(n) \propto n^{h(q)}\]

with \(h(q)\) the generalized scaling (or Hurst) exponent.

DFA

In addition to the various functions illustrated above, the pyActigraphy implements a global function that performs all the necessary steps to carry out a DFA analysis:

help(Fractal.dfa)

Help on method dfa in module pyActigraphy.analysis.fractal:

dfa(ts, n_array, deg=1, overlap=False, log=False) method of builtins.type instance
    Detrended Fluctuation Analysis function

    Compute the q-th order mean squared fluctuations for different segment
    lengths.

    Parameters
    ----------
    ts: pandas.Series
        Input signal.
    n_array: array of int
        Time scales (i.e window sizes). In minutes.
    deg: int, optional
        Degree(s) of the fitting polynomials.
        Default is 1.
    overlap: bool, optional
        If set to True, consecutive windows during segmentation
        overlap by 50%. Default is False.
    log: bool, optional
        If set to True, returned values are log-transformed.
        Default is False.

    Returns
    -------
    q_th_order_msq_fluc: numpy.array
        Array of q-th order mean squared fluctuations.

Let’s first construct an array of time scales (in minutes, as mentioned in the “help” documentation), at which we would like to evaluate the fluctuations:

n_array = np.geomspace(10, 1440, num=50, endpoint=True, dtype=int) # Numbers spaced evenly on a log scale, ranging from an ultradian time scale (10 min.) to a circadian one (1440 min, i.e. 24h)

Then, we calculate the associated fluctuations:

F_n = Fractal.dfa(raw.data,n_array,deg=1)

fig = go.Figure(data=[
    go.Scatter(x=n_array,y=np.log(F_n), name='Data fluctuation',mode='markers+lines')])
fig.update_layout(
    height=800, width=800,
    xaxis=dict(title='Time (min.)',type='log'),
    yaxis=dict(title='log(F(n))')
)

Since human locomoter activity is characterized by a scale-invariant pattern, the log-log plot of \(F_q(n)\) vs \(n\) shows a linear behaviour.

Finally, the generalized Hurst exponent can be extracted using:

Fractal.generalized_hurst_exponent(F_n, n_array, log=False)

(0.996382744180578, 0.004809113976784807)

The generalized Hurst exponent for this recording is:

\[h(q=2) = 0.996 \pm 0.005\]

Multifractal DFA

Briefly, the MFDFA method extends the DFA method by evaluating the local fluctuations at various q-orders (instead of just q=2 for DFA). This results in a series of series of fluctuations.

In pyActigraphy, similarly to the DFA, a MFDFA can be performed with a single function:

help(Fractal.mfdfa)

Help on method mfdfa in module pyActigraphy.analysis.fractal:

mfdfa(ts, n_array, q_array, deg=1, overlap=False, log=False) method of builtins.type instance
    Multifractal Detrended Fluctuation Analysis function

    Compute the q-th order mean squared fluctuations for different segment
    lengths and different index q values.

    Parameters
    ----------
    ts: pandas.Series
        Input signal.
    n_array: array of int
        Time scales (i.e window sizes). In minutes.
    q_array: array of float
        Orders of the mean squares.
    deg: int, optional
        Degree(s) of the fitting polynomials.
        Default is 1.
    overlap: bool, optional
        If set to True, consecutive windows during segmentation
        overlap by 50%. Default is False.
    log: bool, optional
        If set to True, returned values are log-transformed.
        Default is False.

    Returns
    -------
    q_th_order_msq_fluctuations: numpy.array
        (n,q) array of q-th order mean squared fluctuations.

Let’s now define an array of q values:

q_array = [1,2,3,4,5,6]

MF_F_n = Fractal.mfdfa(raw.data,n_array,q_array,deg=1)

fig = go.Figure(data=[
    go.Scatter(x=n_array,y=np.log(MF_F_n[:,q]), name='Data fluctuation (q-th order: {})'.format(q_array[q]),mode='markers+lines') for q in range(len(q_array))])
fig.update_layout(
    height=800, width=800,
    xaxis=dict(title='Time (min.)',type='log'),
    yaxis=dict(title='log(F(n))')
)

mf_h_q = [Fractal.generalized_hurst_exponent(MF_F_n[:,q],n_array) for q in range(len(q_array))]

mf_h_q

[(1.0787894032672178, 0.005090112950947376),
 (0.996382744180578, 0.004809113976784807),
 (0.94984959848717, 0.004939320105785093),
 (0.9173765841355224, 0.005695765261620553),
 (0.8926070849788404, 0.006875225021534225),
 (0.8731126415786645, 0.008108469386111311)]

For q = 2, the result previously obtained with a DFA is recovered.

Et voilà!