Track lost degrees of freedom

[tsalo]

Summary

afni_proc.py does this, which is awesome.

[tsalo]

Here's some code I wrote a while back that can be expanded and adapted for XCP-D:

"""Check the effect of our proposed processing settings on the degrees of freedom."""

import numpy as np

# Constants
TR = 1.725
max_time = 450  # 7:30 total run length
min_time = 180  # 3:00 post-censoring minimum run length
max_volumes = int(np.ceil(max_time / TR))  # 261
min_volumes = int(np.ceil(min_time / TR))  # 105
n_confounds = 36

# Calculate the frequencies that the time series samples
fs = 1 / TR
nyq = 0.5 * fs
spacing = 1 / max_time
n = int(np.ceil(nyq / spacing))
spacing2 = nyq / (n - 1)  # replicate what periodogram gives us
frequencies_hz = np.linspace(0, nyq, n)
# Can roughly reproduce the frequencies here by simulating a time series
# of max_volumes length and running scipy.signal.periodogram.
#
# from scipy import signal
# ts = np.random.random(max_volumes)
# f, p = signal.periodogram(ts, fs=fs)
# assert np.allclose(f, frequencies_hz, atol=0.005)  # True

# Figure out what the change in DOF is from the bandpass filter
highpass, lowpass = 0.01, 0.08
dropped_freqs_idx = np.where((frequencies_hz < highpass) | (frequencies_hz > lowpass))[0]
n_dropped_freqs = dropped_freqs_idx.size

# Calculate the remaining degrees of freedom
# The 2 comes from the df formula for correlations: n - 2
dof = min_volumes - (n_confounds + n_dropped_freqs + 2)
print(dof)  # -33

@erikglee noted that filtfilt is a second-order filter (i.e., it applies the filter forwards and backwards) so I need to modify the effect of filtering on the degrees of freedom.

[tsalo]

Here's a function that seems to work for just calculating the impact of the filter:

def calculate_dof(n_volumes, t_r, high_pass=0, low_pass=np.inf):
    """Calculate the number of degrees of freedom lost by a temporal filter.
    
    Parameters
    ----------
    n_volumes : int
        Number of data points in the time series.
    t_r : float
        Repetition time of the time series, in seconds.
    high_pass : float
        High-pass frequency in Hertz. Default is 0 (no high-pass filter).
    low_pass : float or numpy.inf
        Low-pass frequency in Hertz. Default is np.inf (no low-pass filter).
    
    Returns
    -------
    dof_lost : int
        Number of degrees of freedom lost by applying the filter.
    
    Notes
    -----
    Both Caballero-Gaudes & Reynolds (2017) and Reynolds et al. (preprint)
    say that each frequency removed drops two degrees of freedom.
    """
    import numpy as np

    duration = t_r * n_volumes
    fs = 1 / t_r
    nyq = 0.5 * fs
    spacing = 1 / duration
    n_freqs = int(np.ceil(nyq / spacing))
    frequencies_hz = np.linspace(0, nyq, n_freqs)

    # Figure out what the change in DOF is from the bandpass filter
    dropped_freqs_idx = np.where((frequencies_hz < high_pass) | (frequencies_hz > low_pass))[0]
    n_dropped_freqs = dropped_freqs_idx.size

    # Calculate the lost degrees of freedom
    dof_lost = n_dropped_freqs * 2
    return dof_lost

[tsalo]

I've opened https://neurostars.org/t/effect-of-temporal-filter-order-on-degrees-of-freedom/30818 to see if the experts (Rick Reynolds and Cesar Caballero-Gaudes) know the answer regarding filter order.