[Mne_analysis] 'shrunk' vs ‘empirical’ for the rest alpha power

[Denis A. Engemann]

        External Email - Use Caution        

Hi Elena,

thanks for having shared the data.
My suspicion got confirmed that something went wrong with the rank detection in your data, which sometimes happens with SSS’d data.

Please take a look this little script: https://gist.github.com/dengemann/8206bc00e82853c5b97e951259a7984f

You can see that for the first set of covariance matrices you get a warning informing you about failed rank detection.

As a result, thew shrunk estimator assumes 306 channels and leads to wrong results.

This can be easily seen by looking at the eigenvalue spectra sSecond figure): https://gist.github.com/dengemann/8206bc00e82853c5b97e951259a7984f#file-plot_debug_sss_cov_rank_2020_03_31-py-L21

It seems your rank is 69, so I set this manually and now both covariances are similar:

https://gist.github.com/dengemann/8206bc00e82853c5b97e951259a7984f#file-plot_debug_sss_cov_rank_2020_03_31-py-L32


I predict that your problem will disappear once you fix the rank computation.

Let me know how that goes.
Best wishes,
Denis


 	

> On Mar 31, 2020, at 5:13 PM, Denis A. Engemann <denis-alexander.engemann at inria.fr> wrote:
> 
>        External Email - Use Caution        
> 
> Hi Elena,
> 
> this is hard to tell by just looking at the code.
> 
> Would you be happy to share the (at least parts) of the raw & epochs file with me privately?
> 
> I could then investigate the issue.
> 
> Denis
> 
>> On Mar 31, 2020, at 4:38 PM, Elena Orekhova <orekhova.elena.v at gmail.com> wrote:
>> 
>>        External Email - Use Caution        
>> 
>> 
>> Hi,
>> I try to localize rest alpha band power using sLORETA.
>> 
>> I calculated noise covariance from the empty room data using:
>> 
>> noise_cov = mne.compute_raw_covariance(raw_noise, method=['shrunk', ‘empirical’])
>> 
>> The 'shrunk' was said to be a better choice. However, the result looks weird with the 'shrunk’ and much more normal with ‘empirical’. (This subject has prominent occipital alpha peak in sensors.)
>> 
>> The stc power plots for both methods are attached; the only parameter changed was method=['shrunk'] vs method= [‘empirical’]
>> 
>> Elena
>> 
>> p.s. 
>> The code is below.
>> I use :
>> 
>> Platform:      Darwin-18.7.0-x86_64-i386-64bit
>> Python:        3.6.8 |Anaconda, Inc.| (default, Dec 29 2018, 19:04:46)  [GCC 4.2.1 Compatible Clang 4.0.1 (tags/RELEASE_401/final)]
>> Executable:    /Users/elena/anaconda3/envs/mne/bin/python
>> CPU:           i386: 4 cores
>> Memory:        16.0 GB
>> 
>> mne:           0.20.0
>> numpy:         1.16.2 {blas=mkl_rt, lapack=mkl_rt}
>> scipy:         1.2.1
>> matplotlib:    3.0.3 {backend=Qt5Agg}
>> 
>> sklearn:       0.20.3
>> numba:         Not found
>> nibabel:       2.4.0
>> cupy:          Not found
>> pandas:        0.24.2
>> dipy:          0.16.0
>> mayavi:        4.7.0.dev0 {qt_api=pyqt5, PyQt5=5.9.2}
>> pyvista:       Not found
>> vtk:           8.1.2
>> 
>> ##########
>> epo_rest = mne.read_epochs(raw_fname , proj=False, verbose=None)
>> raw_noise = io.read_raw_fif(raw_fname, preload=True)
>> raw_noise.filter(2, 40, fir_design='firwin') 
>> noise_cov = mne.compute_raw_covariance(raw_noise, method=['shrunk', 'empirical'])
>> inverse_operator = make_inverse_operator(epo_rest.info, fwd, noise_cov, loose=0.2, depth=0.8, verbose=True)  
>> method = 'sLORETA' 
>> lambda2 = 1. / snr ** 2
>> bandwidth =  'hann'
>> stcs_rest = compute_source_psd_epochs(epo_rest[:n_epochs_use], inverse_operator,
>>                                 lambda2=lambda2,
>>                                 method=method, fmin=8, fmax=13,
>>                                 bandwidth=bandwidth,
>>                                 return_generator=True, verbose=True)
>> 
>>    psd_avg = 0.
>>    for i, stc_rest in enumerate(stcs_rest):
>>    psd_avg += stc_rest.data
>>    psd_avg /= n_epochs_use
>>    freqs = stc_rest.times  # the frequencies are stored here
>>    stc_rest.data = psd_avg  
>> 
>> clim=dict(kind='percent', lims=[50, 75, 100] ) 
>> kwargs = dict(initial_time=10.0,  clim=clim, hemi='both', subjects_dir=subjects_dir,time_unit='s',   size=(600, 600))
>> 
>> brain = stc_rest.plot(figure=1, **kwargs)
>> <empirical vs shrunk.pdf>_______________________________________________
>> Mne_analysis mailing list
>> Mne_analysis at nmr.mgh.harvard.edu
>> https://mail.nmr.mgh.harvard.edu/mailman/listinfo/mne_analysis
> 
> 
> _______________________________________________
> Mne_analysis mailing list
> Mne_analysis at nmr.mgh.harvard.edu
> https://mail.nmr.mgh.harvard.edu/mailman/listinfo/mne_analysis