Phase amplitude coupling (PAC) results differ between methods

Support & Discussions
1

Dear MNE Python community,

I want to compute phase-amplitude coupling using pactools. My problem is that I get different comodulograms when using different methods (Tort vs DAR model). I expected to see slight frequency differences but not such large differences. The DAR models seems also very sensitive to the choice of random seed, which is problematic as the output is not reliable.

Please find below a MWE ,which works with a data array (one selected MAG channel showing highest power) that I’m not sure how to share..

Thank a lot,

Best,

Laetitia

from pactools import Comodulogram
import numpy as np
import matplotlib.pyplot as plt

# PAC parameter
low_fq_range = np.arange(1,13,0.2)
high_fq_range = np.arange(13,125,2)
low_fq_width = 2.0  # Hz
high_fq_width = 20 
fs = 250 # Hz

# Read data from one selected channel (max overall power)
sig = np.load('sub-215_highestPowerChannel.npy') # n_epochs*timePoints

# Compute comodulogram for tort, DAR model 
method = 'tort'
estimator_tort = Comodulogram(fs=fs, low_fq_range=low_fq_range, high_fq_range=high_fq_range, random_state=42,low_fq_width=low_fq_width, high_fq_width = high_fq_width, method=method,progress_bar=True, n_jobs = 9)
estimator_tort.fit(sig)

method = 'duprelatour'
estimator_dlt = Comodulogram(fs=fs, low_fq_range=low_fq_range, high_fq_range=high_fq_range,  random_state=42,low_fq_width=low_fq_width, high_fq_width = high_fq_width, method=method,progress_bar=True, n_jobs = 9)
estimator_dlt.fit(sig)

# Plot
estimator_tort.plot(titles=['tort - seed 42'])
estimator_dlt.plot(titles=['duprelatour - seed 42']) # vmin = 0, vmax=0.05)

### The DAR model results heavily depends on which random seed you use
method = 'duprelatour'
seeds = [78, 34]
estimator_dlt_seeds = []
for seed in seeds:
    estimator = Comodulogram(fs=fs, low_fq_range=low_fq_range, high_fq_range=high_fq_range,  random_state=seed,low_fq_width=low_fq_width, high_fq_width = high_fq_width, method=method,progress_bar=True, n_jobs = 9)
    estimator.fit(sig)
    estimator.plot(titles=['duprelatour - seed {}'.format(seed)]) 

plt.show()
  • MNE version: 1.11.0
  • operating system: Windows 10