# Authors: Thomas S. Binns <t.s.binns@outlook.com>
#
# License: BSD (3-clause)
import numpy as np
from mne import BaseEpochs, EpochsArray
from mne.time_frequency import EpochsSpectrum, EpochsTFR
from mne.utils import _validate_type, warn
[docs]
def make_surrogate_data(data, n_shuffles=1000, rng_seed=None, return_generator=True):
"""Create surrogate data for a null hypothesis of connectivity.
Parameters
----------
data : ~mne.time_frequency.EpochsSpectrum | ~mne.time_frequency.EpochsTFR
The Fourier coefficients to create the null hypothesis surrogate data for. Can
be generated from :meth:`mne.Epochs.compute_psd` or
:meth:`mne.Epochs.compute_tfr` with ``output='complex'``.
.. note::
Storing Fourier coefficients in :class:`mne.time_frequency.EpochsSpectrum`
objects requires ``mne >= 1.8``.
n_shuffles : int (default 1000)
The number of surrogate datasets to create.
rng_seed : int | None (default None)
The seed to use for the random number generator. If ``None``, no seed is
specified.
return_generator : bool (default True)
Whether or not to return the surrogate data as a generator object instead of a
list. This allows iterating over the surrogates without having to keep them all
in memory.
Returns
-------
surrogate_data : list of ~mne.time_frequency.EpochsSpectrum or ~mne.time_frequency.EpochsTFR
The surrogate data for the null hypothesis with ``n_shuffles`` entries. Returned
as a generator if ``return_generator=True``.
Notes
-----
.. version-deprecated:: 0.9
This function was renamed to
:func:`~mne_connectivity.make_surrogate_resting_data`.
""" # noqa: E501
warn(
"`make_surrogate_data` is deprecated and will be removed in 1.0. Use "
"`make_surrogate_resting_data` instead.",
FutureWarning,
)
return make_surrogate_resting_data(
data=data,
n_shuffles=n_shuffles,
rng_seed=rng_seed,
return_generator=return_generator,
)
[docs]
def make_surrogate_resting_data(
data, n_shuffles=1000, rng_seed=None, return_generator=True
):
"""Create surrogate resting-state data for a null hypothesis of connectivity.
Parameters
----------
data : ~mne.Epochs | ~mne.time_frequency.EpochsSpectrum | ~mne.time_frequency.EpochsTFR
The time-series or Fourier coefficients to create the null hypothesis surrogate
data for. Fourier coefficients can be generated from
:meth:`mne.Epochs.compute_psd` or :meth:`mne.Epochs.compute_tfr` with
``output='complex'``.
.. version-changed:: 0.9
Added support for :class:`mne.Epochs` objects.
.. note::
Storing Fourier coefficients in :class:`mne.time_frequency.EpochsSpectrum`
objects requires ``mne >= 1.8``.
n_shuffles : int (default 1000)
The number of surrogate datasets to create.
rng_seed : int | None (default None)
The seed to use for the random number generator. If ``None``, no seed is
specified.
return_generator : bool (default True)
Whether or not to return the surrogate data as a generator object instead of a
list. This allows iterating over the surrogates without having to keep them all
in memory.
Returns
-------
surrogate_data : list of ~mne.EpochsArray or ~mne.time_frequency.EpochsSpectrum or ~mne.time_frequency.EpochsTFR
The surrogate data for the null hypothesis with ``n_shuffles`` entries. Returned
as a generator if ``return_generator=True``.
Notes
-----
.. admonition:: Use only with non-evoked data
This approach is only suitable for generating surrogate data for a null
hypothesis of connectivity from non-evoked data (e.g., resting state, inter-trial
period). For working with evoked data, see
:func:`~mne_connectivity.make_surrogate_evoked_data`.
Surrogate data is generated by randomly shuffling the order of epochs, independently
for each channel. This destroys the covariance of the data, such that connectivity
estimates should reflect the null hypothesis of no genuine connectivity between
signals (e.g., only interactions due to background noise)
:footcite:`PellegriniEtAl2023`.
For the surrogate data to properly reflect a null hypothesis, the data which is
shuffled **must not** have a temporal structure that is consistent across epochs.
Examples of this data include evoked potentials, where a stimulus is presented or an
action performed at a set time during each epoch. Such data should not be used for
generating surrogates, as even after shuffling the epochs, it will still show a high
degree of residual connectivity between channels. As a result, connectivity
estimates from your surrogate data will capture genuine interactions, instead of the
desired background noise. Treating these estimates as a null hypothesis will
increase the likelihood of a type II (false negative) error, i.e., that there is no
significant connectivity in your data.
Appropriate data for generating surrogates using this approach includes data from a
resting state, inter-trial period, or similar. Here, a strong temporal consistency
across epochs is not assumed, reducing the chances that genuine connectivity
information is captured in your surrogate connectivity estimates.
Finally, it is important to note that **you should always compare true and
surrogate connectivity estimates from epochs of the same duration**. This will
ensure that spectral information is captured with the same accuracy in both sets of
connectivity estimates. Ideally, **you should also compare true and surrogate
connectivity estimates from the same number of epochs** to avoid biases from noise
(fewer epochs gives noisier estimates) or finite sample sizes (e.g., in coherency,
phase-locking value, etc... :footcite:`VinckEtAl2010`).
.. note::
:func:`~mne_connectivity.make_surrogate_evoked_data` can also be used to
generate surrogates from non-evoked data. However, it is computationally more
expensive, especially if your end goal is to generate connectivity estimates
from non-time-resolved spectral coefficients (e.g., those in
:class:`~mne.time_frequency.EpochsSpectrum` objects).
.. version-changed:: 0.9
This function was renamed from ``make_surrogate_data`` to clarify that it is
intended for generating surrogate data from non-evoked (e.g., resting state)
data, in contrast to the new
:func:`~mne_connectivity.make_surrogate_evoked_data` function.
References
----------
.. footbibliography::
""" # noqa: E501
# Validate inputs
_validate_type(
data,
(BaseEpochs, EpochsSpectrum, EpochsTFR),
"data",
(
"mne.Epochs(Array), mne.time_frequency.EpochsSpectrum(Array), or "
"mne.time_frequency.EpochsTFR(Array)",
),
)
if not isinstance(data, BaseEpochs) and not np.iscomplexobj(data.get_data()):
raise TypeError(
"Values in `data` must be complex-valued when `data` is an EpochsSpectrum "
"or EpochsTFR object."
)
n_epochs, n_chans = data.get_data().shape[:2]
if n_epochs == 1:
raise ValueError("Data must contain more than one epoch for shuffling.")
if n_chans == 1:
raise ValueError("Data must contain more than one channel for shuffling.")
_validate_type(n_shuffles, "int-like", "n_shuffles", "int")
if n_shuffles < 1:
raise ValueError("Number of shuffles must be >= 1.")
_validate_type(return_generator, bool, "return_generator", "bool")
# rng_seed checked by NumPy later
# Make surrogate data and package into Epochs/EpochsSpectrum/EpochsTFR objects
surrogate_data = _shuffle_coefficients(
data, n_shuffles, rng_seed, data_form="resting"
)
if not return_generator:
surrogate_data = list(surrogate_data)
return surrogate_data
[docs]
def make_surrogate_evoked_data(
data, n_shuffles=1000, rng_seed=None, return_generator=True
):
"""Create surrogate evoked data for a null hypothesis of connectivity.
Parameters
----------
data : ~mne.Epochs | ~mne.time_frequency.EpochsTFR
The time-series or Fourier coefficients to create the null hypothesis surrogate
data for. Fourier coefficients can be generated from
:meth:`mne.Epochs.compute_tfr` with ``output='complex'``.
n_shuffles : int (default 1000)
The number of surrogate datasets to create.
rng_seed : int | None (default None)
The seed to use for the random number generator. If ``None``, no seed is
specified.
return_generator : bool (default True)
Whether or not to return the surrogate data as a generator object instead of a
list. This allows iterating over the surrogates without having to keep them all
in memory.
Returns
-------
surrogate_data : list of ~mne.EpochsArray or ~mne.time_frequency.EpochsTFR
The surrogate data for the null hypothesis with ``n_shuffles`` entries. Returned
as a generator if ``return_generator=True``.
Notes
-----
.. admonition:: Suitable for use with evoked or non-evoked data
This approach is suitable for generating surrogate data for a null hypothesis of
connectivity from evoked data. This approach can also work with non-evoked data,
however :func:`~mne_connectivity.make_surrogate_resting_data` is recommended for
this.
Surrogate data is generated by cutting the time series at a random point and
reversing the cut portion, independently for each epoch and channel. This destroys
the covariance of the data, such that connectivity estimates should reflect the null
hypothesis of no genuine connectivity between signals (e.g., only interactions due
to background noise) :footcite:`AruEtAl2015`.
This approach is robust to data containing a temporal structure that is consistent
across epochs. Examples of this data include evoked potentials, where a stimulus is
presented or an action performed at a set time during each epoch.
Finally, it is important to note that **you should always compare true and
surrogate connectivity estimates from epochs of the same duration**. This will
ensure that spectral information is captured with the same accuracy in both sets of
connectivity estimates. Ideally, **you should also compare true and surrogate
connectivity estimates from the same number of epochs** to avoid biases from noise
(fewer epochs gives noisier estimates) or finite sample sizes (e.g., in coherency,
phase-locking value, etc... :footcite:`VinckEtAl2010`).
.. versionadded:: 0.9
References
----------
.. footbibliography::
""" # noqa: E501
# Validate inputs
_validate_type(
data,
(BaseEpochs, EpochsTFR),
"data",
"mne.Epochs(Array) or mne.time_frequency.EpochsTFR(Array)",
)
if not isinstance(data, BaseEpochs) and not np.iscomplexobj(data.get_data()):
raise TypeError(
"Values in `data` must be complex-valued when `data` is an EpochsTFR "
"object."
)
n_chans, n_times = np.asarray(data.get_data().shape)[[1, -1]]
if n_chans == 1:
raise ValueError("Data must contain more than one channel for shuffling.")
if n_times == 1:
raise ValueError("Data must contain more than one timepoint for shuffling.")
_validate_type(n_shuffles, "int-like", "n_shuffles", "int")
if n_shuffles < 1:
raise ValueError("Number of shuffles must be >= 1.")
_validate_type(return_generator, bool, "return_generator", "bool")
# rng_seed checked by NumPy later
# Make surrogate data and package into Epochs/EpochsTFR objects
surrogate_data = _shuffle_coefficients(
data, n_shuffles, rng_seed, data_form="evoked"
)
if not return_generator:
surrogate_data = list(surrogate_data)
return surrogate_data
def _shuffle_coefficients(data, n_shuffles, rng_seed, data_form):
"""Shuffle coefficients to create surrogate data.
Surrogate data for each shuffle is packaged into an Epochs/EpochsSpectrum/EpochsTFR
object, which are together returned as a generator to minimise memory demand.
"""
# Extract data array and Epochs/Spectrum/TFR information
state = None
if isinstance(data, BaseEpochs):
data_class = EpochsArray
defaults = dict(
info=data.info,
events=data.events,
event_id=data.event_id,
tmin=data.tmin,
reject=data.reject,
flat=data.flat,
reject_tmin=data.reject_tmin,
reject_tmax=data.reject_tmax,
baseline=data.baseline,
proj=data.proj,
metadata=data.metadata,
selection=data.selection,
drop_log=data.drop_log,
raw_sfreq=data.info["sfreq"],
)
else:
defaults = dict(
method=None, tmin=None, tmax=None, picks=None, proj=None, n_jobs=None
)
if isinstance(data, EpochsSpectrum):
data_class = EpochsSpectrum
defaults.update(dict(fmin=None, fmax=None, exclude=(), remove_dc=None))
else:
data_class = EpochsTFR
defaults.update(dict(freqs=None, decim=None))
state = data.__getstate__()
data_arr = data.get_data()
# Make surrogate data
rng = np.random.default_rng(rng_seed)
if data_form == "resting":
yield from _shuffle_coefficients_over_epochs(
data_arr, state, defaults, data_class, n_shuffles, rng
)
else:
assert data_form == "evoked", (
f"Invalid `data_form`. Expected 'resting' or 'evoked', got {data_form}. "
"Please contact the MNE-Connectivity developers."
)
yield from _shuffle_coefficients_within_epochs(
data_arr, state, defaults, data_class, n_shuffles, rng
)
def _shuffle_coefficients_over_epochs(
data_arr, state, defaults, data_class, n_shuffles, rng
):
"""Shuffle coefficients over epochs to create surrogate data.
Intended for use with non-evoked (e.g., resting-state) data, as in Pellegrini et al.
(2023), DOI: 10.1016/j.neuroimage.2023.120218.
"""
for _ in range(n_shuffles):
# Shuffle epochs for each channel independently
surrogate_arr = np.zeros_like(data_arr, dtype=data_arr.dtype)
for chan_i in range(data_arr.shape[1]):
surrogate_arr[:, chan_i] = rng.permutation(data_arr[:, chan_i], axis=0)
# Package surrogate data for this shuffle (return as generator)
if data_class == EpochsArray:
yield data_class(data=surrogate_arr, **defaults)
else:
state["data"] = surrogate_arr
yield data_class(inst=state, **defaults)
def _shuffle_coefficients_within_epochs(
data_arr, state, defaults, data_class, n_shuffles, rng
):
"""Shuffle coefficients within epochs to create surrogate data.
Intended for use with (non-)evoked data, as in Aru et al. (2015), DOI:
10.1016/j.conb.2014.08.002.
"""
for _ in range(n_shuffles):
# Shuffle timepoints for each epoch and channel independently
surrogate_arr = np.zeros_like(data_arr, dtype=data_arr.dtype)
cutpoints = rng.integers(1, surrogate_arr.shape[-1], surrogate_arr.shape[:2])
for epoch_i in range(data_arr.shape[0]):
for chan_i in range(data_arr.shape[1]):
surr = np.array_split(
data_arr[epoch_i, chan_i], cutpoints[epoch_i, chan_i], axis=-1
)
surr.reverse()
surrogate_arr[epoch_i, chan_i] = np.concatenate(surr, axis=-1)
# Package surrogate data for this shuffle (return as generator)
if data_class == EpochsArray:
yield data_class(data=surrogate_arr, **defaults)
else:
state["data"] = surrogate_arr
yield data_class(inst=state, **defaults)