(Multi-fractal)Detrended fluctuation analysis

(Multi-fractal)Detrended fluctuation analysis#

Import required libraries:

import circstudio as cs
import numpy as np
import plotly.graph_objects as go
import os
import plotly.io as pio
pio.templates.default = "plotly_white"
pio.renderers.default = "notebook"

Open a sample file:

# Create file path for sample files within the circStudio package
fpath = os.path.join(os.path.dirname(cs.__file__), 'data')

raw = cs.io.read_awd(
    os.path.join(fpath, 'example_01.AWD'), 
    start_time='1918-01-24 08:00:00',
    period='9 days'
)

Theoretical background#

Time series obtained using actigraphy (e.g., activity, light and temperature) often exhibit structure across multiple time scales: fluctuations measured over minutes are not independent from those measured over hours. (Multi-fractal) Detrended Fluctuation Analysis ((MF)DFA) quantifies this scale-invariant organization by asking:

  • How does the typical size of fluctuations change as we measure them over progressively larger time windows?

A key advantage of DFA is that it can detect long-range correlation (“a type of long-term memory in the signal”) structure in non-stationary signals (e.g., drifting baselines), because it explicitly removes local trends before quantifying fluctuations.

In the context of circadian biology, DFA/MF-DFA can help characterize whether activity dynamics show long-range organization, and whether that organization changes with aging, disease, or circadian disruption.

Overview of the method#

DFA/MF-DFA can be divided in three main steps:

  1. Integrate the signal into a profile after mean centering.

  2. Segment the profile into windows of length n and detrend each segment.

  3. Quantify fluctuation amplitude F(n) as a function of window size n

If the signal shows scale-invariant structure, the fluctuation function F(n) follows an approximate powerlaw over a range of scales n: $$F(n) \propto n^h$$

Taking logarithms yields an approximately linear relationship: $$\log F(n) \approx h \log n$$

The slope h is called the Hurst exponent in DFA (or the generalized Hurst exponent h(q) in MF-DFA).

1. Build the DFA profile#

DFA begins by removing the global mean and integrating the deviations: $$Y(t) = \sum_{i=1}^{N} (x_i - \overline{x})$$

This transformation converts the original time series into a cumulative trajectory, making slow trends and persistent structure easier to detect. In circStudio:

profile = cs.Fractal.profile(raw.activity.values)

2. Segment and detrend#

For each window size n, the profile Y(t) is segmented into consecutive non-overlaping windows. Within each segment, a polynomial of degree deg is fitted and removed (linear detrending corresponds to deg=1).

The remaining variability, after removing the local trend, reflects the intrinsic fluctuations of the signal at that specific timescale.

The local fluctuation magnitude is quantified using the mean squared residuals of the polynomial fit. Optionally, with overlap=True, consecutive windows overlap by 50% to improve resolution at the expense of additional computation.

In circStudio, segmentation is perform as follows:

# Segment actigraphy signal with a window size of 1000
segments = cs.Fractal.segmentation(profile, 1000)

Here, we use a window size (number of samples per segment) of 1000 samples. The result is a 2D array in which each row represents one segment of the integrated profile.

Computing local detrended fluctuations#

Next, we compute the mean squared residuals within each segment:

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

Each value in local_fluctuations quantifies how strongly the signal deviates from its local linear trend within a segment (in this case, of length 1000).

Averaging fluctuations (q=2, standard DFA)#

Finally, we average the fluctuation magnitudes across all segments:

cs.Fractal.q_th_order_mean_square(local_fluctuations,q=2)

np.float64(17613.535942976756)

3. Compute the fluctuation function F(n): standard DFA#

Standard DFA corresponds to MF-DFA evaluated at q = 2. In MF-DFA, the parameter q determines how strongly different types of fluctuations influence the overall fluctuation estimate (more details later once we dicuss MF-DFA). A q = 2 means that we will be treating fluctuations in a balanced manner.

Rather than computing the fluctuation amplitude for a single window size, we now repeat the procedure across a range of time scales n. This allows us to evaluate how the magnitude of fluctuations changes as we “zoom out” from short to long temporal windows.

Define the time scales#

We construct a set of window sizes n spanning ultradian to circadian scales:

# Generate array of time scales ranging from 10 min to 14400 min (24h)
n_array = np.geomspace(10, 1440, num=50, endpoint=True, dtype=int)

Using logarithmically spaced scales ensures that both short and long time scales are adequately represented.

Compute the fluctuation function#

We now compute F(n) for each window size n:

f_n = cs.Fractal.dfa(raw.activity, n_array, deg=1)

The output f_n contains the fluctuation amplitude corresponding to each scale in n_array. Grafically, the relationship is approximately linear:

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='log(n)',type='log'),
    yaxis=dict(title='log(F(n))')
)

Estimate the Hurst exponent#

The estimated slope corresponds to the Hurst exponent h, which characterizes how fluctuation amplitudes scale with time.

Conceptually, h reflects the degree of long-range temporal organization (“memory”) in the signal:

  • h = 0.5 indicates random, uncorrelated fluctuations (no long-term structure/”memory”).

  • h > 0.5 indicates persistent dynamics: periods of relatively high (or low) activity tend to be followed by similar periods across different scales.

  • h < 0.5 indicates anti-persistent dynamics: increases are more likely to be followed by decreases, reflecting more alternating behavior.

For actigraphy, h > 0.5 suggests that activity patterns are not random but temporally organized: variability observed over short intervals (e.g., minutes) is statistically related to the variability observed over longer intervals (e.g., hours).

In this sense, the Hurst exponent provides a compact measure of how strongly behavioral activity is structure across time scales.

We can estimate the Hurst exponent directly:

cs.Fractal.generalized_hurst_exponent(f_n, n_array, log=False)

(np.float64(0.996382744180578), np.float64(0.004809113976784806))

The first returned value is the estimated Hurst exponent, and the second is its standard error.

MF-DFA#

MF-DFA generalizes DFA by repeating the analysis across multiple values of q. As mentioned earlier, the role of q is to shift the emphasis of the fluctuation estimate:

  • Larger (positive) q values give more weight to segments with large deviations from the local trend (e.g., peaks, episodes of high activity).

  • Smaller (or negative) q values give more weight to segments with small deviations (quiet periods, subtle variability).

As a result, MF-DFA produces a family of fluctuation curves, F(n,q), rather than a single curve F(n)

Compute MF-DFA fluctuations across multiple q-orders#

# Specify q values (orders)
q_array = [1,2,3,4,5,6]

# Compute MF-DFA fluctuation function F(n,q)
mf_f_n = cs.Fractal.mfdfa(raw.activity, n_array, q_array, deg=1)

mf_f_n is a 2D array with shape (len(n_array), len(q_array)), where each column corresponds to one q-order.

Visualize scaling curves across q#

A useful first check is whether the curves are approximately linear across of scales.

# Visualize curves across q
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='log(n)',type='log'),
    yaxis=dict(title='log(F(n))')
)

Estimate the generalized Hurst exponent h(q)#

To summarize the MF-DFA results, we estimate the slope of each log F(n,q) versus log n relationship. This yields the generalized Hurst exponent function h(q):

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

[(np.float64(1.078789403267218), np.float64(0.005090112950945789)),
 (np.float64(0.996382744180578), np.float64(0.004809113976784806)),
 (np.float64(0.9498495984871702), np.float64(0.0049393201057842465)),
 (np.float64(0.9173765841355221), np.float64(0.005695765261621236)),
 (np.float64(0.8926070849788404), np.float64(0.006875225021534225)),
 (np.float64(0.8731126415786645), np.float64(0.008108469386111311))]

Next steps#

In the next tutorial, we will explore how to smooth actigraphy data using functional linear modeling, allowing us to capture underlying trends while reducing short-term variability.