Analyses

PSD, intra- and inter-brain measures functions

Option	Description
title	analyses.py
authors	Phoebe Chen, Florence Brun, Guillaume Dumas
date	2020-03-18

behav_corr(data, behav, data_name, behav_name, p_thresh, multiple_corr=True, verbose=False)

Correlates neural data with behavioral parameters using appropriate linear correlation methods.

This function first checks for data normality, then applies either Pearson's or Spearman's correlation. For connectivity data, it can apply multiple comparison correction.

Parameters

data : np.ndarray Data to correlate with behavior. Can be either: - 1D array of values (same shape as behav) - 3D array of connectivity values with shape (n_dyads, n_channels, n_channels)

np.ndarray

1D array of behavioral values (e.g., timing, performance scores)

str

Description of the data (used for plot labeling if verbose=True)

str

Description of the behavioral parameter (used for plot labeling if verbose=True)

float

Significance threshold for correlation tests (typically 0.05)

bool, optional

Whether to apply multiple comparison correction (default=True). Uses the FDR-BH method when True.

bool, optional

Whether to generate visualization of the correlation (default=False)

Returns

corr_tuple : namedtuple A named tuple containing: - r: correlation coefficient(s) - For 1D data: a single float - For connectivity data: array of shape (n_channels, n_channels) - pvalue: p-value(s) (original or corrected if multiple_corr=True) - strat: string indicating which strategy was used ('normal', 'non_normal', or information about multiple comparison correction)

Notes

When correlating connectivity matrices with behavioral data, only the significant correlations are returned in the r matrix, with non-significant values set to zero.

Examples

Correlating average alpha power with reaction times

result = behav_corr(alpha_power, reaction_times, 'Alpha Power', 'RT (ms)', 0.05) r_value = result.r p_value = result.pvalue

compute_conn_mvar(complex_signal, mvar_params, ica_params, measure_params, check_stability=True)

Computes connectivity measures based on multivariate autoregressive (MVAR) modeling.

This function fits MVAR models to the data and computes directed connectivity measures based on those models, with optional ICA for source separation.

Parameters

complex_signal : np.ndarray Complex analytic signals with shape (2, n_epochs, n_channels, n_freq_bins, n_times)

dict

Dictionary of MVAR model parameters with keys: - 'mvar_order': int, order of the MVAR model - 'fitting_method': str, method for fitting the MVAR model - 'delta': float, regularization parameter

dict

Dictionary of ICA parameters with keys: - 'method': str, ICA algorithm to use - 'random_state': int, random seed for reproducibility

dict

Dictionary defining the connectivity measure to compute with keys: - 'name': str, name of the connectivity measure (e.g., 'PDC', 'DTF') - 'n_fft': int, number of FFT points

bool, optional

Whether to verify the stability of the MVAR model (default=True) If True, the function will check if the model is stable and prompt for continuation or epoch merging if necessary.

Returns

connectivity : np.ndarray Connectivity measure matrix with shape: (n_epochs, n_freq, n_channels, n_channels, n_fft) or (1, n_freq, n_channels, n_channels, n_fft) if epochs are merged

Notes

MVAR-based connectivity measures provide information about directed (causal) relationships between signals, unlike most of the other connectivity metrics implemented in this module.

The function relies on the MVARICA approach: 1. Apply ICA to separate independent sources 2. Fit an MVAR model to the sources 3. Compute connectivity measures in the source space 4. Back-project to the original signal space

For stability, the number of time samples should be substantially larger than: mvar_order × n_channels^2

Warnings

This function is computationally intensive, especially for high model orders and large numbers of channels.

If the MVAR model is unstable, the function will offer to merge epochs to increase the sample size for more stable estimation.

References

Baccalá, L. A., & Sameshima, K. (2001). Partial directed coherence: a new concept in neural structure determination. Biological cybernetics, 84(6), 463-474.

Examples

mvar_params = {'mvar_order': 5, 'fitting_method': 'least_squares', 'delta': 0.0} ica_params = {'method': 'extended-infomax', 'random_state': 42} measure_params = {'name': 'PDC', 'n_fft': 512} conn = compute_conn_mvar(complex_signal, mvar_params, ica_params, measure_params)

compute_freq_bands(data, sampling_rate, freq_bands, filter_signal=True, **filter_options)

Computes analytic signals for specified frequency bands using FIR filtering and Hilbert transform.

This function performs bandpass filtering for each defined frequency band, then applies the Hilbert transform to obtain complex-valued analytic signals.

Parameters

data : np.ndarray EEG data from two participants with shape (2, n_epochs, n_channels, n_times)

int

Sampling rate of the EEG data in Hz

dict

Dictionary defining frequency bands: {'band_name': [fmin, fmax], ...} e.g., {'alpha': [8, 12], 'beta': [12, 20]}

bool, optional

Whether to apply bandpass filtering (default=True) If False, the Hilbert transform is applied directly to the input signals

**filter_options Additional arguments for mne.filter.filter_data, such as: - filter_length: Length of the FIR filter in samples or 's' for seconds - l_trans_bandwidth: Width of the transition band at the low cut-off frequency - h_trans_bandwidth: Width of the transition band at the high cut-off frequency

Returns

complex_signal : np.ndarray Complex-valued analytic signals with shape: (2, n_epochs, n_channels, n_freq_bands, n_times) where n_freq_bands corresponds to the number of frequency bands in freq_bands

Notes

Unlike compute_single_freq which uses the multitaper method, this function: 1. Applies a bandpass filter to isolate the specified frequency band 2. Uses the Hilbert transform to compute the analytic signal

This approach may be computationally more efficient when analyzing broader frequency bands rather than individual frequencies.

The order of frequency bands in the output corresponds to their order in the freq_bands dictionary.

Examples

bands = {'theta': [4, 8], 'alpha': [8, 13], 'beta': [13, 30]} complex_signals = compute_freq_bands(data, 256, bands, ... filter_length='auto', l_trans_bandwidth=1)

compute_nmPLV(data, sampling_rate, freq_range1, freq_range2, **filter_options)

Computes n:m Phase-Locking Value for cross-frequency coupling between participants.

This function calculates the phase synchronization between two participants when their brain oscillations operate at different frequency ranges with an integer ratio relationship (n:m).

Parameters

data : np.ndarray EEG data from two participants with shape (2, n_epochs, n_channels, n_times)

int

Sampling rate of the EEG data in Hz

List[float]

Frequency range [fmin, fmax] for participant 1

List[float]

Frequency range [fmin, fmax] for participant 2

**filter_options Additional arguments for the underlying filtering functions

Returns

con : np.ndarray n:m PLV connectivity matrix with shape (n_freq, 2n_channels, 2n_channels)

Notes

n:m Phase Locking Value measures synchronization between oscillations at different frequencies when there exists an n:m ratio between them. For example, if participant 1 shows activity at 10 Hz and participant 2 at 20 Hz, this would be a 1:2 relationship.

The function: 1. Computes the ratio between the mean frequencies of the two ranges 2. Applies appropriate phase multiplication to account for the frequency ratio 3. Calculates the phase locking value between the adjusted phases

This measure is useful for studying cross-frequency coupling between brains, where different participants might exhibit coordination at different frequency bands.

References

Palva, J. M., Palva, S., & Kaila, K. (2005). Phase synchrony among neuronal oscillations in the human cortex. Journal of Neuroscience, 25(15), 3962-3972.

Examples

Computing 1:2 PLV between alpha (participant 1) and beta (participant 2)

nm_plv = compute_nmPLV(data, sampling_rate=256, ... freq_range1=[8, 12], freq_range2=[16, 24])

compute_single_freq(data, sampling_rate, freq_range)

Computes analytic signals for each frequency in the specified range using multitaper method.

This function calculates complex-valued time-frequency representations for each integer frequency in the given range.

Parameters

data : np.ndarray EEG data from two participants with shape (2, n_epochs, n_channels, n_times)

int

Sampling rate of the EEG data in Hz

List[float]

A list [fmin, fmax] specifying the frequency range to analyze. Every integer frequency from fmin to fmax will be included.

Returns

complex_signal : np.ndarray Complex-valued analytic signals with shape: (2, n_epochs, n_channels, n_tapers, n_frequencies, n_times)

Notes

This function uses MNE's implementation of the multitaper method with 4 cycles per frequency, which provides a good balance between time and frequency resolution.

Time-frequency decomposition is calculated for each participant separately.

Examples

complex_tf = compute_single_freq(data, sampling_rate=256, freq_range=[8, 30]) print(f"Shape: {complex_tf.shape}") # (2, n_epochs, n_channels, n_tapers, 23, n_times)

compute_sync(complex_signal, mode, epochs_average=True)

Computes frequency-domain connectivity measures from analytic signals.

This function calculates various connectivity metrics between all possible channel pairs based on the input complex-valued signals.

Parameters

complex_signal : np.ndarray Complex analytic signals with shape (2, n_epochs, n_channels, n_freq_bins, n_times)

str

Connectivity measure to compute. Options: - 'envelope_corr': envelope correlation - correlation between signal envelopes - 'pow_corr': power correlation - correlation between signal power - 'plv': phase locking value - consistency of phase differences - 'ccorr': circular correlation coefficient - circular statistic for phase coupling - 'coh': coherence - normalized cross-spectrum - 'imaginary_coh': imaginary coherence - imaginary part of coherence (volume conduction resistant) - 'pli': phase lag index - asymmetry of phase difference distribution - 'wpli': weighted phase lag index - weighted version of PLI with improved properties

bool, optional

If True, connectivity values are averaged across epochs (default) If False, epoch-by-epoch connectivity is preserved

Returns

con : np.ndarray Connectivity matrix with shape: - If epochs_average=True: (n_freq, 2n_channels, 2n_channels) - If epochs_average=False: (n_freq, n_epochs, 2n_channels, 2n_channels)

Notes

Mathematical formulations for each connectivity measure:

• PLV: |⟨e^(i(φₓ-φᵧ))⟩| Measures consistency of phase differences across time

• Envelope correlation: corr(env(x), env(y)) Pearson correlation between signal envelopes

• Coherence: |⟨XY*⟩|²/(⟨|X|²⟩⟨|Y|²⟩) Normalized cross-spectrum

• Imaginary coherence: |Im(⟨XY*⟩)|/√(⟨|X|²⟩⟨|Y|²⟩) Takes only imaginary part which is less affected by volume conduction

• PLI: |⟨sign(Im(XY*))⟩| Quantifies asymmetry in phase difference distribution

• wPLI: |⟨|Im(XY)|sign(Im(XY))⟩|/⟨|Im(XY*)|⟩ Weighted version that downweights phase differences near 0 or π

Raises

ValueError If an unsupported connectivity metric is specified

indices_connectivity_interbrain(epoch_hyper)

Computes indices for interbrain connectivity analysis between EEG channels of two participants.

This function generates all possible pairs of channel indices where the first channel belongs to participant 1 and the second to participant 2, based on a merged epochs object.

Parameters

epoch_hyper : mne.Epochs A merged Epochs object containing data from both participants, with channels ordered such that the first half belongs to participant 1 and the second half to participant 2.

Returns

channels : List[Tuple[int, int]] List of tuples, each containing a pair of channel indices (i, j) where: - i is an index of a channel from participant 1 (in range 0 to n_channels/2 - 1) - j is an index of a channel from participant 2 (in range n_channels/2 to n_channels - 1)

Notes

This function assumes that the channels in epoch_hyper are organized as: [participant1_ch1, ..., participant1_chN, participant2_ch1, ..., participant2_chN]

EOG channels are automatically excluded before generating the pairs.

Examples

interbrain_pairs = indices_connectivity_interbrain(merged_epochs) print(f"Number of interbrain channel pairs: {len(interbrain_pairs)}")

indices_connectivity_intrabrain(epochs)

Computes indices for intrabrain connectivity analysis between all EEG channels.

This function generates all possible pairs of EEG channel indices for a single participant, excluding EOG channels.

Parameters

epochs : mne.Epochs A participant's Epochs object containing channel information

Returns

channels : List[Tuple[int, int]] List of tuples, each containing a pair of channel indices (i, j) where i < j. These indices correspond to all possible pairs of EEG channels for the participant.

Notes

This function automatically removes EOG channels before generating the pairs. The resulting indices can be used as input for connectivity analyses within a single brain.

Examples

channel_pairs = indices_connectivity_intrabrain(participant_epochs) print(f"Number of channel pairs: {len(channel_pairs)}")

pair_connectivity(data, sampling_rate, frequencies, mode, epochs_average=True)

Computes frequency-domain connectivity measures between two participants.

This function is a high-level interface that processes EEG data, computes analytic signals, and calculates connectivity metrics between all possible channel pairs.

Parameters

data : Union[list, np.ndarray] EEG data from two participants with shape (2, n_epochs, n_channels, n_times)

int

Sampling rate of the EEG data in Hz

Union[dict, list]

Specification of frequency bands of interest: - If dict: {'band_name': [fmin, fmax], ...} defining frequency bands e.g., {'alpha': [8, 12], 'beta': [12, 20]} - If list: [fmin, fmax] for each integer frequency in the range e.g., [5, 30] for all integer frequencies from 5 to 30 Hz

str

Connectivity measure to compute. Options: - 'envelope_corr': envelope correlation - 'pow_corr': power correlation - 'plv': phase locking value - 'ccorr': circular correlation coefficient - 'coh': coherence - 'imaginary_coh': imaginary coherence - 'pli': phase lag index - 'wpli': weighted phase lag index

bool, optional

If True, connectivity values are averaged across epochs (default) If False, epoch-by-epoch connectivity is preserved

Returns

result : np.ndarray Connectivity matrix with shape: - If epochs_average=True: (n_freq, 2n_channels, 2n_channels) - If epochs_average=False: (n_freq, n_epochs, 2n_channels, 2n_channels)

The channels are ordered as:
[participant1_ch1, ..., participant1_chN, participant2_ch1, ..., participant2_chN]

To extract only interbrain connectivity values, slice the last two dimensions:
result[:, :, 0:n_channels, n_channels:2*n_channels]

Notes

This function handles the complete process from raw EEG data to connectivity metrics: 1. Computes analytic signals for the specified frequencies 2. Applies the selected connectivity measure

Each connectivity metric has different mathematical properties and interpretations: - PLV measures consistency in phase differences regardless of amplitude - Coherence measures linear relationship between signals in frequency domain - Imaginary coherence is less susceptible to volume conduction - PLI and wPLI are robust against common sources and amplitude effects

Time complexity scales with O(n_channels^2 × n_epochs × n_frequencies × n_times)

Examples

Computing alpha band PLV

alpha_plv = pair_connectivity( ... [subj1_data, subj2_data], ... sampling_rate=256, ... frequencies={'alpha': [8, 12]}, ... mode='plv' ... )

pow(epochs, fmin, fmax, n_fft, n_per_seg, epochs_average)

Computes the Power Spectral Density (PSD) on Epochs using Welch's method.

This function calculates the power spectrum for each channel across the specified frequency range. EOG channels are automatically dropped before computation.

Parameters

epochs : mne.Epochs A participant's Epochs object containing EEG data of shape (n_epochs, n_channels, n_times).

float

Minimum frequency in Hz to include in PSD calculation.

float

Maximum frequency in Hz to include in PSD calculation.

int

Length of FFT used. Must be >= n_per_seg. If larger, the segments will be zero-padded. If n_per_seg is None, n_fft must be <= number of time points.

int or None

Length of each Welch segment (windowed with a Hamming window). If None, n_per_seg is set equal to n_fft.

bool

If True, PSD values are averaged over epochs. If False, PSD won't be averaged (the time course is maintained).

Returns

psd_tuple : namedtuple A named tuple containing: - freq_list: ndarray of frequencies (frequency bins) in Hz - psd: ndarray of PSD values with shape: - If epochs_average=True: (n_channels, n_frequencies) - If epochs_average=False: (n_epochs, n_channels, n_frequencies)

Notes

This function can be iterated on groups and/or conditions: (for epochs in epochs['epochs_%s_%s_%s' % (subj, group, cond_name)]).

The PSD values are computed using Welch's method with Hamming windows and expressed in µV²/Hz.

Examples

from mne import Epochs result = pow(epochs, fmin=1, fmax=40, n_fft=512, n_per_seg=256, epochs_average=True) frequencies = result.freq_list power_values = result.psd

xwt(sig1, sig2, freqs, n_cycles=5.0, mode='xwt')

Performs cross-wavelet transform or wavelet coherence analysis between two signals.

This function computes time-frequency representations of two signals and their cross-spectrum using continuous wavelet transform.

Parameters

sig1 : mne.Epochs EEG data of first participant

mne.Epochs

EEG data of second participant

Union[int, np.ndarray]

Frequencies of interest in Hz

float, optional

Number of cycles in the Morlet wavelet (default=5.0) Controls the time-frequency resolution trade-off

str, optional

Type of analysis to perform (default="xwt"). Options: - 'power': absolute value of cross-wavelet transform - 'phase': phase angles of cross-wavelet transform - 'xwt': raw cross-wavelet transform (complex-valued) - 'wtc': wavelet coherence

Returns

data : np.ndarray Result of the wavelet analysis with shape: (n_chans1, n_chans2, n_epochs, n_freqs, n_samples)

Notes

This function provides different types of wavelet-based analyses:

• Cross-wavelet transform (XWT): Reveals common power and relative phase between two time series in time-frequency space

• Wavelet coherence (WTC): Measures the correlation between two signals in the time-frequency domain normalized between 0 and 1

Unlike Fourier-based methods, wavelet analysis maintains time resolution, allowing for the study of non-stationary signals and transient relationships.

The function automatically checks that both signals have the same sampling rate, number of epochs, channels, and samples.

References

Grinsted, A., Moore, J. C., & Jevrejeva, S. (2004). Application of the cross wavelet transform and wavelet coherence to geophysical time series. Nonlinear processes in geophysics, 11(5/6), 561-566.

Maraun, D., & Kurths, J. (2004). Cross wavelet analysis: significance testing and pitfalls. Nonlinear Processes in Geophysics, 11(4), 505-514.

Examples

Computing wavelet coherence between two EEG signals

coherence = xwt(subj1_epochs, subj2_epochs, freqs=np.arange(4, 40, 1), ... n_cycles=7, mode='wtc')