Note

Go to the end to download the full example code

# Compute multivariate measures of the imaginary part of coherency#

This example demonstrates how multivariate methods based on the imaginary part of coherency [1] can be used to compute connectivity between whole sets of sensors, and how spatial patterns of this connectivity can be interpreted.

The methods in question are: the maximised imaginary part of coherency (MIC); and the multivariate interaction measure (MIM; as well as its extension, the global interaction measure, GIM).

```
# Author: Thomas S. Binns <t.s.binns@outlook.com>
# License: BSD (3-clause)
```

```
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import patheffects as pe
import mne
from mne import EvokedArray, make_fixed_length_epochs
from mne.datasets.fieldtrip_cmc import data_path
from mne_connectivity import 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 and reduced bias).

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 [2]. However, depending on the degree of source mixing, this measure is susceptible to biased estimates of connectivity based on the spatial proximity of sensors [1].

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 [1]. This approach leads to the following methods: the maximised imaginary part of coherency (MIC); and the multivariate interaction measure (MIM).

We start by loading some example MEG data and dividing it into two-second-long epochs.

```
raw = mne.io.read_raw_ctf(data_path() / "SubjectCMC.ds")
raw.pick("mag")
raw.crop(50.0, 110.0).load_data()
raw.notch_filter(50)
raw.resample(100)
epochs = make_fixed_length_epochs(raw, duration=2.0).load_data()
```

```
ds directory : /home/circleci/mne_data/MNE-fieldtrip_cmc-data/SubjectCMC.ds
res4 data read.
hc data read.
Separate EEG position data file not present.
Quaternion matching (desired vs. transformed):
0.33 78.32 0.00 mm <-> 0.33 78.32 0.00 mm (orig : -71.62 40.46 -256.48 mm) diff = 0.000 mm
-0.33 -78.32 -0.00 mm <-> -0.33 -78.32 -0.00 mm (orig : 39.27 -70.16 -258.60 mm) diff = 0.000 mm
114.65 0.00 -0.00 mm <-> 114.65 0.00 -0.00 mm (orig : 64.35 66.64 -262.01 mm) diff = 0.000 mm
Coordinate transformations established.
Polhemus data for 3 HPI coils added
Device coordinate locations for 3 HPI coils added
Picked positions of 4 EEG channels from channel info
4 EEG locations added to Polhemus data.
Measurement info composed.
Finding samples for /home/circleci/mne_data/MNE-fieldtrip_cmc-data/SubjectCMC.ds/SubjectCMC.meg4:
System clock channel is available, checking which samples are valid.
75 x 12000 = 911610 samples from 191 chs
390 samples omitted at the end
Current compensation grade : 0
Removing 5 compensators from info because not all compensation channels were picked.
Reading 0 ... 72000 = 0.000 ... 60.000 secs...
Filtering raw data in 1 contiguous segment
Setting up band-stop filter from 49 - 51 Hz
FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandstop filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 49.38
- Lower transition bandwidth: 0.50 Hz (-6 dB cutoff frequency: 49.12 Hz)
- Upper passband edge: 50.62 Hz
- Upper transition bandwidth: 0.50 Hz (-6 dB cutoff frequency: 50.88 Hz)
- Filter length: 7921 samples (6.601 s)
[Parallel(n_jobs=1)]: Done 17 tasks | elapsed: 0.1s
[Parallel(n_jobs=1)]: Done 71 tasks | elapsed: 0.2s
Not setting metadata
30 matching events found
No baseline correction applied
0 projection items activated
Using data from preloaded Raw for 30 events and 200 original time points ...
0 bad epochs dropped
```

We will focus on connectivity between sensors over the left and right hemispheres, with 75 sensors in the left hemisphere designated as seeds, and 76 sensors in the right hemisphere designated as targets.

```
# left hemisphere sensors
seeds = [idx for idx, ch_info in enumerate(epochs.info["chs"]) if ch_info["loc"][0] < 0]
# right hemisphere sensors
targets = [
idx for idx, ch_info in enumerate(epochs.info["chs"]) if ch_info["loc"][0] > 0
]
multivar_indices = (np.array([seeds]), np.array([targets]))
seed_names = [epochs.info["ch_names"][idx] for idx in seeds]
target_names = [epochs.info["ch_names"][idx] for idx in targets]
# multivariate imaginary part of coherency
(mic, mim) = spectral_connectivity_epochs(
epochs, method=["mic", "mim"], indices=multivar_indices, fmin=5, fmax=30, rank=None
)
# bivariate imaginary part of coherency (for comparison)
bivar_indices = seed_target_indices(seeds, targets)
imcoh = spectral_connectivity_epochs(
epochs, method="imcoh", indices=bivar_indices, fmin=5, fmax=30
)
```

```
Adding metadata with 3 columns
Connectivity computation...
computing connectivity for 1 connections
using t=0.000s..1.990s for estimation (200 points)
frequencies: 5.0Hz..30.0Hz (51 points)
Estimated data ranks:
connection 1 - seeds (75); targets (76)
Using multitaper spectrum estimation with 7 DPSS windows
the following metrics will be computed: MIC, MIM
computing cross-spectral density for epoch 1
computing cross-spectral density for epoch 2
computing cross-spectral density for epoch 3
computing cross-spectral density for epoch 4
computing cross-spectral density for epoch 5
computing cross-spectral density for epoch 6
computing cross-spectral density for epoch 7
computing cross-spectral density for epoch 8
computing cross-spectral density for epoch 9
computing cross-spectral density for epoch 10
computing cross-spectral density for epoch 11
computing cross-spectral density for epoch 12
computing cross-spectral density for epoch 13
computing cross-spectral density for epoch 14
computing cross-spectral density for epoch 15
computing cross-spectral density for epoch 16
computing cross-spectral density for epoch 17
computing cross-spectral density for epoch 18
computing cross-spectral density for epoch 19
computing cross-spectral density for epoch 20
computing cross-spectral density for epoch 21
computing cross-spectral density for epoch 22
computing cross-spectral density for epoch 23
computing cross-spectral density for epoch 24
computing cross-spectral density for epoch 25
computing cross-spectral density for epoch 26
computing cross-spectral density for epoch 27
computing cross-spectral density for epoch 28
computing cross-spectral density for epoch 29
computing cross-spectral density for epoch 30
Computing MIC for connection 1 of 1
0%| | frequency blocks : 0/51 [00:00<?, ?it/s]
2%|▏ | frequency blocks : 1/51 [00:00<00:00, 50.37it/s]
4%|▍ | frequency blocks : 2/51 [00:00<00:00, 53.28it/s]
6%|▌ | frequency blocks : 3/51 [00:00<00:01, 39.04it/s]
8%|▊ | frequency blocks : 4/51 [00:00<00:01, 42.57it/s]
10%|▉ | frequency blocks : 5/51 [00:00<00:01, 44.99it/s]
12%|█▏ | frequency blocks : 6/51 [00:00<00:01, 33.51it/s]
14%|█▎ | frequency blocks : 7/51 [00:00<00:01, 35.87it/s]
16%|█▌ | frequency blocks : 8/51 [00:00<00:01, 37.95it/s]
18%|█▊ | frequency blocks : 9/51 [00:00<00:01, 32.69it/s]
20%|█▉ | frequency blocks : 10/51 [00:00<00:01, 34.50it/s]
22%|██▏ | frequency blocks : 11/51 [00:00<00:01, 36.16it/s]
24%|██▎ | frequency blocks : 12/51 [00:00<00:01, 31.77it/s]
25%|██▌ | frequency blocks : 13/51 [00:00<00:01, 33.24it/s]
27%|██▋ | frequency blocks : 14/51 [00:00<00:01, 34.63it/s]
29%|██▉ | frequency blocks : 15/51 [00:00<00:01, 31.26it/s]
31%|███▏ | frequency blocks : 16/51 [00:00<00:01, 32.53it/s]
33%|███▎ | frequency blocks : 17/51 [00:00<00:01, 33.75it/s]
35%|███▌ | frequency blocks : 18/51 [00:00<00:00, 34.89it/s]
37%|███▋ | frequency blocks : 19/51 [00:00<00:00, 32.24it/s]
39%|███▉ | frequency blocks : 20/51 [00:00<00:00, 33.36it/s]
41%|████ | frequency blocks : 21/51 [00:00<00:00, 34.42it/s]
43%|████▎ | frequency blocks : 22/51 [00:00<00:00, 31.66it/s]
45%|████▌ | frequency blocks : 23/51 [00:00<00:00, 32.67it/s]
47%|████▋ | frequency blocks : 24/51 [00:00<00:00, 33.66it/s]
49%|████▉ | frequency blocks : 25/51 [00:00<00:00, 31.58it/s]
51%|█████ | frequency blocks : 26/51 [00:00<00:00, 32.54it/s]
53%|█████▎ | frequency blocks : 27/51 [00:00<00:00, 33.47it/s]
55%|█████▍ | frequency blocks : 28/51 [00:00<00:00, 31.28it/s]
57%|█████▋ | frequency blocks : 29/51 [00:00<00:00, 32.16it/s]
59%|█████▉ | frequency blocks : 30/51 [00:00<00:00, 33.06it/s]
61%|██████ | frequency blocks : 31/51 [00:00<00:00, 31.03it/s]
63%|██████▎ | frequency blocks : 32/51 [00:00<00:00, 31.89it/s]
65%|██████▍ | frequency blocks : 33/51 [00:00<00:00, 32.75it/s]
67%|██████▋ | frequency blocks : 34/51 [00:01<00:00, 33.61it/s]
69%|██████▊ | frequency blocks : 35/51 [00:01<00:00, 31.60it/s]
71%|███████ | frequency blocks : 36/51 [00:01<00:00, 32.43it/s]
73%|███████▎ | frequency blocks : 37/51 [00:01<00:00, 33.24it/s]
75%|███████▍ | frequency blocks : 38/51 [00:01<00:00, 31.27it/s]
76%|███████▋ | frequency blocks : 39/51 [00:01<00:00, 32.08it/s]
78%|███████▊ | frequency blocks : 40/51 [00:01<00:00, 32.84it/s]
80%|████████ | frequency blocks : 41/51 [00:01<00:00, 30.99it/s]
82%|████████▏ | frequency blocks : 42/51 [00:01<00:00, 31.78it/s]
84%|████████▍ | frequency blocks : 43/51 [00:01<00:00, 32.55it/s]
86%|████████▋ | frequency blocks : 44/51 [00:01<00:00, 30.78it/s]
88%|████████▊ | frequency blocks : 45/51 [00:01<00:00, 31.56it/s]
90%|█████████ | frequency blocks : 46/51 [00:01<00:00, 29.56it/s]
92%|█████████▏| frequency blocks : 47/51 [00:01<00:00, 30.31it/s]
94%|█████████▍| frequency blocks : 48/51 [00:01<00:00, 31.09it/s]
96%|█████████▌| frequency blocks : 49/51 [00:01<00:00, 31.85it/s]
98%|█████████▊| frequency blocks : 50/51 [00:01<00:00, 30.44it/s]
100%|██████████| frequency blocks : 51/51 [00:01<00:00, 31.18it/s]
100%|██████████| frequency blocks : 51/51 [00:01<00:00, 31.98it/s]
Computing MIM for connection 1 of 1
0%| | frequency blocks : 0/51 [00:00<?, ?it/s]
2%|▏ | frequency blocks : 1/51 [00:00<00:00, 54.75it/s]
4%|▍ | frequency blocks : 2/51 [00:00<00:00, 54.89it/s]
6%|▌ | frequency blocks : 3/51 [00:00<00:00, 55.16it/s]
8%|▊ | frequency blocks : 4/51 [00:00<00:01, 35.92it/s]
10%|▉ | frequency blocks : 5/51 [00:00<00:01, 38.96it/s]
12%|█▏ | frequency blocks : 6/51 [00:00<00:01, 40.13it/s]
14%|█▎ | frequency blocks : 7/51 [00:00<00:01, 29.76it/s]
16%|█▌ | frequency blocks : 8/51 [00:00<00:01, 31.96it/s]
18%|█▊ | frequency blocks : 9/51 [00:00<00:01, 28.15it/s]
20%|█▉ | frequency blocks : 10/51 [00:00<00:01, 30.01it/s]
22%|██▏ | frequency blocks : 11/51 [00:00<00:01, 31.72it/s]
24%|██▎ | frequency blocks : 12/51 [00:00<00:01, 29.08it/s]
25%|██▌ | frequency blocks : 13/51 [00:00<00:01, 30.58it/s]
27%|██▋ | frequency blocks : 14/51 [00:00<00:01, 32.02it/s]
29%|██▉ | frequency blocks : 15/51 [00:00<00:01, 29.25it/s]
31%|███▏ | frequency blocks : 16/51 [00:00<00:01, 30.53it/s]
33%|███▎ | frequency blocks : 17/51 [00:00<00:01, 31.77it/s]
35%|███▌ | frequency blocks : 18/51 [00:00<00:01, 29.39it/s]
37%|███▋ | frequency blocks : 19/51 [00:00<00:01, 30.53it/s]
39%|███▉ | frequency blocks : 20/51 [00:00<00:00, 31.66it/s]
41%|████ | frequency blocks : 21/51 [00:00<00:00, 32.78it/s]
43%|████▎ | frequency blocks : 22/51 [00:00<00:00, 30.72it/s]
45%|████▌ | frequency blocks : 23/51 [00:00<00:00, 31.79it/s]
47%|████▋ | frequency blocks : 24/51 [00:00<00:00, 32.81it/s]
49%|████▉ | frequency blocks : 25/51 [00:00<00:00, 30.59it/s]
51%|█████ | frequency blocks : 26/51 [00:00<00:00, 31.58it/s]
53%|█████▎ | frequency blocks : 27/51 [00:00<00:00, 32.55it/s]
55%|█████▍ | frequency blocks : 28/51 [00:00<00:00, 30.71it/s]
57%|█████▋ | frequency blocks : 29/51 [00:00<00:00, 30.99it/s]
59%|█████▉ | frequency blocks : 30/51 [00:00<00:00, 28.74it/s]
61%|██████ | frequency blocks : 31/51 [00:01<00:00, 29.34it/s]
63%|██████▎ | frequency blocks : 32/51 [00:01<00:00, 30.21it/s]
65%|██████▍ | frequency blocks : 33/51 [00:01<00:00, 28.87it/s]
67%|██████▋ | frequency blocks : 34/51 [00:01<00:00, 29.75it/s]
69%|██████▊ | frequency blocks : 35/51 [00:01<00:00, 30.59it/s]
71%|███████ | frequency blocks : 36/51 [00:01<00:00, 28.99it/s]
73%|███████▎ | frequency blocks : 37/51 [00:01<00:00, 29.84it/s]
75%|███████▍ | frequency blocks : 38/51 [00:01<00:00, 30.69it/s]
76%|███████▋ | frequency blocks : 39/51 [00:01<00:00, 31.53it/s]
78%|███████▊ | frequency blocks : 40/51 [00:01<00:00, 30.08it/s]
80%|████████ | frequency blocks : 41/51 [00:01<00:00, 30.90it/s]
82%|████████▏ | frequency blocks : 42/51 [00:01<00:00, 31.71it/s]
84%|████████▍ | frequency blocks : 43/51 [00:01<00:00, 30.07it/s]
86%|████████▋ | frequency blocks : 44/51 [00:01<00:00, 30.88it/s]
88%|████████▊ | frequency blocks : 45/51 [00:01<00:00, 31.68it/s]
90%|█████████ | frequency blocks : 46/51 [00:01<00:00, 30.29it/s]
92%|█████████▏| frequency blocks : 47/51 [00:01<00:00, 31.08it/s]
94%|█████████▍| frequency blocks : 48/51 [00:01<00:00, 31.85it/s]
96%|█████████▌| frequency blocks : 49/51 [00:01<00:00, 30.23it/s]
98%|█████████▊| frequency blocks : 50/51 [00:01<00:00, 31.00it/s]
100%|██████████| frequency blocks : 51/51 [00:01<00:00, 31.74it/s]
100%|██████████| frequency blocks : 51/51 [00:01<00:00, 31.20it/s]
[Connectivity computation done]
Replacing existing metadata with 3 columns
Connectivity computation...
computing connectivity for 5700 connections
using t=0.000s..1.990s for estimation (200 points)
frequencies: 5.0Hz..30.0Hz (51 points)
Using multitaper spectrum estimation with 7 DPSS windows
the following metrics will be computed: Imaginary Coherence
computing cross-spectral density for epoch 1
computing cross-spectral density for epoch 2
computing cross-spectral density for epoch 3
computing cross-spectral density for epoch 4
computing cross-spectral density for epoch 5
computing cross-spectral density for epoch 6
computing cross-spectral density for epoch 7
computing cross-spectral density for epoch 8
computing cross-spectral density for epoch 9
computing cross-spectral density for epoch 10
computing cross-spectral density for epoch 11
computing cross-spectral density for epoch 12
computing cross-spectral density for epoch 13
computing cross-spectral density for epoch 14
computing cross-spectral density for epoch 15
computing cross-spectral density for epoch 16
computing cross-spectral density for epoch 17
computing cross-spectral density for epoch 18
computing cross-spectral density for epoch 19
computing cross-spectral density for epoch 20
computing cross-spectral density for epoch 21
computing cross-spectral density for epoch 22
computing cross-spectral density for epoch 23
computing cross-spectral density for epoch 24
computing cross-spectral density for epoch 25
computing cross-spectral density for epoch 26
computing cross-spectral density for epoch 27
computing cross-spectral density for epoch 28
computing cross-spectral density for epoch 29
computing cross-spectral density for epoch 30
[Connectivity computation done]
```

By averaging across each connection between the seeds and targets, we can see that the bivariate measure of the imaginary part of coherency estimates a strong peak in connectivity between seeds and targets around 13-18 Hz, with a weaker peak around 27 Hz.

```
fig, axis = plt.subplots(1, 1)
axis.plot(imcoh.freqs, np.mean(np.abs(imcoh.get_data()), axis=0), linewidth=2)
axis.set_xlabel("Frequency (Hz)")
axis.set_ylabel("Absolute connectivity (A.U.)")
fig.suptitle("Imaginary part of coherency")
```

```
Text(0.5, 0.98, 'Imaginary part of coherency')
```

## Maximised imaginary part of coherency (MIC)#

For MIC, 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. 7 of [1]):

\(MIC=\frac{\boldsymbol{\alpha}^T \boldsymbol{E \beta}}{\parallel \boldsymbol{\alpha}\parallel \parallel\boldsymbol{\beta}\parallel}\),

where \(\boldsymbol{\alpha}\) and \(\boldsymbol{\beta}\) are the spatial filters for the seeds and targets, respectively, and \(\boldsymbol{E}\) is the imaginary part of the transformed cross-spectral density between the seeds and targets. All elements are frequency-dependent, however this is omitted for readability. MIC is bound between \([-1, 1]\) where the absolute value reflects connectivity strength and the sign reflects the phase angle difference between signals.

MIC can also be computed between identical sets of seeds and targets, allowing connectivity within a single set of signals to be estimated. This is possible as a result of the exclusion of zero phase lag components from the connectivity estimates, which would otherwise return a perfect correlation.

In this instance, we see MIC reveal that in addition to the 13-18 Hz peak, a previously unobserved peak in connectivity around 9 Hz is present. Furthermore, the previous peak around 27 Hz is much less pronounced. This may indicate that the connectivity was the result of some distal interaction exacerbated by strong source mixing, which biased the bivariate connectivity estimate.

```
fig, axis = plt.subplots(1, 1)
axis.plot(mic.freqs, np.abs(mic.get_data()[0]), linewidth=2)
axis.set_xlabel("Frequency (Hz)")
axis.set_ylabel("Absolute connectivity (A.U.)")
fig.suptitle("Maximised imaginary part of coherency")
```

```
Text(0.5, 0.98, 'Maximised imaginary part of coherency')
```

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 \([-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 hemispheres occurs at the left and right posterior and left central regions, based on the areas with the largest absolute values. Using the signs of the 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).

```
# compute average of patterns in desired frequency range
fband = [13, 18]
fband_idx = [mic.freqs.index(freq) for freq in fband]
# patterns have shape [seeds/targets x cons x channels x freqs (x times)]
patterns = np.array(mic.attrs["patterns"])
seed_pattern = patterns[0, :, : len(seeds)]
target_pattern = patterns[1, :, : len(targets)]
# average across frequencies
seed_pattern = np.mean(seed_pattern[0, :, fband_idx[0] : fband_idx[1] + 1], axis=1)
target_pattern = np.mean(target_pattern[0, :, fband_idx[0] : fband_idx[1] + 1], axis=1)
# store the patterns for plotting
seed_info = epochs.copy().pick(seed_names).info
target_info = epochs.copy().pick(target_names).info
seed_pattern = EvokedArray(seed_pattern[:, np.newaxis], seed_info)
target_pattern = EvokedArray(target_pattern[:, np.newaxis], target_info)
# plot the patterns
fig, axes = plt.subplots(1, 4)
seed_pattern.plot_topomap(
times=0,
sensors="m.",
units=dict(mag="A.U."),
cbar_fmt="%.1E",
axes=axes[0:2],
time_format="",
show=False,
)
target_pattern.plot_topomap(
times=0,
sensors="m.",
units=dict(mag="A.U."),
cbar_fmt="%.1E",
axes=axes[2:],
time_format="",
show=False,
)
axes[0].set_position((0.1, 0.1, 0.35, 0.7))
axes[1].set_position((0.4, 0.3, 0.02, 0.3))
axes[2].set_position((0.5, 0.1, 0.35, 0.7))
axes[3].set_position((0.9, 0.3, 0.02, 0.3))
axes[0].set_title("Seed spatial pattern\n13-18 Hz")
axes[2].set_title("Target spatial pattern\n13-18 Hz")
# plot the left hemisphere dipole example
axes[0].plot(
[-0.01, -0.07],
[-0.07, -0.03],
color="lime",
linewidth=2,
path_effects=[pe.Stroke(linewidth=4, foreground="k"), pe.Normal()],
)
plt.show()
```

## Multivariate interaction measure (MIM)#

Although it can be useful to analyse the single, largest connectivity component with MIC, multiple such components exist and can be examined with MIM. MIM can be thought of as an average of all connectivity components between the seeds and targets, and can be useful for an exploration of all available components. It is unnecessary to use the spatial filters of each component explicitly, and instead the desired result can be achieved from \(E\) alone (Eq. 14 of [1]):

\(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
thought of as reflecting the total interaction between the seeds and targets.
MIM can be normalised to lie in the range \([0, 1]\) by dividing the
scores by the number of unique channels in the seeds and targets. Normalised
MIM represents the interaction *per channel*, which can be biased by factors
such as the presence of channels with little to no interaction. In line with
the preferences of the method’s authors [1], since
normalisation alters the interpretability of the results, **normalisation is
not performed by default**.

Here we see MIM reveal the strongest connectivity component to be around 10 Hz, with the higher frequency 13-18 Hz connectivity no longer being so prominent. This suggests that, across all components in the data, there may be more lower frequency connectivity sources than higher frequency sources. Thus, when combining these different components in MIM, the peak around 10 Hz remains, but the 13-18 Hz connectivity is diminished relative to the single, largest connectivity component of MIC.

Looking at the values for normalised MIM, we see it has a maximum of ~0.1. The relatively small connectivity values thus indicate that many of the channels show little to no interaction.

```
fig, axis = plt.subplots(1, 1)
axis.plot(mim.freqs, mim.get_data()[0], linewidth=2)
axis.set_xlabel("Frequency (Hz)")
axis.set_ylabel("Absolute connectivity (A.U.)")
fig.suptitle("Multivariate interaction measure")
n_channels = len(seeds) + len(targets)
normalised_mim = mim.get_data()[0] / n_channels
print(f"Normalised MIM has a maximum value of {normalised_mim.max():.2f}")
```

```
Normalised MIM has a maximum value of 0.10
```

Additionally, the instance where the seeds and targets are identical can be
considered as a special case of MIM: the global interaction measure (GIM; Eq.
15 of [1]). Again, this allows connectivity within a
single set of signals to be estimated. Computing GIM follows from Eq. 14,
however since each interaction is considered twice, correcting the
connectivity by a factor of \(\frac{1}{2}\) is necessary (**the
correction is performed automatically in this implementation**). Like MIM,
GIM can also be > 1, but it can again be normalised to lie in the range
\([0, 1]\) by dividing by the number of unique channels in the seeds and
targets. However, since normalisation alters the interpretability of the
results (i.e. interaction per channel for normalised GIM vs. total
interaction for standard GIM), **GIM is not normalised by default**.

With GIM, we find a broad connectivity peak around 10 Hz, with an additional peak around 20 Hz. The differences observed with GIM highlight the presence of interactions within each hemisphere that are absent for MIC or MIM. Furthermore, the values for normalised GIM are higher than for MIM, with a maximum of ~0.2, again indicating the presence of interactions across channels within each hemisphere.

```
indices = (np.array([[*seeds, *targets]]), np.array([[*seeds, *targets]]))
gim = spectral_connectivity_epochs(
epochs, method="mim", indices=indices, fmin=5, fmax=30, rank=None, verbose=False
)
fig, axis = plt.subplots(1, 1)
axis.plot(gim.freqs, gim.get_data()[0], linewidth=2)
axis.set_xlabel("Frequency (Hz)")
axis.set_ylabel("Connectivity (A.U.)")
fig.suptitle("Global interaction measure")
n_channels = len(seeds) + len(targets)
normalised_gim = gim.get_data()[0] / n_channels
print(f"Normalised GIM has a maximum value of {normalised_gim.max():.2f}")
```

```
Normalised GIM has a maximum value of 0.19
```

## 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 (Eqs. 32 & 33 of [1]). 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 MIC/MIM 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 will project our seed and target data to only the first 25 components of our rank subspace. Results for MIM 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 remaining 50 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 [1]).

```
(mic_red, mim_red) = spectral_connectivity_epochs(
epochs,
method=["mic", "mim"],
indices=multivar_indices,
fmin=5,
fmax=30,
rank=([25], [25]),
)
# subtract mean of scores for comparison
mim_red_meansub = mim_red.get_data()[0] - mim_red.get_data()[0].mean()
mim_meansub = mim.get_data()[0] - mim.get_data()[0].mean()
# 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.set_xlabel("Frequency (Hz)")
axis.set_ylabel("Mean-corrected connectivity (A.U.)")
axis.legend()
fig.suptitle("Multivariate interaction measure (non-normalised)")
# no. channels equal with and without projecting to rank subspace for patterns
assert patterns[0, 0].shape[0] == np.array(mic_red.attrs["patterns"])[0, 0].shape[0]
assert patterns[1, 0].shape[0] == np.array(mic_red.attrs["patterns"])[1, 0].shape[0]
```

```
Replacing existing metadata with 3 columns
Connectivity computation...
computing connectivity for 1 connections
using t=0.000s..1.990s for estimation (200 points)
frequencies: 5.0Hz..30.0Hz (51 points)
Using multitaper spectrum estimation with 7 DPSS windows
the following metrics will be computed: MIC, MIM
computing cross-spectral density for epoch 1
computing cross-spectral density for epoch 2
computing cross-spectral density for epoch 3
computing cross-spectral density for epoch 4
computing cross-spectral density for epoch 5
computing cross-spectral density for epoch 6
computing cross-spectral density for epoch 7
computing cross-spectral density for epoch 8
computing cross-spectral density for epoch 9
computing cross-spectral density for epoch 10
computing cross-spectral density for epoch 11
computing cross-spectral density for epoch 12
computing cross-spectral density for epoch 13
computing cross-spectral density for epoch 14
computing cross-spectral density for epoch 15
computing cross-spectral density for epoch 16
computing cross-spectral density for epoch 17
computing cross-spectral density for epoch 18
computing cross-spectral density for epoch 19
computing cross-spectral density for epoch 20
computing cross-spectral density for epoch 21
computing cross-spectral density for epoch 22
computing cross-spectral density for epoch 23
computing cross-spectral density for epoch 24
computing cross-spectral density for epoch 25
computing cross-spectral density for epoch 26
computing cross-spectral density for epoch 27
computing cross-spectral density for epoch 28
computing cross-spectral density for epoch 29
computing cross-spectral density for epoch 30
Computing MIC for connection 1 of 1
0%| | frequency blocks : 0/51 [00:00<?, ?it/s]
6%|▌ | frequency blocks : 3/51 [00:00<00:00, 146.29it/s]
12%|█▏ | frequency blocks : 6/51 [00:00<00:00, 146.53it/s]
18%|█▊ | frequency blocks : 9/51 [00:00<00:00, 88.56it/s]
24%|██▎ | frequency blocks : 12/51 [00:00<00:00, 98.06it/s]
29%|██▉ | frequency blocks : 15/51 [00:00<00:00, 104.82it/s]
33%|███▎ | frequency blocks : 17/51 [00:00<00:00, 82.91it/s]
39%|███▉ | frequency blocks : 20/51 [00:00<00:00, 89.61it/s]
45%|████▌ | frequency blocks : 23/51 [00:00<00:00, 95.13it/s]
49%|████▉ | frequency blocks : 25/51 [00:00<00:00, 80.98it/s]
55%|█████▍ | frequency blocks : 28/51 [00:00<00:00, 85.53it/s]
61%|██████ | frequency blocks : 31/51 [00:00<00:00, 90.07it/s]
65%|██████▍ | frequency blocks : 33/51 [00:00<00:00, 79.69it/s]
71%|███████ | frequency blocks : 36/51 [00:00<00:00, 84.02it/s]
76%|███████▋ | frequency blocks : 39/51 [00:00<00:00, 88.02it/s]
80%|████████ | frequency blocks : 41/51 [00:00<00:00, 79.57it/s]
86%|████████▋ | frequency blocks : 44/51 [00:00<00:00, 83.24it/s]
92%|█████████▏| frequency blocks : 47/51 [00:00<00:00, 86.72it/s]
96%|█████████▌| frequency blocks : 49/51 [00:00<00:00, 79.18it/s]
100%|██████████| frequency blocks : 51/51 [00:00<00:00, 83.09it/s]
Computing MIM for connection 1 of 1
0%| | frequency blocks : 0/51 [00:00<?, ?it/s]
6%|▌ | frequency blocks : 3/51 [00:00<00:00, 146.83it/s]
12%|█▏ | frequency blocks : 6/51 [00:00<00:00, 147.31it/s]
14%|█▎ | frequency blocks : 7/51 [00:00<00:00, 85.27it/s]
20%|█▉ | frequency blocks : 10/51 [00:00<00:00, 98.42it/s]
25%|██▌ | frequency blocks : 13/51 [00:00<00:00, 107.66it/s]
29%|██▉ | frequency blocks : 15/51 [00:00<00:00, 81.31it/s]
35%|███▌ | frequency blocks : 18/51 [00:00<00:00, 88.90it/s]
41%|████ | frequency blocks : 21/51 [00:00<00:00, 94.92it/s]
45%|████▌ | frequency blocks : 23/51 [00:00<00:00, 79.90it/s]
49%|████▉ | frequency blocks : 25/51 [00:00<00:00, 82.57it/s]
53%|█████▎ | frequency blocks : 27/51 [00:00<00:00, 85.06it/s]
57%|█████▋ | frequency blocks : 29/51 [00:00<00:00, 87.44it/s]
59%|█████▉ | frequency blocks : 30/51 [00:00<00:00, 74.86it/s]
65%|██████▍ | frequency blocks : 33/51 [00:00<00:00, 79.13it/s]
71%|███████ | frequency blocks : 36/51 [00:00<00:00, 83.12it/s]
73%|███████▎ | frequency blocks : 37/51 [00:00<00:00, 73.69it/s]
78%|███████▊ | frequency blocks : 40/51 [00:00<00:00, 77.84it/s]
84%|████████▍ | frequency blocks : 43/51 [00:00<00:00, 81.51it/s]
88%|████████▊ | frequency blocks : 45/51 [00:00<00:00, 74.09it/s]
94%|█████████▍| frequency blocks : 48/51 [00:00<00:00, 77.84it/s]
100%|██████████| frequency blocks : 51/51 [00:00<00:00, 81.42it/s]
100%|██████████| frequency blocks : 51/51 [00:00<00:00, 81.75it/s]
[Connectivity computation done]
```

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 \(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
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
```

## Limitations#

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.

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.

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 [3].

## References#

**Total running time of the script:** (0 minutes 33.700 seconds)