{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n# Compute multivariate measures of the imaginary part of coherency\n\nThis example demonstrates how multivariate methods based on the imaginary part\nof coherency :footcite:`EwaldEtAl2012` can be used to compute connectivity\nbetween whole sets of sensors, and how spatial patterns of this connectivity\ncan be interpreted.\n\nThe methods in question are: the maximised imaginary part of coherency (MIC);\nand the multivariate interaction measure (MIM; as well as its extension, the\nglobal interaction measure, GIM).\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Author: Thomas S. Binns \n# License: BSD (3-clause)\n# sphinx_gallery_thumbnail_number = 3"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import mne\nimport numpy as np\nfrom matplotlib import patheffects as pe\nfrom matplotlib import pyplot as plt\nfrom mne import EvokedArray, make_fixed_length_epochs\nfrom mne.datasets.fieldtrip_cmc import data_path\n\nfrom mne_connectivity import seed_target_indices, spectral_connectivity_epochs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Background\n\nMultivariate forms of signal analysis allow you to simultaneously consider\nthe activity of multiple signals. In the case of connectivity, the\ninteraction between multiple sensors can be analysed at once, producing a\nsingle connectivity spectrum. This approach brings not only practical\nbenefits (e.g. easier interpretability of results from the dimensionality\nreduction), but can also offer methodological improvements (e.g. enhanced\nsignal-to-noise ratio and reduced bias).\n\nA popular bivariate measure of connectivity is the imaginary part of\ncoherency, which looks at the correlation between two signals in the\nfrequency domain and is immune to spurious connectivity arising from volume\nconduction artefacts :footcite:`NolteEtAl2004`. However, in cases where\ninteractions between multiple signals are of interest, computing connectivity\nbetween all possible combinations of signals leads to a very large number of\nresults which is difficult to interpret. A common approach is to average\nresults across these connections, however this risks reducing the\nsignal-to-noise ratio of results and burying interactions that are present\nbetween only a small number of channels.\n\nAdditionally, this bivariate measure is susceptible to biased estimates of\nconnectivity based on the spatial proximity of sensors\n:footcite:`EwaldEtAl2012` depending on the degree of source mixing in the\nsignals.\n\nTo overcome these limitations, spatial filters derived from\neigendecompositions allow connectivity to be analysed in a multivariate\nmanner, removing the source mixing-dependent bias and increasing the\nsignal-to-noise ratio of connectivity estimates :footcite:`EwaldEtAl2012`.\nThis approach goes beyond simply aggregating information across all possible\ncombinations of signals, instead extracting the underlying components of\nconnectivity in a frequency-resolved manner.\n\nThis leads to the following methods: the maximised imaginary part of\ncoherency (MIC); and the multivariate interaction measure (MIM). These\nmethods are similar to the multivariate method based on coherency (CaCoh\n:footcite:`VidaurreEtAl2019`; see :doc:`cacoh` and\n:doc:`compare_coherency_methods`).\n\nWe start by loading some example MEG data and dividing it into\ntwo-second-long epochs.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"raw = mne.io.read_raw_ctf(data_path() / \"SubjectCMC.ds\")\nraw.pick(\"mag\")\nraw.crop(50.0, 110.0).load_data()\nraw.notch_filter(50)\nraw.resample(100)\n\nepochs = make_fixed_length_epochs(raw, duration=2.0).load_data()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We will focus on connectivity between sensors over the left and right\nhemispheres, with 75 sensors in the left hemisphere designated as seeds, and\n76 sensors in the right hemisphere designated as targets.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# left hemisphere sensors\nseeds = [idx for idx, ch_info in enumerate(epochs.info[\"chs\"]) if ch_info[\"loc\"][0] < 0]\n# right hemisphere sensors\ntargets = [\n idx for idx, ch_info in enumerate(epochs.info[\"chs\"]) if ch_info[\"loc\"][0] > 0\n]\n\nmultivar_indices = (np.array([seeds]), np.array([targets]))\n\nseed_names = [epochs.info[\"ch_names\"][idx] for idx in seeds]\ntarget_names = [epochs.info[\"ch_names\"][idx] for idx in targets]\n\n# multivariate imaginary part of coherency\n(mic, mim) = spectral_connectivity_epochs(\n epochs, method=[\"mic\", \"mim\"], indices=multivar_indices, fmin=5, fmax=30, rank=None\n)\n\n# bivariate imaginary part of coherency (for comparison)\nbivar_indices = seed_target_indices(seeds, targets)\nimcoh = spectral_connectivity_epochs(\n epochs, method=\"imcoh\", indices=bivar_indices, fmin=5, fmax=30\n)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"By averaging across each connection between the seeds and targets, we can see\nthat the bivariate measure of the imaginary part of coherency estimates a\nstrong peak in connectivity between seeds and targets around 13-18 Hz, with a\nweaker peak around 27 Hz.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig, axis = plt.subplots(1, 1)\naxis.plot(imcoh.freqs, np.mean(np.abs(imcoh.get_data()), axis=0), linewidth=2)\naxis.set_xlabel(\"Frequency (Hz)\")\naxis.set_ylabel(\"Absolute connectivity (A.U.)\")\nfig.suptitle(\"Imaginary part of coherency\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Maximised imaginary part of coherency (MIC)\n\nFor MIC, a set of spatial filters are found that will maximise the estimated\nconnectivity between the seed and target signals. These maximising filters\ncorrespond to the eigenvectors with the largest eigenvalue, derived from an\neigendecomposition of information from the cross-spectral density (Eq. 7 of\n:footcite:`EwaldEtAl2012`):\n\n$\\textrm{MIC}=\\Large{\\frac{\\boldsymbol{\\alpha}^T \\boldsymbol{E \\beta}}\n{\\parallel\\boldsymbol{\\alpha}\\parallel \\parallel\\boldsymbol{\\beta}\n\\parallel}}$,\n\nwhere $\\boldsymbol{\\alpha}$ and $\\boldsymbol{\\beta}$ are the\nspatial filters for the seeds and targets, respectively, and\n$\\boldsymbol{E}$ is the imaginary part of the transformed\ncross-spectral density between the seeds and targets. All elements are\nfrequency-dependent, however this is omitted for readability. MIC is bound\nbetween $[-1, 1]$ where the absolute value reflects connectivity\nstrength and the sign reflects the phase angle difference between signals.\n\nMIC can also be computed between identical sets of seeds and targets,\nallowing connectivity within a single set of signals to be estimated. This is\npossible as a result of the exclusion of zero phase lag components from the\nconnectivity estimates, which would otherwise return a perfect correlation.\n\nIn this instance, we see MIC reveal that in addition to the 13-18 Hz peak, a\npreviously unobserved peak in connectivity around 9 Hz is present.\nFurthermore, the previous peak around 27 Hz is much less pronounced. This may\nindicate that the connectivity was the result of some distal interaction\nexacerbated by strong source mixing, which biased the bivariate connectivity\nestimate.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig, axis = plt.subplots(1, 1)\naxis.plot(mic.freqs, np.abs(mic.get_data()[0]), linewidth=2)\naxis.set_xlabel(\"Frequency (Hz)\")\naxis.set_ylabel(\"Absolute connectivity (A.U.)\")\nfig.suptitle(\"Maximised imaginary part of coherency\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Furthermore, spatial patterns of connectivity can be constructed from the\nspatial filters to give a picture of the location of the channels involved in\nthe connectivity :footcite:`HaufeEtAl2014`. This information is stored under\n``attrs['patterns']`` of the connectivity class, with one value per frequency\nfor each channel in the seeds and targets. As with MIC, the absolute value of\nthe patterns reflect the strength, however the sign differences can be used\nto visualise the orientation of the underlying dipole sources. The spatial\npatterns are **not** bound between $[-1, 1]$.\n\nHere, we average across the patterns in the 13-18 Hz range. Plotting the\npatterns shows that the greatest connectivity between the left and right\nhemispheres occurs at the left and right posterior and left central regions,\nbased on the areas with the largest absolute values. Using the signs of the\nvalues, we can infer the existence of a dipole source between the central and\nposterior regions of the left hemisphere accounting for the connectivity\ncontributions (represented on the plot as a green line).\n\n**Note:** The spatial patterns are not a substitute for source\nreconstruction. If you need the spatial patterns in source space, you should\nperform source reconstruction before computing connectivity (see e.g.\n:doc:`mne_inverse_coherence_epochs`).\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# compute average of patterns in desired frequency range\nfband = [13, 18]\nfband_idx = [mic.freqs.index(freq) for freq in fband]\n\n# patterns have shape [seeds/targets x cons x channels x freqs (x times)]\npatterns = np.array(mic.attrs[\"patterns\"])\nseed_pattern = patterns[0, :, : len(seeds)]\ntarget_pattern = patterns[1, :, : len(targets)]\n# average across frequencies\nseed_pattern = np.mean(seed_pattern[0, :, fband_idx[0] : fband_idx[1] + 1], axis=1)\ntarget_pattern = np.mean(target_pattern[0, :, fband_idx[0] : fband_idx[1] + 1], axis=1)\n\n# store the patterns for plotting\nseed_info = epochs.copy().pick(seed_names).info\ntarget_info = epochs.copy().pick(target_names).info\nseed_pattern = EvokedArray(seed_pattern[:, np.newaxis], seed_info)\ntarget_pattern = EvokedArray(target_pattern[:, np.newaxis], target_info)\n\n# plot the patterns\nfig, axes = plt.subplots(1, 4)\nseed_pattern.plot_topomap(\n times=0,\n sensors=\"m.\",\n units=dict(mag=\"A.U.\"),\n cbar_fmt=\"%.1E\",\n axes=axes[0:2],\n time_format=\"\",\n show=False,\n)\ntarget_pattern.plot_topomap(\n times=0,\n sensors=\"m.\",\n units=dict(mag=\"A.U.\"),\n cbar_fmt=\"%.1E\",\n axes=axes[2:],\n time_format=\"\",\n show=False,\n)\naxes[0].set_position((0.1, 0.1, 0.35, 0.7))\naxes[1].set_position((0.4, 0.3, 0.02, 0.3))\naxes[2].set_position((0.5, 0.1, 0.35, 0.7))\naxes[3].set_position((0.9, 0.3, 0.02, 0.3))\naxes[0].set_title(\"Seed spatial pattern\\n13-18 Hz\")\naxes[2].set_title(\"Target spatial pattern\\n13-18 Hz\")\n\n# plot the left hemisphere dipole example\naxes[0].plot(\n [-0.01, -0.07],\n [-0.07, -0.03],\n color=\"lime\",\n linewidth=2,\n path_effects=[pe.Stroke(linewidth=4, foreground=\"k\"), pe.Normal()],\n)\n\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Multivariate interaction measure (MIM)\n\nAlthough it can be useful to analyse the single, largest connectivity\ncomponent with MIC, multiple such components exist and can be examined with\nMIM. MIM can be thought of as an average of all connectivity components\nbetween the seeds and targets, and can be useful for an exploration of all\navailable components. It is unnecessary to use the spatial filters of each\ncomponent explicitly, and instead the desired result can be achieved from\n$E$ alone (Eq. 14 of :footcite:`EwaldEtAl2012`):\n\n$\\textrm{MIM}=tr(\\boldsymbol{EE}^T)$,\n\nwhere again the frequency dependence is omitted.\n\nUnlike MIC, MIM is positive-valued and can be > 1. Without normalisation, MIM\ncan be thought of as reflecting the total interaction between the seeds and\ntargets. MIM can be normalised to lie in the range $[0, 1]$ by dividing\nthe scores by the number of unique channels in the seeds and targets.\nNormalised MIM represents the interaction *per channel*, which can be biased\nby factors such as the presence of channels with little to no interaction. In\nline with the preferences of the method's authors :footcite:`EwaldEtAl2012`,\nsince normalisation alters the interpretability of the results,\n**normalisation is not performed by default**.\n\nHere we see MIM reveal the strongest connectivity component to be around 10\nHz, with the higher frequency 13-18 Hz connectivity no longer being so\nprominent. This suggests that, across all components in the data, there may\nbe more lower frequency connectivity sources than higher frequency sources.\nThus, when combining these different components in MIM, the peak around 10 Hz\nremains, but the 13-18 Hz connectivity is diminished relative to the single,\nlargest connectivity component of MIC.\n\nLooking at the values for normalised MIM, we see it has a maximum of ~0.1.\nThe relatively small connectivity values thus indicate that many of the\nchannels show little to no interaction.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig, axis = plt.subplots(1, 1)\naxis.plot(mim.freqs, mim.get_data()[0], linewidth=2)\naxis.set_xlabel(\"Frequency (Hz)\")\naxis.set_ylabel(\"Absolute connectivity (A.U.)\")\nfig.suptitle(\"Multivariate interaction measure\")\n\nn_channels = len(seeds) + len(targets)\nnormalised_mim = mim.get_data()[0] / n_channels\nprint(f\"Normalised MIM has a maximum value of {normalised_mim.max():.2f}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Additionally, the instance where the seeds and targets are identical can be\nconsidered as a special case of MIM: the global interaction measure (GIM; Eq.\n15 of :footcite:`EwaldEtAl2012`). Again, this allows connectivity within a\nsingle set of signals to be estimated. Computing GIM follows from Eq. 14,\nhowever since each interaction is considered twice, correcting the\nconnectivity by a factor of $\\frac{1}{2}$ is necessary (**the\ncorrection is performed automatically in this implementation**). Like MIM,\nGIM can also be > 1, but it can again be normalised to lie in the range\n$[0, 1]$ by dividing by the number of unique channels in the seeds and\ntargets. However, since normalisation alters the interpretability of the\nresults (i.e. interaction per channel for normalised GIM vs. total\ninteraction for standard GIM), **GIM is not normalised by default**.\n\nWith GIM, we find a broad connectivity peak around 10 Hz, with an additional\npeak around 20 Hz. The differences observed with GIM highlight the presence\nof interactions within each hemisphere that are absent for MIC or MIM.\nFurthermore, the values for normalised GIM are higher than for MIM, with a\nmaximum of ~0.2, again indicating the presence of interactions across\nchannels within each hemisphere.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"indices = (np.array([[*seeds, *targets]]), np.array([[*seeds, *targets]]))\ngim = spectral_connectivity_epochs(\n epochs, method=\"mim\", indices=indices, fmin=5, fmax=30, rank=None, verbose=False\n)\n\nfig, axis = plt.subplots(1, 1)\naxis.plot(gim.freqs, gim.get_data()[0], linewidth=2)\naxis.set_xlabel(\"Frequency (Hz)\")\naxis.set_ylabel(\"Connectivity (A.U.)\")\nfig.suptitle(\"Global interaction measure\")\n\nn_channels = len(seeds) + len(targets)\nnormalised_gim = gim.get_data()[0] / n_channels\nprint(f\"Normalised GIM has a maximum value of {normalised_gim.max():.2f}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Handling high-dimensional data\n\nAn important issue to consider when using these multivariate methods is\noverfitting, which risks biasing connectivity estimates to maximise noise in\nthe data. This risk can be reduced by performing a preliminary dimensionality\nreduction prior to estimating the connectivity with a singular value\ndecomposition (Eqs. 32 & 33 of :footcite:`EwaldEtAl2012`). The degree of this\ndimensionality reduction can be specified using the ``rank`` argument, which\nby default will not perform any dimensionality reduction (assuming your data\nis full rank; see below if not). Choosing an expected rank of the data\nrequires *a priori* knowledge about the number of components you expect to\nobserve in the data.\n\nWhen comparing MIC/MIM scores across recordings, **it is highly recommended\nto estimate connectivity from the same number of channels (or equally from\nthe same degree of rank subspace projection)** to avoid biases in\nconnectivity estimates. Bias can be avoided by specifying a consistent rank\nsubspace to project to using the ``rank`` argument, standardising your\nconnectivity estimates regardless of changes in e.g. the number of channels\nacross recordings. Note that this does not refer to the number of seeds and\ntargets *within* a connection being identical, rather to the number of seeds\nand targets *across* connections.\n\nHere, we will project our seed and target data to only the first 25\ncomponents of our rank subspace. Results for MIM show that the general\nspectral pattern of connectivity is retained in the rank subspace-projected\ndata, suggesting that a fair degree of redundant connectivity information is\ncontained in the remaining 50 components of the seed and target data. We also\nassert that the spatial patterns of MIC are returned in the original sensor\nspace despite this rank subspace projection, being reconstructed using the\nproducts of the singular value decomposition (Eqs. 46 & 47 of\n:footcite:`EwaldEtAl2012`).\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"(mic_red, mim_red) = spectral_connectivity_epochs(\n epochs,\n method=[\"mic\", \"mim\"],\n indices=multivar_indices,\n fmin=5,\n fmax=30,\n rank=([25], [25]),\n)\n\n# subtract mean of scores for comparison\nmim_red_meansub = mim_red.get_data()[0] - mim_red.get_data()[0].mean()\nmim_meansub = mim.get_data()[0] - mim.get_data()[0].mean()\n\n# compare standard and rank subspace-projected MIM\nfig, axis = plt.subplots(1, 1)\naxis.plot(mim.freqs, mim_meansub, linewidth=2, label=\"standard MIM\")\naxis.plot(mim_red.freqs, mim_red_meansub, linewidth=2, label=\"rank subspace (25) MIM\")\naxis.set_xlabel(\"Frequency (Hz)\")\naxis.set_ylabel(\"Mean-corrected connectivity (A.U.)\")\naxis.legend()\nfig.suptitle(\"Multivariate interaction measure (non-normalised)\")\n\n# no. channels equal with and without projecting to rank subspace for patterns\nassert patterns[0, 0].shape[0] == np.array(mic_red.attrs[\"patterns\"])[0, 0].shape[0]\nassert patterns[1, 0].shape[0] == np.array(mic_red.attrs[\"patterns\"])[1, 0].shape[0]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the case that your data is not full rank and ``rank`` is left as ``None``,\nan automatic rank computation is performed and an appropriate degree of\ndimensionality reduction will be enforced. The rank of the data is determined\nby computing the singular values of the data and finding those within a\nfactor of $1e^{-6}$ relative to the largest singular value.\n\nWhilst unlikely, there may be scenarios in which this threshold is too\nlenient. In these cases, you should inspect the singular values of your data\nto identify an appropriate degree of dimensionality reduction to perform,\nwhich you can then specify manually using the ``rank`` argument. The code\nbelow shows one possible approach for finding an appropriate rank of\nclose-to-singular data with a more conservative threshold.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# gets the singular values of the data\ns = np.linalg.svd(raw.get_data(), compute_uv=False)\n# finds how many singular values are 'close' to the largest singular value\nrank = np.count_nonzero(s >= s[0] * 1e-4) # 1e-4 is the 'closeness' criteria, which is\n# a hyper-parameter"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Limitations\n\nThese multivariate methods offer many benefits in the form of dimensionality\nreduction, signal-to-noise ratio improvements, and invariance to\nestimate-biasing source mixing; however, no method is perfect. Important\nconsiderations must be taken into account when choosing methods based on the\nimaginary part of coherency such as MIC or MIM versus those based on\ncoherency/coherence, such as CaCoh.\n\nIn short, if you want to examine connectivity between signals from the same\nmodality or from different modalities using a shared reference, you should\nconsider using MIC and MIM to avoid spurious connectivity estimates stemming\nfrom e.g. volume conduction artefacts.\n\nOn the other hand, if you want to examine connectivity between signals from\ndifferent modalities using different references, CaCoh is a more appropriate\nmethod than MIC/MIM. This is because volume conduction artefacts are of less\nconcern, and CaCoh does not risk biasing connectivity estimates towards\ninteractions with particular phase lags like MIC/MIM.\n\nThese scenarios are described in more detail in the\n:doc:`compare_coherency_methods` example.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## References\n.. footbibliography::\n\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.12"
}
},
"nbformat": 4,
"nbformat_minor": 0
}