.. only:: html
.. note::
:class: sphx-glr-download-link-note
Click :ref:`here ` to download the full example code
.. rst-class:: sphx-glr-example-title
.. _sphx_glr_auto_examples_plot_surrogate_analysis.py:
Surrogate analysis
------------------
This example shows how to estimate a significance threshold in a comodulogram.
A comodulogram shows the estimated PAC metric on a grid of frequency bands.
In absence of PAC, a PAC metric will return values close to zero, but not
exactly zero. To estimate if a value is significantly far from zero, we use a
surrogate analysis.
In a surrogate analysis, we recompute the comodulogram many times, adding each
time a random time shift to remove any possible coupling between the
components. Nota that these time shifts have to be far from zero to effectively
remove a potential coupling. These comodulograms give us an estimation of the
fluctuation of the metric in the absence of PAC.
To derive a significance level from the list of comodulograms, we discuss here
two different methods:
- Computing a z-score on each couple of frequency, and thresholding at z = 4.
- Computing a threshold at a given p-value, over a distribution of comodulogram
maxima.
.. code-block:: default
import numpy as np
import matplotlib.pyplot as plt
from pactools import Comodulogram
from pactools import simulate_pac
Let's first create an artificial signal with PAC.
.. code-block:: default
fs = 200. # Hz
high_fq = 50.0 # Hz
low_fq = 5.0 # Hz
low_fq_width = 1.0 # Hz
n_points = 1000
noise_level = 0.4
signal = simulate_pac(n_points=n_points, fs=fs, high_fq=high_fq, low_fq=low_fq,
low_fq_width=low_fq_width, noise_level=noise_level,
random_state=0)
Then, let's define the range of low frequency, and the PAC metric used.
.. code-block:: default
low_fq_range = np.linspace(1, 10, 50)
method = 'duprelatour' # or 'tort', 'ozkurt', 'penny', 'colgin', ...
We also choose the number of comodulograms computed in the surrogate
analysis. A good rule of thumb is 10 / p_value. Example: 10 / 0.05 = 200.
.. code-block:: default
n_surrogates = 200
As a surrogate analysis recquires to compute many comodulograms, the
computation can be slow. If you have multiple cores in your CPU, you can
leverage them using the parameter `n_jobs` > 1.
.. code-block:: default
n_jobs = 1
To compute the comodulogram, we need to instanciate a `Comodulogram` object,
then call the method `fit`. Adding the surrogate analysis is as simple as
adding the `n_surrogates` parameter.
.. code-block:: default
estimator = Comodulogram(fs=fs, low_fq_range=low_fq_range,
low_fq_width=low_fq_width, method=method,
n_surrogates=n_surrogates, progress_bar=True,
n_jobs=n_jobs)
estimator.fit(signal)
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
[ ] 0% | 0.00 sec | comodulogram: DAR(10, 1)
[ ] 2% | 2.92 sec | comodulogram: DAR(10, 1)
[. ] 4% | 5.85 sec | comodulogram: DAR(10, 1)
[.. ] 6% | 8.78 sec | comodulogram: DAR(10, 1)
[... ] 8% | 11.76 sec | comodulogram: DAR(10, 1)
[.... ] 10% | 14.80 sec | comodulogram: DAR(10, 1)
[.... ] 12% | 17.78 sec | comodulogram: DAR(10, 1)
[..... ] 14% | 20.66 sec | comodulogram: DAR(10, 1)
[...... ] 16% | 23.62 sec | comodulogram: DAR(10, 1)
[....... ] 18% | 26.81 sec | comodulogram: DAR(10, 1)
[........ ] 20% | 30.32 sec | comodulogram: DAR(10, 1)
[........ ] 22% | 35.13 sec | comodulogram: DAR(10, 1)
[......... ] 24% | 40.81 sec | comodulogram: DAR(10, 1)
[.......... ] 26% | 44.83 sec | comodulogram: DAR(10, 1)
[........... ] 28% | 48.69 sec | comodulogram: DAR(10, 1)
[............ ] 30% | 51.91 sec | comodulogram: DAR(10, 1)
[............ ] 32% | 55.63 sec | comodulogram: DAR(10, 1)
[............. ] 34% | 58.93 sec | comodulogram: DAR(10, 1)
[.............. ] 36% | 62.19 sec | comodulogram: DAR(10, 1)
[............... ] 38% | 65.18 sec | comodulogram: DAR(10, 1)
[................ ] 40% | 68.08 sec | comodulogram: DAR(10, 1)
[................ ] 42% | 71.60 sec | comodulogram: DAR(10, 1)
[................. ] 44% | 75.72 sec | comodulogram: DAR(10, 1)
[.................. ] 46% | 79.59 sec | comodulogram: DAR(10, 1)
[................... ] 48% | 82.82 sec | comodulogram: DAR(10, 1)
[.................... ] 50% | 86.43 sec | comodulogram: DAR(10, 1)
[.................... ] 52% | 89.83 sec | comodulogram: DAR(10, 1)
[..................... ] 54% | 93.50 sec | comodulogram: DAR(10, 1)
[...................... ] 56% | 97.07 sec | comodulogram: DAR(10, 1)
[....................... ] 58% | 100.72 sec | comodulogram: DAR(10, 1)
[........................ ] 60% | 103.99 sec | comodulogram: DAR(10, 1)
[........................ ] 62% | 107.60 sec | comodulogram: DAR(10, 1)
[......................... ] 64% | 110.97 sec | comodulogram: DAR(10, 1)
[.......................... ] 66% | 114.53 sec | comodulogram: DAR(10, 1)
[........................... ] 68% | 118.17 sec | comodulogram: DAR(10, 1)
[............................ ] 70% | 121.28 sec | comodulogram: DAR(10, 1)
[............................ ] 72% | 125.06 sec | comodulogram: DAR(10, 1)
[............................. ] 74% | 128.19 sec | comodulogram: DAR(10, 1)
[.............................. ] 76% | 131.29 sec | comodulogram: DAR(10, 1)
[............................... ] 78% | 134.33 sec | comodulogram: DAR(10, 1)
[................................ ] 80% | 137.74 sec | comodulogram: DAR(10, 1)
[................................ ] 82% | 141.35 sec | comodulogram: DAR(10, 1)
[................................. ] 84% | 144.52 sec | comodulogram: DAR(10, 1)
[.................................. ] 86% | 147.83 sec | comodulogram: DAR(10, 1)
[................................... ] 88% | 151.08 sec | comodulogram: DAR(10, 1)
[.................................... ] 90% | 154.34 sec | comodulogram: DAR(10, 1)
[.................................... ] 92% | 157.52 sec | comodulogram: DAR(10, 1)
[..................................... ] 94% | 160.84 sec | comodulogram: DAR(10, 1)
[...................................... ] 96% | 164.17 sec | comodulogram: DAR(10, 1)
[....................................... ] 98% | 167.66 sec | comodulogram: DAR(10, 1)
[........................................] 100% | 170.84 sec | comodulogram: DAR(10, 1)
[........................................] 100% | 170.84 sec | comodulogram: DAR(10, 1)
Then we plot the significance level on top of the comodulogram.
Here we present two methods.
The z-score method presented here considers independently each pair of
frequency of the comodulogram. For each pair, we compute the mean `mu` and
standard deviation `sigma` of the PAC values computed over the surrogates
signals. We then transform the original PAC values `PAC` (non time-shifted)
into z-scores `Z`: Z = (PAC - mu) / sigma
This procedure, used for example in [Canolty et al, 2006], suffers from
multiple-testing issues, and also assumes that the distribution of PAC values
is Gaussian.
The other method presented here considers the ditribution of comodulogram
maxima. For each surrogate comodulogram, we compute the maximum PAC value.
From the obtained empirical distribution of maxima, we compute the
95-percentile, which corresponds to a p-value of 0.05.
This method does not assume the distribution to be Gaussian, nor suffers from
multiple-testing issues. This is the default method in the `plot` method.
.. code-block:: default
fig, axs = plt.subplots(1, 2, figsize=(10, 4))
z_score = 4.
estimator.plot(contour_method='z_score', contour_level=z_score,
titles=['With a z-score on each couple of frequency'],
axs=[axs[0]])
p_value = 0.05
estimator.plot(contour_method='comod_max', contour_level=p_value,
titles=['With a p-value on the distribution of maxima'],
axs=[axs[1]])
plt.show()
.. image:: /auto_examples/images/sphx_glr_plot_surrogate_analysis_001.png
:class: sphx-glr-single-img
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
/home/tom/work/github/pactools/examples/plot_surrogate_analysis.py:109: UserWarning: Matplotlib is currently using agg, which is a non-GUI backend, so cannot show the figure.
plt.show()
References
[Canolty et al, 2006]
Canolty, R. T., Edwards, E., Dalal, S. S., Soltani, M., Nagarajan,
S. S., Kirsch, H. E., ... & Knight, R. T. (2006). High gamma power is
phase-locked to theta oscillations in human neocortex. science,
313(5793), 1626-1628.
.. rst-class:: sphx-glr-timing
**Total running time of the script:** ( 2 minutes 51.227 seconds)
.. _sphx_glr_download_auto_examples_plot_surrogate_analysis.py:
.. only :: html
.. container:: sphx-glr-footer
:class: sphx-glr-footer-example
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download Python source code: plot_surrogate_analysis.py `
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download Jupyter notebook: plot_surrogate_analysis.ipynb `
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery `_