diff --git a/CITATION.cff b/CITATION.cff index 96be26d1..b6dc8467 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -20,8 +20,15 @@ authors: - given-names: "Alexandre" family-names: "Gramfort" orcid: https://orcid.org/0000-0001-9791-4404 + - given-names: "Thomas Samuel" + family-names: "Binns" + orcid: https://orcid.org/0000-0003-0657-0891 + - given-names: "Mohammad" + family-names: "Orabe" + orcid: https://orcid.org/0009-0004-7177-799X + title: "mne-connectivity" -version: 0.2.0 -date-released: 2022-01-13 +version: 0.6.0 +date-released: 2024-12-06 url: "https://github.com/mne-tools/mne-connectivity" diff --git a/doc/authors.inc b/doc/authors.inc index 21a3fad4..61ead35b 100644 --- a/doc/authors.inc +++ b/doc/authors.inc @@ -14,3 +14,5 @@ .. _Thomas Binns: https://github.com/tsbinns .. _Tien Nguyen: https://github.com/nguyen-td .. _Richard Köhler: https://github.com/richardkoehler +.. _Mohammad Orabe: https://github.com/orabe +.. _Mina Jamshidi: https://github.com/minajamshidi diff --git a/doc/references.bib b/doc/references.bib index 14271fd5..473e605a 100644 --- a/doc/references.bib +++ b/doc/references.bib @@ -82,6 +82,16 @@ @article{HaufeEtAl2013 year = {2013} } +@article{HaufeEtAl2014, + author={Haufe, Stefan and Meinecke, Frank and G{\"o}rgen, Kai and D{\"a}hne, Sven and Haynes, John-Dylan and Blankertz, Benjamin and Bie{\ss}mann, Felix}, + doi = {10.1016/j.neuroimage.2013.10.067}, + journal={NeuroImage}, + pages={96--110}, + title={On the interpretation of weight vectors of linear models in multivariate neuroimaging}, + volume={87}, + year={2014}, +} + @article{HippEtAl2012, author = {Hipp, Joerg F and Hawellek, David J and Corbetta, Maurizio and Siegel, Markus and Engel, Andreas K}, doi = {10.1038/nn.3101}, @@ -268,6 +278,16 @@ @article{VinckEtAl2015 year = {2015} } +@article{ViriyopaseEtAl2012, + author={Viriyopase, Atthaphon and Bojak, Ingo and Zeitler, Magteld and Gielen, Stan}, + doi={10.3389/fncom.2012.00049}, + journal={Frontiers in Computational Neuroscience}, + pages={49}, + title={When long-range zero-lag synchronization is feasible in cortical networks}, + volume={6}, + year={2012} +} + @article{Whittle1963, author = {Whittle, Peter}, doi = {10.1093/biomet/50.1-2.129}, diff --git a/doc/whats_new.rst b/doc/whats_new.rst index 0199231e..ada7d447 100644 --- a/doc/whats_new.rst +++ b/doc/whats_new.rst @@ -24,7 +24,7 @@ Version 0.7 (in dev) Enhancements ~~~~~~~~~~~~ -- +- Add support for a new multivariate connectivity method (canonical coherence; ``cacoh``) in :func:`mne_connectivity.spectral_connectivity_epochs` and :func:`mne_connectivity.spectral_connectivity_time` by `Thomas Binns`_ and `Mohammad Orabe`_ and `Mina Jamshidi`_ (:pr:`163`). Bug ~~~ @@ -38,6 +38,9 @@ API Authors ~~~~~~~ +* `Thomas Binns`_ +* `Mohammad Orabe`_ +* `Mina Jamshidi`_ :doc:`Find out what was new in previous releases ` diff --git a/examples/cacoh.py b/examples/cacoh.py new file mode 100644 index 00000000..9f3d2838 --- /dev/null +++ b/examples/cacoh.py @@ -0,0 +1,433 @@ +""" +======================================== +Compute multivariate coherency/coherence +======================================== + +This example demonstrates how canonical coherency (CaCoh) +:footcite:`VidaurreEtAl2019` - a multivariate method based on coherency - can +be used to compute connectivity between whole sets of sensors, alongside +spatial patterns of the connectivity. +""" + +# Authors: Thomas S. Binns +# Mohammad Orabe +# License: BSD (3-clause) + +# %% +import numpy as np +from matplotlib import pyplot as plt + +from mne_connectivity import ( + make_signals_in_freq_bands, + seed_target_indices, + spectral_connectivity_epochs, +) + +############################################################################### +# Background +# ---------- +# +# Multivariate forms of signal analysis allow you to simultaneously consider +# the activity of multiple signals. In the case of connectivity, the +# interaction between multiple sensors can be analysed at once, producing a +# single connectivity spectrum. This approach brings not only practical +# benefits (e.g. easier interpretability of results from the dimensionality +# reduction), but can also offer methodological improvements (e.g. enhanced +# signal-to-noise ratio). +# +# A popular bivariate measure of connectivity is coherency/coherence, which +# looks at the correlation between two signals in the frequency domain. +# However, in cases where interactions between multiple signals are of +# interest, computing connectivity between all possible combinations of signals +# leads to a very large number of results which is difficult to interpret. A +# common approach is to average results across these connections, however this +# risks reducing the signal-to-noise ratio of results and burying interactions +# that are present between only a small number of channels. +# +# Canonical coherency (CaCoh) is a multivariate form of coherency that uses +# eigendecomposition-derived spatial filters to extract the underlying +# components of connectivity in a frequency-resolved manner +# :footcite:`VidaurreEtAl2019`. This approach goes beyond simply aggregating +# information across all possible combinations of signals. +# +# It is similar to multivariate methods based on the imaginary part of +# coherency (MIC & MIM :footcite:`EwaldEtAl2012`; see :doc:`mic_mim` and +# :doc:`compare_coherency_methods`). + + +############################################################################### +# Data Simulation +# --------------- +# +# To demonstrate the CaCoh method, we will use some simulated data consisting +# of two sets of interactions between signals in a given frequency range: +# +# - 5 seeds and 3 targets interacting in the 10-12 Hz frequency range. +# - 5 seeds and 3 targets interacting in the 23-25 Hz frequency range. +# +# We can consider the seeds and targets to be signals of different modalities, +# e.g. cortical EEG signals and subcortical LFP signals, cortical EEG signals +# and muscular EMG signals, etc.... We use the +# :func:`~mne_connectivity.make_signals_in_freq_bands` function to simulate +# these signals. + +# %% + +# Generate simulated data +data_10_12 = make_signals_in_freq_bands( + n_seeds=5, + n_targets=3, + freq_band=(10, 12), # 10-12 Hz interaction + rng_seed=42, +) + +data_23_25 = make_signals_in_freq_bands( + n_seeds=5, + n_targets=3, + freq_band=(23, 25), # 23-25 Hz interaction + rng_seed=44, +) + +# Combine data into a single object +data = data_10_12.add_channels([data_23_25]) + +############################################################################### +# Computing CaCoh +# --------------- +# +# Having simulated the signals, we can create the indices for computing +# connectivity between all seeds and all targets in a single multivariate +# connection (see :doc:`handling_ragged_arrays` for more information), after +# which we compute connectivity. +# +# For CaCoh, a set of spatial filters are found that will maximise the +# estimated connectivity between the seed and target signals. These maximising +# filters correspond to the eigenvectors with the largest eigenvalue, derived +# from an eigendecomposition of information from the cross-spectral density +# (Eq. 8 of :footcite:`VidaurreEtAl2019`): +# +# :math:`\textrm{CaCoh}=\Large{\frac{\boldsymbol{a}^T\boldsymbol{D}(\Phi) +# \boldsymbol{b}}{\sqrt{\boldsymbol{a}^T\boldsymbol{a}\boldsymbol{b}^T +# \boldsymbol{b}}}}` +# +# where: :math:`\boldsymbol{D}(\Phi)` is the cross-spectral density between +# seeds and targets transformed for a given phase angle :math:`\Phi`; and +# :math:`\boldsymbol{a}` and :math:`\boldsymbol{b}` are eigenvectors for the +# seeds and targets, such that :math:`\boldsymbol{a}^T\boldsymbol{D}(\Phi) +# \boldsymbol{b}` maximises coherency between the seeds and targets. All +# elements are frequency-dependent, however this is omitted for readability. +# +# CaCoh is complex-valued in the range :math:`[-1, 1]` where the sign reflects +# the phase angle of the interaction (akin to coherency). Taking the absolute +# value is akin to taking the coherence, which is the magnitude of the +# interaction regardless of phase angle. + +# %% + +# Generate connectivity indices +seeds = [0, 1, 2, 3, 4, 8, 9, 10, 11, 12] +targets = [5, 6, 7, 13, 14, 15] +multivar_indices = ([seeds], [targets]) + +# Compute CaCoh +cacoh = spectral_connectivity_epochs( + data, method="cacoh", indices=multivar_indices, sfreq=100, fmin=3, fmax=35 +) +print(f"Results shape: {cacoh.get_data().shape} (connections x frequencies)") + +# Get absolute CaCoh +cacoh_abs = np.abs(cacoh.get_data())[0] + +############################################################################### +# As you can see below, using CaCoh we have summarised the most relevant +# connectivity information from our 10 seed channels and 6 target channels as a +# single spectrum of connectivity values. This lower-dimensional representation +# of signal interactions is much more interpretable when analysing connectivity +# in complex systems such as the brain. + +# %% + +# Plot CaCoh +fig, axis = plt.subplots(1, 1) +axis.plot(cacoh.freqs, cacoh_abs, linewidth=2) +axis.set_xlabel("Frequency (Hz)") +axis.set_ylabel("Connectivity (A.U.)") +fig.suptitle("CaCoh") + +############################################################################### +# Note that we plot the absolute values of the results (coherence) rather than +# the complex values (coherency). The absolute value of connectivity will +# generally be of most interest. However, information such as the phase of +# interaction can only be extracted from the complex-valued results, e.g. with +# the :func:`numpy.angle` function. + +# %% + +# Plot phase of connectivity +fig, axis = plt.subplots(1, 1) +axis.plot(cacoh.freqs, np.angle(cacoh.get_data()[0]), linewidth=2) +axis.set_xlabel("Frequency (Hz)") +axis.set_ylabel("Phase of connectivity (radians)") +fig.suptitle("CaCoh") + +############################################################################### +# CaCoh versus coherence +# ---------------------- +# +# To further demonstrate the signal-to-noise ratio benefits of CaCoh, below we +# compute connectivity between each seed and target using bivariate coherence. +# With our 10 seeds and 6 targets, this gives us a total of 60 unique +# connections which is very difficult to interpret without aggregating some +# information. A common approach is to simply average across these connections, +# which we do below. + +# %% + +# Define bivariate connectivity indices +bivar_indices = seed_target_indices(seeds, targets) + +# Compute bivariate coherence +coh = spectral_connectivity_epochs( + data, method="coh", indices=bivar_indices, sfreq=100, fmin=3, fmax=35 +) +print(f"Original results shape: {coh.get_data().shape} (connections x frequencies)") + +# Average results across connections +coh_mean = np.mean(coh.get_data(), axis=0) +print(f"Averaged results shape: {coh_mean.shape} (connections x frequencies)") + +############################################################################### +# Plotting the bivariate and multivariate results together, we can see that +# coherence still captures the interactions at 10-12 Hz and 23-25 Hz, however +# the scale of the connectivity is much smaller. This reflects the fact that +# CaCoh is able to capture the relevant components of interactions between +# multiple signals, regardless of whether they are present in all channels. + +# %% + +# Plot CaCoh & Coh +fig, axis = plt.subplots(1, 1) +axis.plot(cacoh.freqs, cacoh_abs - np.min(cacoh_abs), linewidth=2, label="CaCoh") +axis.plot(coh.freqs, coh_mean - np.min(coh_mean), linewidth=2, label="Coh") +axis.set_xlabel("Frequency (Hz)") +axis.set_ylabel("Baseline-corrected connectivity (A.U.)") +axis.legend() +fig.suptitle("CaCoh vs. coherence") + +############################################################################### +# The ability of multivariate connectivity methods to capture the underlying +# components of connectivity is extremely useful when dealing with data from +# a large number of channels, with inter-channel interactions at distinct +# frequencies, a problem explored in more detail in the :doc:`mic_mim` example. + +############################################################################### +# Extracting spatial information from CaCoh +# ----------------------------------------- +# +# Whilst a lower-dimensional representation of connectivity information is +# useful, we lose information about which channels are involved in the +# connectivity. Thankfully, this information can be recovered by constructing +# spatial patterns of connectivity from the spatial filters +# :footcite:`HaufeEtAl2014`. +# +# The spatial patterns are stored under ``attrs['patterns']`` of the +# connectivity class, with one value per frequency for each channel in the +# seeds and targets. The patterns can be positive- and negative-valued. Sign +# differences of the patterns can be used to visualise the orientation of +# underlying dipole sources, whereas their absolute value reflects the strength +# of a channel's contribution to the connectivity component. The spatial +# patterns are **not** bound between :math:`[-1, 1]`. +# +# Averaging across the patterns in the 10-12 Hz and 23-25 Hz ranges, we can see +# how it is possible to identify which channels are contributing to +# connectivity at different frequencies. + +# %% + +freqs = cacoh.freqs +fbands = ((10, 12), ((23, 25))) + +fig, axes = plt.subplots(1, 2) + +# patterns have shape [seeds/targets x cons x channels x freqs (x times)] +patterns = np.abs(np.array(cacoh.attrs["patterns"])) +seed_pattern = patterns[0, :, : len(seeds)] +target_pattern = patterns[1, :, : len(targets)] + +vmin = np.nanmin(patterns) +vmax = np.nanmax(patterns) + +for axis, fband in zip(axes, fbands): + # average across frequencies + seed_pattern_fband = np.mean( + seed_pattern[0, :, freqs.index(fband[0]) : freqs.index(fband[1]) + 1], axis=1 + ) + target_pattern_fband = np.mean( + target_pattern[0, :, freqs.index(fband[0]) : freqs.index(fband[1]) + 1], axis=1 + ) + + # combine into a single array + pattern_fband = np.concatenate((seed_pattern_fband, target_pattern_fband), axis=0) + + # plot the pattern + mesh = axis.pcolormesh( + np.flip(np.expand_dims(pattern_fband, 1)), vmin=vmin, vmax=vmax + ) + axis.set_yticks([1.5, 4.5, 8.5, 13.5]) + axis.set_xticks([0.5]) + axis.set_xticklabels([f"{fband[0]}-{fband[1]}"]) + +# Label axes +fig.suptitle("Spatial patterns") +axes[0].set_yticklabels( + [ + "Targets\n(23-25 Hz)", + "Targets\n(10-12 Hz)", + "Seeds\n(23-25 Hz)", + "Seeds\n(10-12 Hz)", + ], + rotation=45, + va="center", +) +axes[0].set_ylabel("Channels") +axes[1].get_yaxis().set_visible(False) +fig.text(0.47, 0.02, "Frequency band (Hz)", ha="center") + +# Set colourbar +fig.subplots_adjust(right=0.8) +cbar_axis = fig.add_axes([0.85, 0.15, 0.02, 0.7]) +fig.colorbar(mesh, cax=cbar_axis) +cbar_axis.set_ylabel("Contribution to connectivity (A.U.)") +cbar_axis.set_yticks([vmin, vmax]) +cbar_axis.set_yticklabels(["Low", "High"]) + +plt.show() + +############################################################################### +# For an example on interpreting spatial filters with real data, see the +# :doc:`mic_mim` example. + +############################################################################### +# Handling high-dimensional data +# ------------------------------ +# +# An important issue to consider when using these multivariate methods is +# overfitting, which risks biasing connectivity estimates to maximise noise in +# the data. This risk can be reduced by performing a preliminary dimensionality +# reduction prior to estimating the connectivity with a singular value +# decomposition (Eq. 15 of :footcite:`VidaurreEtAl2019`). The degree of this +# dimensionality reduction can be specified using the ``rank`` argument, which +# by default will not perform any dimensionality reduction (assuming your data +# is full rank; see below if not). Choosing an expected rank of the data +# requires *a priori* knowledge about the number of components you expect to +# observe in the data. +# +# When comparing CaCoh scores across recordings, **it is highly recommended +# to estimate connectivity from the same number of channels (or equally from +# the same degree of rank subspace projection)** to avoid biases in +# connectivity estimates. Bias can be avoided by specifying a consistent rank +# subspace to project to using the ``rank`` argument, standardising your +# connectivity estimates regardless of changes in e.g. the number of channels +# across recordings. Note that this does not refer to the number of seeds and +# targets *within* a connection being identical, rather to the number of seeds +# and targets *across* connections. +# +# Here, we project our seed and target data to only the first 2 components of +# our rank subspace. Results show that the general spectral pattern of +# connectivity is retained in the rank subspace-projected data, suggesting that +# a fair degree of redundant connectivity information is contained in the +# excluded components of the seed and target data. +# +# We also assert that the spatial patterns of MIC are returned in the original +# sensor space despite this rank subspace projection, being reconstructed using +# the products of the singular value decomposition (Eqs. 46 & 47 of +# :footcite:`EwaldEtAl2012`). + +# %% + +# Compute CaCoh following rank subspace projection +cacoh_red = spectral_connectivity_epochs( + data, + method="cacoh", + indices=multivar_indices, + sfreq=100, + fmin=3, + fmax=35, + rank=([2], [2]), +) + +# compare standard and rank subspace-projected CaCoh +fig, axis = plt.subplots(1, 1) +axis.plot(cacoh.freqs, np.abs(cacoh.get_data()[0]), linewidth=2, label="standard CaCoh") +axis.plot( + cacoh_red.freqs, + np.abs(cacoh_red.get_data()[0]), + linewidth=2, + label="rank subspace (2) CaCoh", +) +axis.set_xlabel("Frequency (Hz)") +axis.set_ylabel("Connectivity (A.U.)") +axis.legend(loc="lower right") +fig.suptitle("CaCoh") + +# no. channels equal with and without projecting to rank subspace for patterns +assert patterns[0, 0].shape[0] == np.array(cacoh_red.attrs["patterns"])[0, 0].shape[0] +assert patterns[1, 0].shape[0] == np.array(cacoh_red.attrs["patterns"])[1, 0].shape[0] + +############################################################################### +# See :doc:`mic_mim` for an example of applying the rank subspace +# projection to real data with a large number of channels. +# +# In the case that your data is not full rank and ``rank`` is left as ``None``, +# an automatic rank computation is performed and an appropriate degree of +# dimensionality reduction will be enforced. The rank of the data is determined +# by computing the singular values of the data and finding those within a +# factor of :math:`1e^{-6}` relative to the largest singular value. +# +# Whilst unlikely, there may be scenarios in which this threshold may be too +# lenient. In these cases, you should inspect the singular values of your data +# to identify an appropriate degree of dimensionality reduction to perform, +# which you can then specify manually using the ``rank`` argument. The code +# below shows one possible approach for finding an appropriate rank of +# close-to-singular data with a more conservative threshold. + +# %% + +# gets the singular values of the data across epochs +s = np.linalg.svd(data, compute_uv=False).min(axis=0) +# finds how many singular values are 'close' to the largest singular value +rank = np.count_nonzero(s >= s[0] * 1e-4) # 1e-4 is the 'closeness' criteria, which is +# a hyper-parameter + +############################################################################### +# Limitations +# ----------- +# +# Multivariate methods offer many benefits in the form of dimensionality +# reduction and signal-to-noise ratio improvements. However, no method is +# perfect. When we simulated the data, we mentioned how we considered the seeds +# and targets to be signals of different modalities. This is an important +# factor in whether CaCoh should be used over methods based solely on the +# imaginary part of coherency such as MIC and MIM. +# +# In short, if you want to examine connectivity between signals from the same +# modality, you should consider using another method instead of CaCoh. Rather, +# methods based on the imaginary part of coherency such as MIC and MIM should +# be used to avoid spurious connectivity estimates stemming from e.g. volume +# conduction artefacts. +# +# On the other hand, if you want to examine connectivity between signals from +# different modalities, CaCoh is a more appropriate method than MIC/MIM. This +# is because volume conduction artefacts are of less concern, and CaCoh does +# not risk biasing connectivity estimates towards interactions with particular +# phase lags like MIC/MIM. +# +# These scenarios are described in more detail in the +# :doc:`compare_coherency_methods` example. + +############################################################################### +# References +# ---------- +# .. footbibliography:: + +# %% diff --git a/examples/compare_coherency_methods.py b/examples/compare_coherency_methods.py new file mode 100644 index 00000000..70848b6b --- /dev/null +++ b/examples/compare_coherency_methods.py @@ -0,0 +1,530 @@ +""" +===================================== +Comparison of coherency-based methods +===================================== + +This example demonstrates the distinct forms of information captured by +coherency-based connectivity methods, and highlights the scenarios in which +these different methods should be applied. +""" + +# Authors: Thomas S. Binns +# Mohammad Orabe +# License: BSD (3-clause) + +# %% +import numpy as np +from matplotlib import pyplot as plt + +from mne_connectivity import ( + make_signals_in_freq_bands, + seed_target_indices, + spectral_connectivity_epochs, +) + +############################################################################### +# An introduction to coherency-based connectivity methods +# ------------------------------------------------------- +# +# MNE-Connectivity supports several methods based on coherency. These are: +# +# - coherency (Cohy) +# - coherence (Coh; absolute coherency) +# - imaginary part of coherency (ImCoh) +# - canonical coherency (CaCoh) +# - maximised imaginary part of coherency (MIC) +# - multivariate interaction measure (MIM) +# +# | +# +# All of these methods centre on Cohy, a complex-valued estimate of the +# correlation between signals in the frequency domain. It is an undirected +# measure of connectivity, being invariant to the direction of information flow +# between signals. +# +# The common approach for handling these complex-valued coherency scores is to +# either take their absolute values (Coh) or their imaginary values (ImCoh +# :footcite:`NolteEtAl2004`). +# +# In addition to these traditional bivariate connectivity measures (i.e. +# between two signals), advanced multivariate measures (i.e. between groups of +# signals) have also been developed based on Cohy (CaCoh +# :footcite:`VidaurreEtAl2019`; can take the absolute value for a multivariate +# form of Coh; see :doc:`cacoh`) or ImCoh (MIC & MIM :footcite:`EwaldEtAl2012`; +# see :doc:`mic_mim`). +# +# Despite their similarities, there are distinct scenarios in which these +# different methods are most appropriate, as we will show in this example. + +############################################################################### +# Zero and non-zero time-lag interactions +# --------------------------------------- +# +# The key difference between Cohy/Coh and ImCoh is how information about zero +# time-lag interactions is captured. +# +# We generally assume that communication within the brain involves some delay +# in the flow of information (i.e. a non-zero time-lag). This reflects the time +# taken for: the propagation of action potentials along axons; the release of +# neurotransmitters from presynaptic terminals and binding to receptors on +# postsynaptic terminals; etc... +# +# In contrast, interactions with no delay (i.e. a zero time-lag) are often +# considered to reflect non-physiological activity, such as volume conduction +# - the propagation of electrical activity through the brain's conductive +# tissue from a single source to multiple electrodes simultaneously +# :footcite:`NolteEtAl2004`. Such interactions therefore do not reflect +# genuine, physiological communication between brain regions. Naturally, +# having a method that can discard spurious zero time-lag connectivity +# estimates is very desirable. +# +# **Note:** Not all zero time-lag interactions are necessarily +# non-physiological :footcite:`ViriyopaseEtAl2012`. +# +# To demonstrate the differences in how Cohy/Coh and ImCoh handle zero time-lag +# interactions, we simulate two sets of data with: +# +# 1. A non-zero time-lag interaction at 10-12 Hz. +# 2. A zero time-lag interaction at 23-25 Hz. + +# %% + +# Generate simulated data +data_delay = make_signals_in_freq_bands( + n_seeds=3, + n_targets=3, + freq_band=(10, 12), # 10-12 Hz interaction + connection_delay=2, # samples; non-zero time-lag + rng_seed=42, +) + +data_no_delay = make_signals_in_freq_bands( + n_seeds=3, + n_targets=3, + freq_band=(23, 25), # 23-25 Hz interaction + connection_delay=0, # samples; zero time-lag + rng_seed=44, +) + +# Combine data into a single object +data = data_delay.add_channels([data_no_delay]) + +############################################################################### +# We compute the connectivity of these simulated signals using CaCoh (a +# multivariate form of Cohy/Coh) and MIC (a multivariate form of ImCoh). + +# %% + +# Generate connectivity indices +seeds = [0, 1, 2, 6, 7, 8] +targets = [3, 4, 5, 9, 10, 11] +bivar_indices = seed_target_indices(seeds, targets) +multivar_indices = ([seeds], [targets]) + +# Compute CaCoh & MIC +(cacoh, mic) = spectral_connectivity_epochs( + data, method=["cacoh", "mic"], indices=multivar_indices, sfreq=100, fmin=3, fmax=35 +) + +############################################################################### +# As you can see, both CaCoh and MIC capture the non-zero time-lag interaction +# at 10-12 Hz, however only CaCoh captures the zero time-lag interaction at +# 23-25 Hz. + +# %% + +# Plot CaCoh & MIC +fig, axis = plt.subplots(1, 1) +axis.plot(cacoh.freqs, np.abs(cacoh.get_data()[0]), linewidth=2, label="CaCoh") +axis.plot( + mic.freqs, np.abs(mic.get_data()[0]), linewidth=2, label="MIC", linestyle="--" +) +axis.set_xlabel("Frequency (Hz)") +axis.set_ylabel("Connectivity (A.U.)") +axis.annotate("Non-zero\ntime-lag\ninteraction", xy=(13.5, 0.85)) +axis.annotate("Zero\ntime-lag\ninteraction", xy=(26.5, 0.85)) +axis.legend(loc="upper left") +fig.suptitle("CaCoh vs. MIC\nNon-zero & zero time-lags") + + +# %% + + +def plot_connectivity_circle(): + """Plot a circle with radius 1, real and imag. axes, and angles marked.""" + fig, axis = plt.subplots(1, 1) + t = np.linspace(0, 2 * np.pi, 100) + axis.plot(np.cos(t), np.sin(t), color="k", linewidth=0.1) + axis.plot([-1, 1], [0, 0], color="k", linestyle="--") + axis.plot([0, 0], [-1, 1], color="k", linestyle="--") + axis.axis("off") + + fontdict = {"fontsize": 10} + qpi = np.pi / 4 + axis.text(1, 0, " 0°", ha="left", va="center", fontdict=fontdict) + axis.text(np.pi / 4, np.pi / 4, "45°", ha="center", va="center", fontdict=fontdict) + axis.text(0, 1, "90°", ha="center", va="bottom", fontdict=fontdict) + axis.text(-qpi, qpi, "135°", ha="center", va="center", fontdict=fontdict) + axis.text(-1, 0, "180°", ha="right", va="center", fontdict=fontdict) + axis.text(-qpi, -qpi, "-135°", ha="center", va="center", fontdict=fontdict) + axis.text(0, -1, "-90°", ha="center", va="top", fontdict=fontdict) + axis.text(qpi, -qpi, "-45°", ha="center", va="center", fontdict=fontdict) + + fontdict = {"fontsize": 12} + axis.text(1.15, 0, " Real", ha="left", va="center", fontdict=fontdict) + axis.text(0, 1.15, "Imaginary", ha="center", va="bottom", fontdict=fontdict) + axis.text(0, 0, "0 ", ha="right", va="top", fontdict=fontdict) + axis.text(-1, 0, "-1", ha="left", va="top", fontdict=fontdict) + axis.text(1, 0, "+1", ha="right", va="top", fontdict=fontdict) + axis.text(0, -1, "-1 ", ha="right", va="bottom", fontdict=fontdict) + axis.text(0, 1, "+1 ", ha="right", va="top", fontdict=fontdict) + + axis.set_aspect("equal") + + return fig, axis + + +############################################################################### +# The different interactions (not) captured by CaCoh and MIC can be understood +# by visualising the complex values of the interactions. + +# %% + +# Get complex connectivity values at frequency bands +freqs = cacoh.freqs +cacoh_10_12 = np.mean(cacoh.get_data()[0, freqs.index(10) : freqs.index(12) + 1]) +cacoh_23_25 = np.mean(cacoh.get_data()[0, freqs.index(23) : freqs.index(25) + 1]) + +# Plot complex connectivity values +fig, axis = plot_connectivity_circle() +axis.quiver( + 0, + 0, + np.real(cacoh_10_12), + np.imag(cacoh_10_12), + units="xy", + scale=1, + linewidth=2, + color="C2", + label="10-12 Hz", +) +axis.quiver( + 0, + 0, + np.real(cacoh_23_25), + np.imag(cacoh_23_25), + units="xy", + scale=1, + linewidth=2, + color="C3", + label="23-25 Hz", + zorder=99, +) +axis.legend(loc="upper right", bbox_to_anchor=[1.1, 1.1]) + +############################################################################### +# Above, we plot the complex-valued CaCoh scores for the 10-12 Hz and 23-25 Hz +# interactions as vectors with origin :math:`(0, 0)` bound within a circle of +# radius 1 (reflecting the fact that coherency scores span the set of complex +# values in the range :math:`[-1, 1]`). +# +# The circumference of the circle spans the range :math:`(-\pi, \pi]`. The real +# axis corresponds to vectors with angles of 0° (:math:`0\pi`; positive +# values) or 180° (:math:`\pi`; negative values). The imaginary axis +# corresponds to vectors with angles of 90° (:math:`\frac{1}{2}\pi`; positive +# values) or -90° (:math:`-\frac{1}{2}\pi`; negative values). +# +# Zero time-lag interactions have angles of 0° and 180° (i.e. no phase +# difference), corresponding to a non-zero real component, but a zero-valued +# imaginary component. We see this nicely for the 23-25 Hz interaction, which +# has an angle of ~0°. Taking the absolute CaCoh value shows us the magnitude +# of this interaction to be ~0.9. However, first projecting this information to +# the imaginary axis gives us a magnitude of ~0. +# +# In contrast, non-zero time-lag interactions do not lie on the real axis (i.e. +# a phase difference), corresponding to non-zero real and imaginary components. +# We see this nicely for the 10-12 Hz interaction, which has an angle of ~-75°. +# Taking the absolute CaCoh value shows us the magnitude of this interaction to +# be ~0.9, which is also seen when first projecting this information to the +# imaginary axis. +# +# This distinction is why connectivity methods utilising information from both +# real and imaginary components (Cohy, Coh, CaCoh) capture both zero and +# non-zero time-lag interactions, whereas methods using only the imaginary +# component (ImCoh, MIC, MIM) capture only non-zero time-lag interactions. +# +# The ability to capture these different interactions is not a feature specific +# to multivariate connectivity methods, as shown below for the bivariate +# methods Coh and ImCoh. + +# %% + +# Compute Coh & ImCoh +(coh, imcoh) = spectral_connectivity_epochs( + data, method=["coh", "imcoh"], indices=bivar_indices, sfreq=100, fmin=3, fmax=35 +) + +coh_mean = np.mean(coh.get_data(), axis=0) +imcoh_mean = np.mean(np.abs(imcoh.get_data()), axis=0) + +coh_mean_subbed = coh_mean - np.mean(coh_mean) +imcoh_mean_subbed = imcoh_mean - np.mean(imcoh_mean) + +# Plot Coh & ImCoh +fig, axis = plt.subplots(1, 1) +axis.plot(coh.freqs, coh_mean_subbed, linewidth=2, label="Coh") +axis.plot(imcoh.freqs, imcoh_mean_subbed, linewidth=2, label="ImCoh", linestyle="--") +axis.set_xlabel("Frequency (Hz)") +axis.set_ylabel("Mean-corrected connectivity (A.U.)") +axis.annotate("Non-zero\ntime-lag\ninteraction", xy=(13, 0.25)) +axis.annotate("Zero\ntime-lag\ninteraction", xy=(25, 0.25)) +axis.legend(loc="upper left") +fig.suptitle("Coh vs. ImCoh\nNon-zero & zero time-lags") + +############################################################################### +# When different coherency-based methods are most appropriate +# ----------------------------------------------------------- +# +# With this information, we can define situations under which these different +# approaches are most appropriate. +# +# | +# +# **In situations where non-physiological zero time-lag interactions are +# assumed, methods based on only the imaginary part of coherency (ImCoh, MIC, +# MIM) should be used.** Examples of situations include: +# +# - Connectivity between channels of a single modality. +# - Connectivity between channels of different modalities where the same +# reference is used. +# +# Note that this applies not only to sensor-space signals, but also to +# source-space signals where remnants of these non-physiological interactions +# may remain even after source reconstruction. +# +# | +# +# **In situations where non-physiological zero time-lag interactions are not +# assumed, methods based on real and imaginary parts of coherency (Cohy, Coh, +# CaCoh) should be used.** An example includes: +# +# - Connectivity between channels of different modalities where different +# references are used. +# +# | +# +# Equally, when there are no non-physiological zero time-lag interactions, one +# should not used methods based on only the imaginary part of coherency. There +# are two key reasons: +# +# **1. Discarding physiological zero time-lag interactions** +# +# First, not all zero time-lag interactions are non-physiological +# :footcite:`ViriyopaseEtAl2012`. Accordingly, methods based on only the +# imaginary part of coherency may lead to information about genuine +# connectivity being lost. +# +# In situations where non-physiological zero time-lag interactions are present, +# the potential loss of physiological information is generally acceptable to +# avoid spurious connectivity estimates. However, unnecessarily discarding this +# information can of course be detrimental. +# +# **2. Biasing interactions based on the angle of interaction** +# +# Depending on their angles, two non-zero time-lag interactions may have the +# same magnitude in the complex plane, but different magnitudes when projected +# to the imaginary axis. +# +# This is demonstrated below, where we simulate 2 interactions with non-zero +# time-lags at 10-12 Hz and 23-25 Hz. Computing the connectivity, we see how +# both interactions have a similar magnitude (~0.9), but different angles +# (~-45° for 10-12 Hz; ~-90° for 23-25 Hz). + +# %% + +# Generate simulated data +data_10_12 = make_signals_in_freq_bands( + n_seeds=3, + n_targets=3, + freq_band=(10, 12), # 10-12 Hz interaction + connection_delay=1, # samples + rng_seed=40, +) + +data_23_25 = make_signals_in_freq_bands( + n_seeds=3, + n_targets=3, + freq_band=(23, 25), # 23-25 Hz interaction + connection_delay=1, # samples + rng_seed=42, +) + +# Combine data into a single array +data = data_10_12.add_channels([data_23_25]) + +# Compute CaCoh & MIC +(cacoh, mic) = spectral_connectivity_epochs( + data, method=["cacoh", "mic"], indices=multivar_indices, sfreq=100, fmin=3, fmax=35 +) + +# Get complex connectivity values at frequency bands +freqs = cacoh.freqs +cacoh_10_12 = np.mean(cacoh.get_data()[0, freqs.index(10) : freqs.index(12) + 1]) +cacoh_23_25 = np.mean(cacoh.get_data()[0, freqs.index(23) : freqs.index(25) + 1]) + +# Plot complex connectivity values +fig, axis = plot_connectivity_circle() +axis.quiver( + 0, + 0, + np.real(cacoh_10_12), + np.imag(cacoh_10_12), + units="xy", + scale=1, + linewidth=2, + color="C2", + label="10-12 Hz", +) +axis.quiver( + 0, + 0, + np.real(cacoh_23_25), + np.imag(cacoh_23_25), + units="xy", + scale=1, + linewidth=2, + color="C3", + label="23-25 Hz", + zorder=99, +) +axis.legend(loc="upper right", bbox_to_anchor=[1.1, 1.1]) + +############################################################################### +# Plotting the connectivity values for CaCoh and MIC, we see how the 10-12 Hz +# and 23-25 Hz interactions have a similar magnitude for CaCoh, whereas the MIC +# scores for the 10-12 Hz interaction are lower than for the 23-25 Hz +# interaction. +# +# This difference reflects the fact that as the angle of interaction deviates +# from :math:`\pm` 90°, less information will be represented in the imaginary +# part of coherency. Accordingly, considering only the imaginary part of +# coherency can bias connectivity estimates based on the angle of interaction. + +# %% + +# Plot CaCoh & MIC +fig, axis = plt.subplots(1, 1) +axis.plot(cacoh.freqs, np.abs(cacoh.get_data()[0]), linewidth=2, label="CaCoh") +axis.plot( + mic.freqs, np.abs(mic.get_data()[0]), linewidth=2, label="MIC", linestyle="--" +) +axis.set_xlabel("Frequency (Hz)") +axis.set_ylabel("Connectivity (A.U.)") +axis.annotate("$\pm$45°\ninteraction", xy=(12.5, 0.9)) +axis.annotate("$\pm$90°\ninteraction", xy=(26.5, 0.9)) +axis.legend(loc="upper left") +fig.suptitle("CaCoh vs. MIC\n$\pm$45° & $\pm$90° interactions") + +############################################################################### +# In situations where non-physiological zero time-lag interactions are present, +# this phase angle-dependent bias is generally acceptable to avoid spurious +# connectivity estimates. However in situations where non-physiological zero +# time-lag interactions are not present, such a bias is clearly problematic. +# +# | +# +# Again, these considerations are not specific to multivariate methods, as +# shown below with Coh and ImCoh. + +# %% + +# Compute Coh & ImCoh +(coh, imcoh) = spectral_connectivity_epochs( + data, method=["coh", "imcoh"], indices=bivar_indices, sfreq=100, fmin=3, fmax=35 +) + +coh_mean = np.mean(coh.get_data(), axis=0) +imcoh_mean = np.mean(np.abs(imcoh.get_data()), axis=0) + +coh_mean_subbed = coh_mean - np.mean(coh_mean) +imcoh_mean_subbed = imcoh_mean - np.mean(imcoh_mean) + +# Plot Coh & ImCoh +fig, axis = plt.subplots(1, 1) +axis.plot(coh.freqs, coh_mean_subbed, linewidth=2, label="Coh") +axis.plot(imcoh.freqs, imcoh_mean_subbed, linewidth=2, label="ImCoh", linestyle="--") +axis.set_xlabel("Frequency (Hz)") +axis.set_ylabel("Mean-corrected connectivity (A.U.)") +axis.annotate("$\pm$45°\ninteraction", xy=(13, 0.25)) +axis.annotate("$\pm$90°\ninteraction", xy=(26.5, 0.25)) +axis.legend(loc="upper left") +fig.suptitle("Coh vs. ImCoh\n$\pm$45° & $\pm$90° interactions") + +############################################################################### +# Bivariate vs. multivariate coherency methods +# -------------------------------------------- +# +# As we have seen, coherency-based methods can be bivariate (Cohy, Coh, ImCoh) +# and multivariate (CaCoh, MIC, MIM). Whilst both forms capture the same +# information, there are several benefits to using multivariate methods when +# investigating connectivity between many signals. +# +# The multivariate methods can be used to capture the most relevant +# interactions between two groups of signals, representing this information in +# the component, rather than signal space. +# +# The dimensionality reduction associated with these methods offers: a much +# easier interpretation of the results; a higher signal-to-noise ratio compared +# to e.g. averaging bivariate connectivity estimates across multiple pairs of +# signals; and even reduced bias in what information is captured +# :footcite:`EwaldEtAl2012`. +# +# Furthermore, despite the dimensionality reduction of multivariate methods it +# is still possible to investigate the topographies of connectivity, with +# spatial patterns of connectivity being returned alongside the connectivity +# values themselves :footcite:`HaufeEtAl2014`. +# +# More information about the multivariate coherency-based methods can be found +# in the following examples: +# +# - CaCoh - :doc:`cacoh` +# - MIC & MIM - :doc:`mic_mim` + +############################################################################### +# Alternative approaches to computing connectivity +# ------------------------------------------------ +# +# Coherency-based methods are only some of the many approaches available in +# MNE-Connectivity for studying interactions between signals. Other +# non-directed measures include those based on the phase-lag index +# :footcite:`StamEtAl2007,VinckEtAl2011` (see also :doc:`dpli_wpli_pli`) and +# phase locking value :footcite:`LachauxEtAl1999,BrunaEtAl2018`. +# +# Furthermore, directed measures of connectivity which determine the direction +# of information flow are also available, including a variant of the phase-lag +# index :footcite:`StamEtAl2012` (see also :doc:`dpli_wpli_pli`), the phase +# slope index :footcite:`NolteEtAl2008` (see also +# :func:`mne_connectivity.phase_slope_index`), and Granger causality +# :footcite:`BarnettSeth2015,WinklerEtAl2016` (see also +# :doc:`granger_causality`). + +############################################################################### +# Conclusion +# ---------- +# +# Altogether, there are clear scenarious in which different coherency-based +# methods are appropriate. +# +# Methods based on the imaginary part of coherency alone (ImCoh, MIC, MIM) +# should be used when non-physiological zero time-lag interactions are present. +# +# In contrast, methods based on the real and imaginary parts of coherency +# (Cohy, Coh, CaCoh) should be used when non-physiological zero time-lag +# interactions are absent. + +############################################################################### +# References +# ---------- +# .. footbibliography:: + +# %% diff --git a/examples/mic_mim.py b/examples/mic_mim.py index 244e045b..67c6e6d3 100644 --- a/examples/mic_mim.py +++ b/examples/mic_mim.py @@ -42,17 +42,32 @@ # A popular bivariate measure of connectivity is the imaginary part of # coherency, which looks at the correlation between two signals in the # frequency domain and is immune to spurious connectivity arising from volume -# conduction artefacts :footcite:`NolteEtAl2004`. However, depending on the -# degree of source mixing, this measure is susceptible to biased estimates of +# conduction artefacts :footcite:`NolteEtAl2004`. However, in cases where +# interactions between multiple signals are of interest, computing connectivity +# between all possible combinations of signals leads to a very large number of +# results which is difficult to interpret. A common approach is to average +# results across these connections, however this risks reducing the +# signal-to-noise ratio of results and burying interactions that are present +# between only a small number of channels. +# +# Additionally, this bivariate measure is susceptible to biased estimates of # connectivity based on the spatial proximity of sensors -# :footcite:`EwaldEtAl2012`. +# :footcite:`EwaldEtAl2012` depending on the degree of source mixing in the +# signals. +# +# To overcome this limitation, spatial filters derived from eigendecompositions +# allows connectivity to be analysed in a multivariate manner, removing the +# source mixing-dependent bias and increase the signal-to-noise ratio of +# connectivity estimates :footcite:`EwaldEtAl2012`. This approach goes beyond +# simply aggregating information across all possible combinations of signals, +# extracting the underlying components of connectivity in a frequency-resolved +# manner. # -# To overcome this limitation, spatial filters can be used to estimate -# connectivity free from this source mixing-dependent bias, which additionally -# increases the signal-to-noise ratio and allows signals to be analysed in a -# multivariate manner :footcite:`EwaldEtAl2012`. This approach leads to the -# following methods: the maximised imaginary part of coherency (MIC); and the -# multivariate interaction measure (MIM). +# This leads to the following methods: the maximised imaginary part of +# coherency (MIC); and the multivariate interaction measure (MIM). These +# methods are similar to the multivariate method based on coherency (CaCoh +# :footcite:`VidaurreEtAl2019`; see :doc:`cacoh` and +# :doc:`compare_coherency_methods`). # # We start by loading some example MEG data and dividing it into # two-second-long epochs. @@ -121,8 +136,9 @@ # eigendecomposition of information from the cross-spectral density (Eq. 7 of # :footcite:`EwaldEtAl2012`): # -# :math:`MIC=\frac{\boldsymbol{\alpha}^T \boldsymbol{E \beta}}{\parallel -# \boldsymbol{\alpha}\parallel \parallel\boldsymbol{\beta}\parallel}`, +# :math:`\textrm{MIC}=\frac{\boldsymbol{\alpha}^T \boldsymbol{E \beta}} +# {\parallel\boldsymbol{\alpha}\parallel \parallel\boldsymbol{\beta} +# \parallel}`, # # where :math:`\boldsymbol{\alpha}` and :math:`\boldsymbol{\beta}` are the # spatial filters for the seeds and targets, respectively, and @@ -155,13 +171,13 @@ ############################################################################### # Furthermore, spatial patterns of connectivity can be constructed from the -# spatial filters to give a picture of the location of the sources involved in -# the connectivity. This information is stored under ``attrs['patterns']`` of -# the connectivity class, with one value per frequency for each channel in the -# seeds and targets. As with MIC, the absolute value of the patterns reflect -# the strength, however the sign differences can be used to visualise the -# orientation of the underlying dipole sources. The spatial patterns are -# **not** bound between :math:`[-1, 1]`. +# spatial filters to give a picture of the location of the channels involved in +# the connectivity :footcite:`HaufeEtAl2014`. This information is stored under +# ``attrs['patterns']`` of the connectivity class, with one value per frequency +# for each channel in the seeds and targets. As with MIC, the absolute value of +# the patterns reflect the strength, however the sign differences can be used +# to visualise the orientation of the underlying dipole sources. The spatial +# patterns are **not** bound between :math:`[-1, 1]`. # # Here, we average across the patterns in the 13-18 Hz range. Plotting the # patterns shows that the greatest connectivity between the left and right @@ -170,6 +186,11 @@ # values, we can infer the existence of a dipole source between the central and # posterior regions of the left hemisphere accounting for the connectivity # contributions (represented on the plot as a green line). +# +# **Note:** The spatial patterns are not a substitute for source +# reconstruction. If you need the spatial patterns in source space, you should +# perform source reconstruction before computing connectivity (see e.g. +# :doc:`mne_inverse_coherence_epochs`). # %% @@ -242,7 +263,7 @@ # component explicitly, and instead the desired result can be achieved from # :math:`E` alone (Eq. 14 of :footcite:`EwaldEtAl2012`): # -# :math:`MIM=tr(\boldsymbol{EE}^T)`, +# :math:`\textrm{MIM}=tr(\boldsymbol{EE}^T)`, # # where again the frequency dependence is omitted. Unlike MIC, MIM is # positive-valued and can be > 1. Without normalisation, MIM can be @@ -371,8 +392,8 @@ # compare standard and rank subspace-projected MIM fig, axis = plt.subplots(1, 1) -axis.plot(mim_red.freqs, mim_red_meansub, linewidth=2, label="rank subspace (25) MIM") axis.plot(mim.freqs, mim_meansub, linewidth=2, label="standard MIM") +axis.plot(mim_red.freqs, mim_red_meansub, linewidth=2, label="rank subspace (25) MIM") axis.set_xlabel("Frequency (Hz)") axis.set_ylabel("Mean-corrected connectivity (A.U.)") axis.legend() @@ -402,7 +423,8 @@ # gets the singular values of the data s = np.linalg.svd(raw.get_data(), compute_uv=False) # finds how many singular values are 'close' to the largest singular value -rank = np.count_nonzero(s >= s[0] * 1e-4) # 1e-4 is the 'closeness' criteria +rank = np.count_nonzero(s >= s[0] * 1e-4) # 1e-4 is the 'closeness' criteria, which is +# a hyper-parameter ############################################################################### @@ -411,27 +433,24 @@ # # These multivariate methods offer many benefits in the form of dimensionality # reduction, signal-to-noise ratio improvements, and invariance to -# estimate-biasing source mixing; however, no method is perfect. The immunity -# of the imaginary part of coherency to volume conduction comes from the fact -# that these artefacts have zero phase lag, and hence a zero-valued imaginary -# component. By projecting the complex-valued coherency to the imaginary axis, -# signals of a given magnitude with phase lag differences close to 90° or 270° -# see their contributions to the connectivity estimate increased relative to -# comparable signals with phase lag differences close to 0° or 180°. Therefore, -# the imaginary part of coherency is biased towards connectivity involving 90° -# and 270° phase lag difference components. +# estimate-biasing source mixing; however, no method is perfect. Important +# considerations must be taken into account when choosing methods based on the +# imaginary part of coherency such as MIC or MIM versus those based on +# coherency/coherence, such as CaCoh. +# +# In short, if you want to examine connectivity between signals from the same +# modality, you should consider using MIC and MIM to avoid spurious +# connectivity estimates stemming from e.g. volume conduction artefacts. # -# Whilst this is not a limitation specific to the multivariate extension of -# this measure, these multivariate methods can introduce further bias: when -# maximising the imaginary part of coherency, components with phase lag -# differences close to 90° and 270° will likely give higher connectivity -# estimates, and so may be prioritised by the spatial filters. +# On the other hand, if you want to examine connectivity between signals from +# different modalities, CaCoh is a more appropriate method than MIC/MIM. This +# is because voilume conduction artefacts are of less concern, and CaCoh does +# not risk biasing connectivity estimates towards interactions with particular +# phase lags like MIC/MIM. # -# Such a limitation should be kept in mind when estimating connectivity using -# these methods. Possible sanity checks can involve comparing the spectral -# profiles of MIC/MIM to coherence and the imaginary part of coherency -# computed on the same data, as well as comparing to other multivariate -# measures, such as canonical coherence :footcite:`VidaurreEtAl2019`. +# These scenarios are described in more detail in the +# :doc:`compare_coherency_methods` example. + ############################################################################### # References diff --git a/mne_connectivity/base.py b/mne_connectivity/base.py index fdb03884..e38fd2c2 100644 --- a/mne_connectivity/base.py +++ b/mne_connectivity/base.py @@ -725,7 +725,7 @@ def get_data(self, output="compact"): if output == "raveled": data = self._data else: - if self.method in ["mic", "mim", "gc", "gc_tr"]: + if self.method in ["cacoh", "mic", "mim", "gc", "gc_tr"]: # multivariate results cannot be returned in a dense form as a # single set of results would correspond to multiple entries in # the matrix, and there could also be cases where multiple diff --git a/mne_connectivity/spectral/epochs.py b/mne_connectivity/spectral/epochs.py index aa718455..1cfab73e 100644 --- a/mne_connectivity/spectral/epochs.py +++ b/mne_connectivity/spectral/epochs.py @@ -644,12 +644,13 @@ def spectral_connectivity_epochs( %(names)s method : str | list of str Connectivity measure(s) to compute. These can be ``['coh', 'cohy', - 'imcoh', 'mic', 'mim', 'plv', 'ciplv', 'ppc', 'pli', 'dpli', 'wpli', - 'wpli2_debiased', 'gc', 'gc_tr']``. These are: + 'imcoh', 'cacoh', 'mic', 'mim', 'plv', 'ciplv', 'ppc', 'pli', 'dpli', + 'wpli', 'wpli2_debiased', 'gc', 'gc_tr']``. These are: * %(coh)s * %(cohy)s * %(imcoh)s + * %(cacoh)s * %(mic)s * %(mim)s * %(plv)s @@ -663,8 +664,8 @@ def spectral_connectivity_epochs( * %(gc)s * %(gc_tr)s - Multivariate methods (``['mic', 'mim', 'gc', 'gc_tr]``) cannot be - called with the other methods. + Multivariate methods (``['cacoh', 'mic', 'mim', 'gc', 'gc_tr']``) + cannot be called with the other methods. indices : tuple of array | None Two arrays with indices of connections for which to compute connectivity. If a bivariate method is called, each array for the seeds @@ -726,7 +727,7 @@ def spectral_connectivity_epochs( Two arrays with the rank to project the seed and target data to, respectively, using singular value decomposition. If None, the rank of the data is computed and projected to. Only used if ``method`` contains - any of ``['mic', 'mim', 'gc', 'gc_tr']``. + any of ``['cacoh', 'mic', 'mim', 'gc', 'gc_tr']``. block_size : int How many connections to compute at once (higher numbers are faster but require more memory). @@ -832,24 +833,42 @@ def spectral_connectivity_epochs( C = ---------------------- sqrt(E[Sxx] * E[Syy]) + 'cacoh' : Canonical Coherency (CaCoh) :footcite:`VidaurreEtAl2019` + given by: + + :math:`\textrm{CaCoh}=\Large{\frac{\boldsymbol{a}^T\boldsymbol{D} + (\Phi)\boldsymbol{b}}{\sqrt{\boldsymbol{a}^T\boldsymbol{a} + \boldsymbol{b}^T\boldsymbol{b}}}}` + + where: :math:`\boldsymbol{D}(\Phi)` is the cross-spectral density + between seeds and targets transformed for a given phase angle + :math:`\Phi`; and :math:`\boldsymbol{a}` and :math:`\boldsymbol{b}` + are eigenvectors for the seeds and targets, such that + :math:`\boldsymbol{a}^T\boldsymbol{D}(\Phi)\boldsymbol{b}` + maximises coherency between the seeds and targets. Taking the + absolute value of the results gives maximised coherence. + 'mic' : Maximised Imaginary part of Coherency (MIC) :footcite:`EwaldEtAl2012` given by: - :math:`MIC=\Large{\frac{\boldsymbol{\alpha}^T \boldsymbol{E \beta}} - {\parallel\boldsymbol{\alpha}\parallel \parallel\boldsymbol{\beta} - \parallel}}` + :math:`\textrm{MIC}=\Large{\frac{\boldsymbol{\alpha}^T + \boldsymbol{E \beta}}{\parallel\boldsymbol{\alpha}\parallel + \parallel\boldsymbol{\beta}\parallel}}` where: :math:`\boldsymbol{E}` is the imaginary part of the transformed cross-spectral density between seeds and targets; and :math:`\boldsymbol{\alpha}` and :math:`\boldsymbol{\beta}` are eigenvectors for the seeds and targets, such that - :math:`\boldsymbol{\alpha}^T \boldsymbol{E \beta}` maximises - connectivity between the seeds and targets. + :math:`\boldsymbol{\alpha}^T \boldsymbol{E \beta}` maximises the + imaginary part of coherency between the seeds and targets. 'mim' : Multivariate Interaction Measure (MIM) :footcite:`EwaldEtAl2012` given by: - :math:`MIM=tr(\boldsymbol{EE}^T)` + :math:`\textrm{MIM}=tr(\boldsymbol{EE}^T)` + + where :math:`\boldsymbol{E}` is the imaginary part of the + transformed cross-spectral density between seeds and targets. 'plv' : Phase-Locking Value (PLV) :footcite:`LachauxEtAl1999` given by:: @@ -893,7 +912,7 @@ def spectral_connectivity_epochs( :math:`GC = ln\Large{(\frac{\lvert\boldsymbol{S}_{tt}\rvert}{\lvert \boldsymbol{S}_{tt}-\boldsymbol{H}_{ts}\boldsymbol{\Sigma}_{ss - \lvert t}\boldsymbol{H}_{ts}^*\rvert}})`, + \lvert t}\boldsymbol{H}_{ts}^*\rvert}})` where: :math:`s` and :math:`t` represent the seeds and targets, respectively; :math:`\boldsymbol{H}` is the spectral transfer diff --git a/mne_connectivity/spectral/epochs_multivariate.py b/mne_connectivity/spectral/epochs_multivariate.py index 2881cc00..1e4e8e3b 100644 --- a/mne_connectivity/spectral/epochs_multivariate.py +++ b/mne_connectivity/spectral/epochs_multivariate.py @@ -4,6 +4,8 @@ # Thomas S. Binns # Tien D. Nguyen # Richard M. Köhler +# Mohammad Orabe +# Mina Jamshidi Idaji # # License: BSD (3-clause) @@ -103,6 +105,7 @@ class _EpochMeanMultivariateConEstBase(_AbstractConEstBase): n_steps = None patterns = None + con_scores_dtype = np.float64 def __init__(self, n_signals, n_cons, n_freqs, n_times, n_jobs=1): self.n_signals = n_signals @@ -114,10 +117,14 @@ def __init__(self, n_signals, n_cons, n_freqs, n_times, n_jobs=1): # include time dimension, even when unused for indexing flexibility if n_times == 0: self.csd_shape = (n_signals**2, n_freqs) - self.con_scores = np.zeros((n_cons, n_freqs, 1)) + self.con_scores = np.zeros( + (n_cons, n_freqs, 1), dtype=self.con_scores_dtype + ) else: self.csd_shape = (n_signals**2, n_freqs, n_times) - self.con_scores = np.zeros((n_cons, n_freqs, n_times)) + self.con_scores = np.zeros( + (n_cons, n_freqs, n_times), dtype=self.con_scores_dtype + ) # allocate space for accumulation of CSD self._acc = np.zeros(self.csd_shape, dtype=np.complex128) @@ -170,10 +177,14 @@ def reshape_csd(self): class _MultivariateCohEstBase(_EpochMeanMultivariateConEstBase): - """Base estimator for multivariate imag. part of coherency methods. - - See Ewald et al. (2012). NeuroImage. DOI: 10.1016/j.neuroimage.2011.11.084 - for equation references. + """Base estimator for multivariate coherency methods. + + See: + - Imaginary part of coherency, i.e. multivariate imaginary part of + coherency (MIC) and multivariate interaction measure (MIM): Ewald et al. + (2012). NeuroImage. DOI: 10.1016/j.neuroimage.2011.11.084 + - Coherency/coherence, i.e. canonical coherency (CaCoh): Vidaurre et al. + (2019). NeuroImage. DOI: 10.1016/j.neuroimage.2019.116009 """ name: Optional[str] = None @@ -185,8 +196,8 @@ def __init__(self, n_signals, n_cons, n_freqs, n_times, n_jobs=1): ) def compute_con(self, indices, ranks, n_epochs=1): - """Compute multivariate imag. part of coherency between signals.""" - assert self.name in ["MIC", "MIM"], ( + """Compute multivariate coherency methods.""" + assert self.name in ["CaCoh", "MIC", "MIM"], ( "the class name is not recognised, please contact the " "mne-connectivity developers" ) @@ -196,7 +207,7 @@ def compute_con(self, indices, ranks, n_epochs=1): times = np.arange(n_times) freqs = np.arange(self.n_freqs) - if self.name == "MIC": + if self.name in ["CaCoh", "MIC"]: self.patterns = np.full( (2, self.n_cons, indices[0].shape[1], self.n_freqs, n_times), np.nan ) @@ -213,20 +224,14 @@ def compute_con(self, indices, ranks, n_epochs=1): C = csd[np.ix_(times, freqs, con_idcs, con_idcs)] - # Eqs. 32 & 33 + # Eqs. 32 & 33 of Ewald et al.; Eq. 15 of Vidaurre et al. C_bar, U_bar_aa, U_bar_bb = self._csd_svd( C, seed_idcs, seed_rank, target_rank ) - # Eqs. 3 & 4 - E = self._compute_e(C_bar, n_seeds=U_bar_aa.shape[3]) - - if self.name == "MIC": - self._compute_mic( - E, C, seed_idcs, target_idcs, n_times, U_bar_aa, U_bar_bb, con_i - ) - else: - self._compute_mim(E, seed_idcs, target_idcs, con_i) + self._compute_con_daughter( + seed_idcs, target_idcs, C, C_bar, U_bar_aa, U_bar_bb, con_i + ) con_i += 1 @@ -243,7 +248,7 @@ def _csd_svd(self, csd, seed_idcs, seed_rank, target_rank): C_bb = csd[..., n_seeds:, n_seeds:] C_ba = csd[..., n_seeds:, :n_seeds] - # Eq. 32 + # Eqs. 32 (Ewald et al.) & 15 (Vidaurre et al.) if seed_rank != n_seeds: U_aa = np.linalg.svd(np.real(C_aa), full_matrices=False)[0] U_bar_aa = U_aa[..., :seed_rank] @@ -260,7 +265,7 @@ def _csd_svd(self, csd, seed_idcs, seed_rank, target_rank): np.identity(n_targets), (n_times, self.n_freqs) + (n_targets, n_targets) ) - # Eq. 33 + # Eq. 33 (Ewald et al.) C_bar_aa = np.matmul(U_bar_aa.transpose(0, 1, 3, 2), np.matmul(C_aa, U_bar_aa)) C_bar_ab = np.matmul(U_bar_aa.transpose(0, 1, 3, 2), np.matmul(C_ab, U_bar_bb)) C_bar_bb = np.matmul(U_bar_bb.transpose(0, 1, 3, 2), np.matmul(C_bb, U_bar_bb)) @@ -273,20 +278,29 @@ def _csd_svd(self, csd, seed_idcs, seed_rank, target_rank): return C_bar, U_bar_aa, U_bar_bb - def _compute_e(self, csd, n_seeds): - """Compute E from the CSD.""" - C_r = np.real(csd) + def _compute_con_daughter( + self, seed_idcs, target_idcs, C, C_bar, U_bar_aa, U_bar_bb, con_i + ): + """Compute multivariate coherency for one connection. + + An empty method to be implemented by subclasses. + """ + + def _compute_t(self, C_r, n_seeds): + """Compute transformation matrix, T, for frequencies (in parallel). - parallel, parallel_compute_t, _ = parallel_func( - _mic_mim_compute_t, self.n_jobs, verbose=False + Eq. 3 of Ewald et al.; part of Eq. 9 of Vidaurre et al. + """ + parallel, parallel_invsqrtm, _ = parallel_func( + _invsqrtm, self.n_jobs, verbose=False ) # imag. part of T filled when data is rank-deficient - T = np.zeros(csd.shape, dtype=np.complex128) + T = np.zeros(C_r.shape, dtype=np.complex128) for block_i in ProgressBar(range(self.n_steps), mesg="frequency blocks"): freqs = self._get_block_indices(block_i, self.n_freqs) T[:, freqs] = np.array( - parallel(parallel_compute_t(C_r[:, f], T[:, f], n_seeds) for f in freqs) + parallel(parallel_invsqrtm(C_r[:, f], T[:, f], n_seeds) for f in freqs) ).transpose(1, 0, 2, 3) if not np.isreal(T).all() or not np.isfinite(T).all(): @@ -296,20 +310,96 @@ def _compute_e(self, csd, n_seeds): "using full rank data or specify an appropriate rank for the " "seeds and targets that is less than or equal to their ranks" ) - T = np.real(T) # make T real if check passes + + return np.real(T) # make T real if check passes + + def reshape_results(self): + """Remove time dimension from results, if necessary.""" + if self.n_times == 0: + self.con_scores = self.con_scores[..., 0] + if self.patterns is not None: + self.patterns = self.patterns[..., 0] + + +def _invsqrtm(C, T, n_seeds): + """Compute inverse sqrt of CSD over times (used for CaCoh, MIC, & MIM). + + Parameters + ---------- + C : np.ndarray, shape=(n_times, n_channels, n_channels) + CSD for a single frequency and all times (n_times=1 if the mode is not + time-frequency resolved, e.g. multitaper). + T : np.ndarray, shape=(n_times, n_channels, n_channels) + Empty array to store the inverse square root of the CSD in. + n_seeds : int + Number of seed channels for the connection. + + Returns + ------- + T : np.ndarray, shape=(n_times, n_channels, n_channels) + Inverse square root of the CSD. Name comes from Ewald et al. (2012). + + Notes + ----- + Kept as a standalone function to allow for parallelisation over CSD + frequencies. + + See Eq. 3 of Ewald et al. (2012). NeuroImage. DOI: + 10.1016/j.neuroimage.2011.11.084. + """ + for time_i in range(C.shape[0]): + T[time_i, :n_seeds, :n_seeds] = sp.linalg.fractional_matrix_power( + C[time_i, :n_seeds, :n_seeds], -0.5 + ) + T[time_i, n_seeds:, n_seeds:] = sp.linalg.fractional_matrix_power( + C[time_i, n_seeds:, n_seeds:], -0.5 + ) + + return T + + +class _MultivariateImCohEstBase(_MultivariateCohEstBase): + """Base estimator for multivariate imag. part of coherency methods. + + See Ewald et al. (2012). NeuroImage. DOI: 10.1016/j.neuroimage.2011.11.084 + for equation references. + """ + + def _compute_con_daughter( + self, seed_idcs, target_idcs, C, C_bar, U_bar_aa, U_bar_bb, con_i + ): + """Compute multivariate imag. part of coherency for one connection.""" + assert self.name in ["MIC", "MIM"], ( + "the class name is not recognised, please contact the " + "mne-connectivity developers" + ) + + # Eqs. 3 & 4 + E = self._compute_e(C_bar, n_seeds=U_bar_aa.shape[3]) + + if self.name == "MIC": + self._compute_mic(E, C, seed_idcs, target_idcs, U_bar_aa, U_bar_bb, con_i) + else: + self._compute_mim(E, seed_idcs, target_idcs, con_i) + + def _compute_e(self, C, n_seeds): + """Compute E from the CSD.""" + C_r = np.real(C) + + # Eq. 3 + T = self._compute_t(C_r, n_seeds) # Eq. 4 - D = np.matmul(T, np.matmul(csd, T)) + D = np.matmul(T, np.matmul(C, T)) # E as imag. part of D between seeds and targets return np.imag(D[..., :n_seeds, n_seeds:]) - def _compute_mic( - self, E, C, seed_idcs, target_idcs, n_times, U_bar_aa, U_bar_bb, con_i - ): - """Compute MIC and the associated spatial patterns.""" + def _compute_mic(self, E, C, seed_idcs, target_idcs, U_bar_aa, U_bar_bb, con_i): + """Compute MIC & spatial patterns for one connection.""" n_seeds = len(seed_idcs) n_targets = len(target_idcs) + n_times = C.shape[0] times = np.arange(n_times) freqs = np.arange(self.n_freqs) @@ -372,7 +462,7 @@ def _compute_mic( ).T def _compute_mim(self, E, seed_idcs, target_idcs, con_i): - """Compute MIM (a.k.a. GIM if seeds == targets).""" + """Compute MIM (a.k.a. GIM if seeds == targets) for one connection.""" # Eq. 14 self.con_scores[con_i] = ( np.matmul(E, E.transpose(0, 1, 3, 2)).trace(axis1=2, axis2=3).T @@ -384,37 +474,182 @@ def _compute_mim(self, E, seed_idcs, target_idcs, con_i): ): self.con_scores[con_i] *= 0.5 - def reshape_results(self): - """Remove time dimension from results, if necessary.""" - if self.n_times == 0: - self.con_scores = self.con_scores[..., 0] - if self.patterns is not None: - self.patterns = self.patterns[..., 0] +class _MICEst(_MultivariateImCohEstBase): + """Multivariate imaginary part of coherency (MIC) estimator.""" -def _mic_mim_compute_t(C, T, n_seeds): - """Compute T for a single frequency (used for MIC and MIM).""" - for time_i in range(C.shape[0]): - T[time_i, :n_seeds, :n_seeds] = sp.linalg.fractional_matrix_power( - C[time_i, :n_seeds, :n_seeds], -0.5 + name = "MIC" + + +class _MIMEst(_MultivariateImCohEstBase): + """Multivariate interaction measure (MIM) estimator.""" + + name = "MIM" + + +class _CaCohEst(_MultivariateCohEstBase): + """Canonical coherence (CaCoh) estimator. + + See Vidaurre et al. (2019). NeuroImage. DOI: + 10.1016/j.neuroimage.2019.116009 for equation references. + """ + + name = "CaCoh" + con_scores_dtype = np.complex128 # CaCoh is complex-valued + + def _compute_con_daughter( + self, seed_idcs, target_idcs, C, C_bar, U_bar_aa, U_bar_bb, con_i + ): + """Compute CaCoh & spatial patterns for one connection.""" + assert self.name == "CaCoh", ( + "the class name is not recognised, please contact the " + "mne-connectivity developers" ) - T[time_i, n_seeds:, n_seeds:] = sp.linalg.fractional_matrix_power( - C[time_i, n_seeds:, n_seeds:], -0.5 + n_seeds = len(seed_idcs) + n_targets = len(target_idcs) + + rank_seeds = U_bar_aa.shape[3] # n_seeds after SVD + + C_bar_ab = C_bar[..., :rank_seeds, rank_seeds:] + + # Same as Eq. 3 of Ewald et al. (2012) + T = self._compute_t(np.real(C_bar), n_seeds=rank_seeds) + T_aa = T[..., :rank_seeds, :rank_seeds] # left term in Eq. 9 + T_bb = T[..., rank_seeds:, rank_seeds:] # right term in Eq. 9 + + max_coh, max_phis = self._first_optimise_phi(C_bar_ab, T_aa, T_bb) + + max_coh, max_phis = self._final_optimise_phi( + C_bar_ab, T_aa, T_bb, max_coh, max_phis ) - return T + # Store connectivity score as complex value + self.con_scores[con_i] = (max_coh * np.exp(-1j * max_phis)).T + + self._compute_patterns( + max_phis, + C, + C_bar_ab, + T_aa, + T_bb, + U_bar_aa, + U_bar_bb, + n_seeds, + n_targets, + con_i, + ) + def _first_optimise_phi(self, C_ab, T_aa, T_bb): + """Find the rough angle, phi, at which coherence is maximised.""" + n_iters = 5 -class _MICEst(_MultivariateCohEstBase): - """Multivariate imaginary part of coherency (MIC) estimator.""" + # starting phi values to optimise over (in radians) + phis = np.linspace(np.pi / n_iters, np.pi, n_iters) + phis_coh = np.zeros((n_iters, *C_ab.shape[:2])) + for iter_i, iter_phi in enumerate(phis): + phi = np.full(C_ab.shape[:2], fill_value=iter_phi) + phis_coh[iter_i] = self._compute_cacoh(phi, C_ab, T_aa, T_bb) - name = "MIC" + return np.max(phis_coh, axis=0), phis[np.argmax(phis_coh, axis=0)] + def _final_optimise_phi(self, C_ab, T_aa, T_bb, max_coh, max_phis): + """Fine-tune the angle at which coherence is maximised. -class _MIMEst(_MultivariateCohEstBase): - """Multivariate interaction measure (MIM) estimator.""" + Uses a 2nd order Taylor expansion to approximate change in coherence + w.r.t. phi, and determining the next phi to evaluate coherence on (over + a total of 10 iterations). - name = "MIM" + Depending on how the new phi affects coherence, the step size for the + subsequent iteration is adjusted, like that in the Levenberg-Marquardt + algorithm. + + Each time-freq. entry of coherence has its own corresponding phi. + """ + n_iters = 10 # sufficient for (close to) exact solution + delta_phi = 1e-6 + mus = np.ones_like(max_phis) # optimisation step size + + for _ in range(n_iters): + # 2nd order Taylor expansion around phi + coh_plus = self._compute_cacoh(max_phis + delta_phi, C_ab, T_aa, T_bb) + coh_minus = self._compute_cacoh(max_phis - delta_phi, C_ab, T_aa, T_bb) + f_prime = (coh_plus - coh_minus) / (2 * delta_phi) + f_prime_prime = (coh_plus + coh_minus - 2 * max_coh) / (delta_phi**2) + + # determine new phi to test + phis = max_phis + (-f_prime / (f_prime_prime - mus)) + # bound phi in range [-pi, pi] + phis = np.mod(phis + np.pi / 2, np.pi) - np.pi / 2 + + coh = self._compute_cacoh(phis, C_ab, T_aa, T_bb) + + # find where new phi increases coh & update these values + greater_coh = coh > max_coh + max_coh[greater_coh] = coh[greater_coh] + max_phis[greater_coh] = phis[greater_coh] + + # update step size + mus[greater_coh] /= 2 + mus[~greater_coh] *= 2 + + return max_coh, phis + + def _compute_cacoh(self, phis, C_ab, T_aa, T_bb): + """Compute the maximum coherence for a given set of phis.""" + # from numerator of Eq. 5 + # for a given CSD entry, projects it onto a span with angle phi, such + # that the magnitude of the projected line is captured in the real part + C_ab = np.real(np.exp(-1j * np.expand_dims(phis, axis=(2, 3))) * C_ab) + + # Eq. 9; T_aa/bb is sqrt(inv(real(C_aa/bb))) + D = np.matmul(T_aa, np.matmul(C_ab, T_bb)) + + # Eq. 12 + a = np.linalg.eigh(np.matmul(D, D.transpose(0, 1, 3, 2)))[1][..., -1] + b = np.linalg.eigh(np.matmul(D.transpose(0, 1, 3, 2), D))[1][..., -1] + + # Eq. 8 + numerator = np.einsum( + "ijk,ijk->ij", a, np.matmul(D, np.expand_dims(b, axis=3))[..., 0] + ) + denominator = np.sqrt( + np.einsum("ijk,ijk->ij", a, a) * np.einsum("ijk,ijk->ij", b, b) + ) + + return np.abs(numerator / denominator) + + def _compute_patterns( + self, + phis, + C, + C_bar_ab, + T_aa, + T_bb, + U_bar_aa, + U_bar_bb, + n_seeds, + n_targets, + con_i, + ): + """Compute CaCoh spatial patterns for the optimised phi.""" + C_bar_ab = np.real(np.exp(-1j * np.expand_dims(phis, axis=(2, 3))) * C_bar_ab) + D = np.matmul(T_aa, np.matmul(C_bar_ab, T_bb)) + a = np.linalg.eigh(np.matmul(D, D.transpose(0, 1, 3, 2)))[1][..., -1] + b = np.linalg.eigh(np.matmul(D.transpose(0, 1, 3, 2), D))[1][..., -1] + + # Eq. 7 rearranged - multiply both sides by sqrt(inv(real(C_aa/bb))) + alpha = np.matmul(T_aa, np.expand_dims(a, axis=3)) # filter for seeds + beta = np.matmul(T_bb, np.expand_dims(b, axis=3)) # filter for targets + + # Eq. 14; U_bar inclusion follows Eqs. 46 & 47 of Ewald et al. (2012) + # seed spatial patterns + self.patterns[0, con_i, :n_seeds] = ( + np.matmul(np.real(C[..., :n_seeds, :n_seeds]), np.matmul(U_bar_aa, alpha)) + )[..., 0].T + # target spatial patterns + self.patterns[1, con_i, :n_targets] = ( + np.matmul(np.real(C[..., n_seeds:, n_seeds:]), np.matmul(U_bar_bb, beta)) + )[..., 0].T class _GCEstBase(_EpochMeanMultivariateConEstBase): @@ -811,12 +1046,13 @@ class _GCTREst(_GCEstBase): # map names to estimator types _CON_METHOD_MAP_MULTIVARIATE = { + "cacoh": _CaCohEst, "mic": _MICEst, "mim": _MIMEst, "gc": _GCEst, "gc_tr": _GCTREst, } -_multivariate_methods = ["mic", "mim", "gc", "gc_tr"] +_multivariate_methods = ["cacoh", "mic", "mim", "gc", "gc_tr"] _gc_methods = ["gc", "gc_tr"] -_patterns_methods = ["mic"] # methods with spatial patterns +_patterns_methods = ["cacoh", "mic"] # methods with spatial patterns diff --git a/mne_connectivity/spectral/tests/data/README.md b/mne_connectivity/spectral/tests/data/README.md index ea9da2bd..4c62ebb6 100644 --- a/mne_connectivity/spectral/tests/data/README.md +++ b/mne_connectivity/spectral/tests/data/README.md @@ -1,7 +1,9 @@ -Author: Thomas S. Binns +Authors: +- Thomas S. Binns +- Mohammad Orabe The files found here are used for the regression test of the multivariate -connectivity methods for MIC, MIM, GC, and TRGC +connectivity methods for CaCoh, MIC, MIM, GC, and TRGC (`test_multivariate_spectral_connectivity_epochs_regression()` of `test_spectral.py`). @@ -9,8 +11,8 @@ connectivity methods for MIC, MIM, GC, and TRGC data with 15 epochs and 200 timepoints per epoch. Connectivity was computed in MATLAB using the original implementations of these methods and saved as a dictionary in `example_multivariate_matlab_results.pkl`. A publicly-available -implementation of the methods in MATLAB can be found here: -https://github.com/sccn/roiconnect. +implementation of the methods in MATLAB (except CaCoh) can be found here: +https://github.com/sccn/roiconnect. As the MNE code for computing the cross-spectral density matrix is not available in MATLAB, the CSD matrix was computed using MNE and then loaded into diff --git a/mne_connectivity/spectral/tests/data/example_multivariate_matlab_results.pkl b/mne_connectivity/spectral/tests/data/example_multivariate_matlab_results.pkl index 7b1de665..85bc2f84 100644 Binary files a/mne_connectivity/spectral/tests/data/example_multivariate_matlab_results.pkl and b/mne_connectivity/spectral/tests/data/example_multivariate_matlab_results.pkl differ diff --git a/mne_connectivity/spectral/tests/test_spectral.py b/mne_connectivity/spectral/tests/test_spectral.py index 8c493cb3..eba8120a 100644 --- a/mne_connectivity/spectral/tests/test_spectral.py +++ b/mne_connectivity/spectral/tests/test_spectral.py @@ -466,7 +466,7 @@ def test_spectral_connectivity(method, mode): assert out_lens[0] == 10 -@pytest.mark.parametrize("method", ["mic", "mim", "gc"]) +@pytest.mark.parametrize("method", ["cacoh", "mic", "mim", "gc"]) def test_spectral_connectivity_epochs_multivariate(method): """Test over-epoch multivariate connectivity methods.""" mode = "multitaper" # stick with single mode in interest of time @@ -519,7 +519,7 @@ def test_spectral_connectivity_epochs_multivariate(method): freqs.index(fend + trans_bandwidth * 2) + 1, ) - if method in ["mic", "mim"]: + if method in ["cacoh", "mic", "mim"]: lower_t = 0.2 upper_t = 0.5 @@ -570,27 +570,25 @@ def test_spectral_connectivity_epochs_multivariate(method): assert np.allclose(trgc[0, : bidx[0]].mean(), 0, atol=lower_t) assert np.allclose(trgc[0, bidx[1] :].mean(), 0, atol=lower_t) - # check all-to-all conn. computed for MIC/MIM when no indices given - if method in ["mic", "mim"]: + # check all-to-all conn. computed for CaCoh/MIC/MIM when no indices given + if method in ["cacoh", "mic", "mim"]: con = spectral_connectivity_epochs( data, method=method, mode=mode, indices=None, sfreq=sfreq ) assert con.indices is None assert con.n_nodes == n_signals - if method == "mic": + if method in ["cacoh", "mic"]: assert np.array(con.attrs["patterns"]).shape[2] == n_signals # check ragged indices padded correctly - ragged_indices = (np.array([[0]]), np.array([[1, 2]])) + ragged_indices = ([[0]], [[1, 2]]) con = spectral_connectivity_epochs( data, method=method, mode=mode, indices=ragged_indices, sfreq=sfreq ) - assert np.all( - np.array(con.indices) == np.array([np.array([[0, -1]]), np.array([[1, 2]])]) - ) + assert np.all(np.array(con.indices) == np.array([[[0, -1]], [[1, 2]]])) - # check shape of MIC patterns - if method == "mic": + # check shape of CaCoh/MIC patterns + if method in ["cacoh", "mic"]: for mode in ["multitaper", "cwt_morlet"]: con = spectral_connectivity_epochs( data, @@ -636,7 +634,7 @@ def test_spectral_connectivity_epochs_multivariate(method): assert np.shape(con.attrs["patterns"][1][0])[1] == len(fmin) # check patterns shape matches input data, not rank - rank = (np.array([1]), np.array([1])) + rank = ([1], [1]) con = spectral_connectivity_epochs( data, method=method, @@ -649,7 +647,7 @@ def test_spectral_connectivity_epochs_multivariate(method): assert np.shape(con.attrs["patterns"][1][0])[0] == n_targets # check patterns padded correctly - ragged_indices = (np.array([[0]]), np.array([[1, 2]])) + ragged_indices = ([[0]], [[1, 2]]) con = spectral_connectivity_epochs( data, method=method, mode=mode, indices=ragged_indices, sfreq=sfreq ) @@ -684,8 +682,8 @@ def test_multivariate_spectral_connectivity_epochs_regression(): fpath = os.path.dirname(os.path.realpath(__file__)) data = pd.read_pickle(os.path.join(fpath, "data", "example_multivariate_data.pkl")) sfreq = 100 - indices = (np.array([[0, 1]]), np.array([[2, 3]])) - methods = ["mic", "mim", "gc", "gc_tr"] + indices = ([[0, 1]], [[2, 3]]) + methods = ["cacoh", "mic", "mim", "gc", "gc_tr"] con = spectral_connectivity_epochs( data, method=methods, @@ -704,9 +702,15 @@ def test_multivariate_spectral_connectivity_epochs_regression(): n_jobs=1, ) - # should take the absolute of the MIC scores, as the MATLAB implementation - # returns the absolute values. - mne_results = {this_con.method: np.abs(this_con.get_data()) for this_con in con} + mne_results = {} + for this_con in con: + # must take the absolute of the MIC scores, as the MATLAB + # implementation returns the absolute values. + if this_con.method == "mic": + mne_results[this_con.method] = np.abs(this_con.get_data()) + else: + mne_results[this_con.method] = this_con.get_data() + matlab_results = pd.read_pickle( os.path.join(fpath, "data", "example_multivariate_matlab_results.pkl") ) @@ -715,7 +719,8 @@ def test_multivariate_spectral_connectivity_epochs_regression(): @pytest.mark.parametrize( - "method", ["mic", "mim", "gc", "gc_tr", ["mic", "mim", "gc", "gc_tr"]] + "method", + ["cacoh", "mic", "mim", "gc", "gc_tr", ["cacoh", "mic", "mim", "gc", "gc_tr"]], ) @pytest.mark.parametrize("mode", ["multitaper", "fourier", "cwt_morlet"]) def test_multivar_spectral_connectivity_epochs_error_catch(method, mode): @@ -731,14 +736,14 @@ def test_multivar_spectral_connectivity_epochs_error_catch(method, mode): rng_seed=0, ) - indices = (np.array([[0, 1]]), np.array([[2, 3]])) + indices = ([[0, 1]], [[2, 3]]) cwt_freqs = np.arange(10, 25 + 1) # check bad indices without nested array caught with pytest.raises( TypeError, match="multivariate indices must contain array-likes" ): - non_nested_indices = (np.array([0, 1]), np.array([2, 3])) + non_nested_indices = ([0, 1], [2, 3]) spectral_connectivity_epochs( data, method=method, mode=mode, indices=non_nested_indices, gc_n_lags=10 ) @@ -747,7 +752,7 @@ def test_multivar_spectral_connectivity_epochs_error_catch(method, mode): with pytest.raises( ValueError, match="multivariate indices cannot contain repeated" ): - repeated_indices = (np.array([[0, 1, 1]]), np.array([[2, 2, 3]])) + repeated_indices = ([[0, 1, 1]], [[2, 2, 3]]) spectral_connectivity_epochs( data, method=method, mode=mode, indices=repeated_indices, gc_n_lags=10 ) @@ -763,7 +768,7 @@ def test_multivar_spectral_connectivity_epochs_error_catch(method, mode): ) # check bad rank args caught - too_low_rank = (np.array([0]), np.array([0])) + too_low_rank = ([0], [0]) with pytest.raises(ValueError, match="ranks for seeds and targets must be"): spectral_connectivity_epochs( data, @@ -773,7 +778,7 @@ def test_multivar_spectral_connectivity_epochs_error_catch(method, mode): rank=too_low_rank, cwt_freqs=cwt_freqs, ) - too_high_rank = (np.array([3]), np.array([3])) + too_high_rank = ([3], [3]) with pytest.raises(ValueError, match="ranks for seeds and targets must be"): spectral_connectivity_epochs( data, @@ -793,7 +798,7 @@ def test_multivar_spectral_connectivity_epochs_error_catch(method, mode): rank=too_few_rank, cwt_freqs=cwt_freqs, ) - too_much_rank = (np.array([2, 2]), np.array([2, 2])) + too_much_rank = ([2, 2], [2, 2]) with pytest.raises(ValueError, match="rank argument must have shape"): spectral_connectivity_epochs( data, @@ -824,9 +829,9 @@ def test_multivar_spectral_connectivity_epochs_error_catch(method, mode): gc_n_lags=10, cwt_freqs=cwt_freqs, ) - assert rank_con.attrs["rank"] == (np.array([1]), np.array([1])) + assert rank_con.attrs["rank"] == ([1], [1]) - if method in ["mic", "mim"]: + if method in ["cacoh", "mic", "mim"]: # check rank-deficient transformation matrix caught with pytest.raises(RuntimeError, match="the transformation matrix"): spectral_connectivity_epochs( @@ -835,7 +840,7 @@ def test_multivar_spectral_connectivity_epochs_error_catch(method, mode): mode=mode, indices=indices, sfreq=sfreq, - rank=(np.array([2]), np.array([2])), + rank=([2], [2]), cwt_freqs=cwt_freqs, ) @@ -863,7 +868,7 @@ def test_multivar_spectral_connectivity_epochs_error_catch(method, mode): ) # check intersecting indices caught - bad_indices = (np.array([[0, 1]]), np.array([[0, 2]])) + bad_indices = ([[0, 1]], [[0, 2]]) with pytest.raises( ValueError, match="seed and target indices must not intersect" ): @@ -891,12 +896,12 @@ def test_multivar_spectral_connectivity_epochs_error_catch(method, mode): mode=mode, indices=indices, sfreq=sfreq, - rank=(np.array([2]), np.array([2])), + rank=([2], [2]), cwt_freqs=cwt_freqs, ) -@pytest.mark.parametrize("method", ["mic", "mim", "gc", "gc_tr"]) +@pytest.mark.parametrize("method", ["cacoh", "mic", "mim", "gc", "gc_tr"]) def test_multivar_spectral_connectivity_parallel(method): """Test multivar. freq.-domain connectivity methods run in parallel.""" data = make_signals_in_freq_bands( @@ -909,7 +914,7 @@ def test_multivar_spectral_connectivity_parallel(method): rng_seed=0, ) - indices = (np.array([[0, 1]]), np.array([[2, 3]])) + indices = ([[0, 1]], [[2, 3]]) spectral_connectivity_epochs( data, method=method, mode="multitaper", indices=indices, gc_n_lags=10, n_jobs=2 @@ -941,9 +946,9 @@ def test_multivar_spectral_connectivity_flipped_indices(): # if we're not careful, when finding the channels we need to compute the # CSD for, we might accidentally reorder the connectivity indices - indices = (np.array([[0, 1]]), np.array([[2, 3]])) - flipped_indices = (np.array([[2, 3]]), np.array([[0, 1]])) - concat_indices = (np.array([[0, 1], [2, 3]]), np.array([[2, 3], [0, 1]])) + indices = ([[0, 1]], [[2, 3]]) + flipped_indices = ([[2, 3]], [[0, 1]]) + concat_indices = ([[0, 1], [2, 3]], [[2, 3], [0, 1]]) # we test on GC since this is a directed connectivity measure method = "gc" @@ -1050,7 +1055,9 @@ def test_epochs_tmin_tmax(kind): assert len(w) == 1 # just one even though there were multiple epochs -@pytest.mark.parametrize("method", ["coh", "mic", "mim", "plv", "pli", "wpli", "ciplv"]) +@pytest.mark.parametrize( + "method", ["coh", "cacoh", "mic", "mim", "plv", "pli", "wpli", "ciplv"] +) @pytest.mark.parametrize("mode", ["cwt_morlet", "multitaper"]) @pytest.mark.parametrize("data_option", ["sync", "random"]) def test_spectral_connectivity_time_phaselocked(method, mode, data_option): @@ -1058,7 +1065,7 @@ def test_spectral_connectivity_time_phaselocked(method, mode, data_option): data.""" rng = np.random.default_rng(0) n_epochs = 5 - n_channels = 3 + n_channels = 4 n_times = 1000 sfreq = 250 data = np.zeros((n_epochs, n_channels, n_times)) @@ -1080,7 +1087,14 @@ def test_spectral_connectivity_time_phaselocked(method, mode, data_option): ) data[i, c] = np.squeeze(np.sin(x)) - multivar_methods = ["mic", "mim"] + multivar_methods = ["cacoh", "mic", "mim"] + + if method == "cacoh": + # CaCoh within set of signals will always be 1, so need to specify + # distinct seeds and targets + indices = ([[0, 1]], [[2, 3]]) + else: + indices = None # the frequency band should contain the frequency at which there is a # hypothesized "connection" @@ -1090,22 +1104,23 @@ def test_spectral_connectivity_time_phaselocked(method, mode, data_option): con = spectral_connectivity_time( data, freqs, + indices=indices, method=method, mode=mode, sfreq=sfreq, fmin=freq_band_low_limit, fmax=freq_band_high_limit, n_jobs=1, - faverage=True if method != "mic" else False, - average=True if method != "mic" else False, + faverage=method not in ["cacoh", "mic"], + average=method not in ["cacoh", "mic"], sm_times=0, ) con_matrix = con.get_data() - # MIC values can be pos. and neg., so must be averaged after taking the - # absolute values for the test to work + # CaCoh/MIC values can be pos. and neg., so must be averaged after taking + # the absolute values for the test to work if method in multivar_methods: - if method == "mic": + if method in ["cacoh", "mic"]: con_matrix = np.mean(np.abs(con_matrix), axis=(0, 2)) assert con.shape == (n_epochs, 1, len(con.freqs)) else: @@ -1145,7 +1160,7 @@ def test_spectral_connectivity_time_delayed(): trans_bandwidth = 2.0 # Hz delay = 5 # samples (non-zero delay needed for GC to be >> 0) - indices = (np.array([[0, 1]]), np.array([[2, 3]])) + indices = ([[0, 1]], [[2, 3]]) # 20-30 Hz connectivity fstart, fend = 20.0, 30.0 @@ -1409,7 +1424,7 @@ def test_spectral_connectivity_time_padding(method, mode, padding): ) -@pytest.mark.parametrize("method", ["mic", "mim", "gc", "gc_tr"]) +@pytest.mark.parametrize("method", ["cacoh", "mic", "mim", "gc", "gc_tr"]) @pytest.mark.parametrize("average", [True, False]) @pytest.mark.parametrize("faverage", [True, False]) def test_multivar_spectral_connectivity_time_shapes(method, average, faverage): @@ -1425,7 +1440,7 @@ def test_multivar_spectral_connectivity_time_shapes(method, average, faverage): rng_seed=0, ) - indices = (np.array([[0, 1]]), np.array([[2, 3]])) + indices = ([[0, 1]], [[2, 3]]) n_cons = len(indices[0]) freqs = np.arange(10, 25 + 1) @@ -1449,13 +1464,13 @@ def test_multivar_spectral_connectivity_time_shapes(method, average, faverage): ) assert con.shape == tuple(con_shape) - # check shape of MIC patterns are correct - if method == "mic": + # check shape of CaCoh/MIC patterns are correct + if method in ["cacoh", "mic"]: for indices_type in ["full", "ragged"]: if indices_type == "full": - indices = (np.array([[0, 1]]), np.array([[2, 3]])) + indices = ([[0, 1]], [[2, 3]]) else: - indices = (np.array([[0, 1]]), np.array([[2]])) + indices = ([[0, 1]], [[2]]) max_n_chans = 2 patterns_shape = [n_cons, max_n_chans] if faverage: @@ -1483,13 +1498,10 @@ def test_multivar_spectral_connectivity_time_shapes(method, average, faverage): assert not np.any(np.isnan(patterns[0, ..., :, :])) assert not np.any(np.isnan(patterns[0, ..., 0, :])) assert np.all(np.isnan(patterns[1, ..., 1, :])) # padded entry - assert np.all( - np.array(con.indices) - == np.array((np.array([[0, 1]]), np.array([[2, -1]]))) - ) + assert np.all(np.array(con.indices) == np.array(([[0, 1]], [[2, -1]]))) -@pytest.mark.parametrize("method", ["mic", "mim", "gc", "gc_tr"]) +@pytest.mark.parametrize("method", ["cacoh", "mic", "mim", "gc", "gc_tr"]) @pytest.mark.parametrize("mode", ["multitaper", "cwt_morlet"]) def test_multivar_spectral_connectivity_time_error_catch(method, mode): """Test error catching for time-resolved multivar. connectivity methods.""" @@ -1506,7 +1518,7 @@ def test_multivar_spectral_connectivity_time_error_catch(method, mode): rng_seed=0, ) - indices = (np.array([[0, 1]]), np.array([[2, 3]])) + indices = ([[0, 1]], [[2, 3]]) freqs = np.arange(10, 25 + 1) # test type-checking of data @@ -1517,7 +1529,7 @@ def test_multivar_spectral_connectivity_time_error_catch(method, mode): with pytest.raises( TypeError, match="multivariate indices must contain array-likes" ): - non_nested_indices = (np.array([0, 1]), np.array([2, 3])) + non_nested_indices = ([0, 1], [2, 3]) spectral_connectivity_time( data, freqs, method=method, mode=mode, indices=non_nested_indices ) @@ -1526,7 +1538,7 @@ def test_multivar_spectral_connectivity_time_error_catch(method, mode): with pytest.raises( ValueError, match="multivariate indices cannot contain repeated" ): - repeated_indices = (np.array([[0, 1, 1]]), np.array([[2, 2, 3]])) + repeated_indices = ([[0, 1, 1]], [[2, 2, 3]]) spectral_connectivity_time( data, freqs, method=method, mode=mode, indices=repeated_indices ) @@ -1539,12 +1551,12 @@ def test_multivar_spectral_connectivity_time_error_catch(method, mode): ) # check bad rank args caught - too_low_rank = (np.array([0]), np.array([0])) + too_low_rank = ([0], [0]) with pytest.raises(ValueError, match="ranks for seeds and targets must be"): spectral_connectivity_time( data, freqs, method=method, indices=indices, mode=mode, rank=too_low_rank ) - too_high_rank = (np.array([3]), np.array([3])) + too_high_rank = ([3], [3]) with pytest.raises(ValueError, match="ranks for seeds and targets must be"): spectral_connectivity_time( data, freqs, method=method, indices=indices, mode=mode, rank=too_high_rank @@ -1554,20 +1566,20 @@ def test_multivar_spectral_connectivity_time_error_catch(method, mode): spectral_connectivity_time( data, freqs, method=method, indices=indices, mode=mode, rank=too_few_rank ) - too_much_rank = (np.array([2, 2]), np.array([2, 2])) + too_much_rank = ([2, 2], [2, 2]) with pytest.raises(ValueError, match="rank argument must have shape"): spectral_connectivity_time( data, freqs, method=method, indices=indices, mode=mode, rank=too_much_rank ) - # check all-to-all conn. computed for MIC/MIM when no indices given - if method in ["mic", "mim"]: + # check all-to-all conn. computed for CaCoh/MIC/MIM when no indices given + if method in ["cacoh", "mic", "mim"]: con = spectral_connectivity_time( data, freqs, method=method, indices=None, mode=mode ) assert con.indices is None assert con.n_nodes == n_signals - if method == "mic": + if method in ["cacoh", "mic"]: assert np.array(con.attrs["patterns"]).shape[3] == n_signals if method in ["gc", "gc_tr"]: @@ -1578,7 +1590,7 @@ def test_multivar_spectral_connectivity_time_error_catch(method, mode): ) # check intersecting indices caught - bad_indices = (np.array([[0, 1]]), np.array([[0, 2]])) + bad_indices = ([[0, 1]], [[0, 2]]) with pytest.raises( ValueError, match="seed and target indices must not intersect" ): @@ -1631,12 +1643,12 @@ def test_multivar_save_load(tmp_path): tmp_file = os.path.join(tmp_path, "foo_mvc.nc") - non_ragged_indices = (np.array([[0, 1]]), np.array([[2, 3]])) - ragged_indices = (np.array([[0, 1]]), np.array([[2]])) + non_ragged_indices = ([[0, 1]], [[2, 3]]) + ragged_indices = ([[0, 1]], [[2]]) for indices in [non_ragged_indices, ragged_indices]: con = spectral_connectivity_epochs( epochs, - method=["mic", "mim", "gc", "gc_tr"], + method=["cacoh", "mic", "mim", "gc", "gc_tr"], indices=indices, fmin=10, fmax=30, @@ -1657,7 +1669,7 @@ def test_multivar_save_load(tmp_path): @pytest.mark.parametrize("method", ["coh", "plv", "pli", "wpli", "ciplv"]) -@pytest.mark.parametrize("indices", [None, (np.array([0, 1]), np.array([2, 3]))]) +@pytest.mark.parametrize("indices", [None, ([0, 1], [2, 3])]) def test_spectral_connectivity_indices_roundtrip_io(tmp_path, method, indices): """Test that indices values and type is maintained after saving. @@ -1699,8 +1711,8 @@ def test_spectral_connectivity_indices_roundtrip_io(tmp_path, method, indices): assert con.indices is None and read_con.indices is None -@pytest.mark.parametrize("method", ["mic", "mim", "gc", "gc_tr"]) -@pytest.mark.parametrize("indices", [None, (np.array([[0, 1]]), np.array([[2, 3]]))]) +@pytest.mark.parametrize("method", ["cacoh", "mic", "mim", "gc", "gc_tr"]) +@pytest.mark.parametrize("indices", [None, ([[0, 1]], [[2, 3]])]) def test_multivar_spectral_connectivity_indices_roundtrip_io(tmp_path, method, indices): """Test that indices values and type is maintained after saving. diff --git a/mne_connectivity/spectral/time.py b/mne_connectivity/spectral/time.py index 58b49cb1..412f46af 100644 --- a/mne_connectivity/spectral/time.py +++ b/mne_connectivity/spectral/time.py @@ -70,10 +70,12 @@ def spectral_connectivity_time( Only the frequencies within the range specified by ``fmin`` and ``fmax`` are used. method : str | list of str - Connectivity measure(s) to compute. These can be ``['coh', 'mic', - 'mim', 'plv', 'ciplv', 'pli', 'wpli', 'gc', 'gc_tr']``. These are: + Connectivity measure(s) to compute. These can be ``['coh', 'cacoh', + 'mic', 'mim', 'plv', 'ciplv', 'pli', 'wpli', 'gc', 'gc_tr']``. These + are: * %(coh)s + * %(cacoh)s * %(mic)s * %(mim)s * %(plv)s @@ -83,8 +85,8 @@ def spectral_connectivity_time( * %(gc)s * %(gc_tr)s - Multivariate methods (``['mic', 'mim', 'gc', 'gc_tr]``) cannot be - called with the other methods. + Multivariate methods (``['cacoh', 'mic', 'mim', 'gc', 'gc_tr']``) + cannot be called with the other methods. average : bool Average connectivity scores over epochs. If ``True``, output will be an instance of :class:`SpectralConnectivity`, otherwise @@ -152,7 +154,7 @@ def spectral_connectivity_time( Two arrays with the rank to project the seed and target data to, respectively, using singular value decomposition. If `None`, the rank of the data is computed and projected to. Only used if ``method`` - contains any of ``['mic', 'mim', 'gc', 'gc_tr']``. + contains any of ``['cacoh', 'mic', 'mim', 'gc', 'gc_tr']``. decim : int To reduce memory usage, decimation factor after time-frequency decomposition. Returns ``tfr[…, ::decim]``. @@ -253,24 +255,42 @@ def spectral_connectivity_time( C = --------------------- sqrt(E[Sxx] * E[Syy]) + 'cacoh' : Canonical Coherency (CaCoh) :footcite:`VidaurreEtAl2019` + given by: + + :math:`\textrm{CaCoh}=\Large{\frac{\boldsymbol{a}^T\boldsymbol{D} + (\Phi)\boldsymbol{b}}{\sqrt{\boldsymbol{a}^T\boldsymbol{a} + \boldsymbol{b}^T\boldsymbol{b}}}}` + + where: :math:`\boldsymbol{D}(\Phi)` is the cross-spectral density + between seeds and targets transformed for a given phase angle + :math:`\Phi`; and :math:`\boldsymbol{a}` and :math:`\boldsymbol{b}` + are eigenvectors for the seeds and targets, such that + :math:`\boldsymbol{a}^T\boldsymbol{D}(\Phi)\boldsymbol{b}` + maximises coherency between the seeds and targets. Taking the + absolute value of the results gives maximised coherence. + 'mic' : Maximised Imaginary part of Coherency (MIC) :footcite:`EwaldEtAl2012` given by: - :math:`MIC=\Large{\frac{\boldsymbol{\alpha}^T \boldsymbol{E \beta}} - {\parallel\boldsymbol{\alpha}\parallel \parallel\boldsymbol{\beta} - \parallel}}` + :math:`\textrm{MIC}=\Large{\frac{\boldsymbol{\alpha}^T + \boldsymbol{E \beta}}{\parallel\boldsymbol{\alpha}\parallel + \parallel\boldsymbol{\beta}\parallel}}` where: :math:`\boldsymbol{E}` is the imaginary part of the transformed cross-spectral density between seeds and targets; and :math:`\boldsymbol{\alpha}` and :math:`\boldsymbol{\beta}` are eigenvectors for the seeds and targets, such that - :math:`\boldsymbol{\alpha}^T \boldsymbol{E \beta}` maximises - connectivity between the seeds and targets. + :math:`\boldsymbol{\alpha}^T \boldsymbol{E \beta}` maximises the + imaginary part of coherency between the seeds and targets. 'mim' : Multivariate Interaction Measure (MIM) :footcite:`EwaldEtAl2012` given by: - :math:`MIM=tr(\boldsymbol{EE}^T)` + :math:`\textrm{MIM}=tr(\boldsymbol{EE}^T)` + + where :math:`\boldsymbol{E}` is the imaginary part of the + transformed cross-spectral density between seeds and targets. 'plv' : Phase-Locking Value (PLV) :footcite:`LachauxEtAl1999` given by:: @@ -300,7 +320,7 @@ def spectral_connectivity_time( :math:`GC = ln\Large{(\frac{\lvert\boldsymbol{S}_{tt}\rvert}{\lvert \boldsymbol{S}_{tt}-\boldsymbol{H}_{ts}\boldsymbol{\Sigma}_{ss - \lvert t}\boldsymbol{H}_{ts}^*\rvert}})`, + \lvert t}\boldsymbol{H}_{ts}^*\rvert}})` where: :math:`s` and :math:`t` represent the seeds and targets, respectively; :math:`\boldsymbol{H}` is the spectral transfer @@ -541,7 +561,12 @@ def spectral_connectivity_time( conn = dict() conn_patterns = dict() for m in method: - conn[m] = np.zeros((n_epochs, n_cons, n_freqs)) + # CaCoh complex-valued, all other methods real-valued + if m == "cacoh": + con_scores_dtype = np.complex128 + else: + con_scores_dtype = np.float64 + conn[m] = np.zeros((n_epochs, n_cons, n_freqs), dtype=con_scores_dtype) # prevent allocating memory for a huge array if not required if m in _patterns_methods: # patterns shape of [epochs x seeds/targets x cons x channels x freqs] diff --git a/mne_connectivity/utils/docs.py b/mne_connectivity/utils/docs.py index 3971ea72..fe56b1d9 100644 --- a/mne_connectivity/utils/docs.py +++ b/mne_connectivity/utils/docs.py @@ -70,7 +70,8 @@ docdict["coh"] = "'coh' : Coherence" docdict["cohy"] = "'cohy' : Coherency" -docdict["imcoh"] = "'imcoh' : Imaginary part of coherency" +docdict["imcoh"] = "'imcoh' : Imaginary part of Coherency" +docdict["cacoh"] = "'cacoh' : Canonical Coherency (CaCoh)" docdict["mic"] = "'mic' : Maximised Imaginary part of Coherency (MIC)" docdict["mim"] = "'mim' : Multivariate Interaction Measure (MIM)" docdict["plv"] = "'plv' : Phase-Locking Value (PLV)"