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

Wed Apr 1 05:22:59 EDT 2020

        External Email - Use Caution        

Dear Denis,

Thank you for your help!

Unfortunately, defining the rank as rank=69  does not work:

noise_covs2 = mne.compute_raw_covariance(
    raw_noise, method=methods,rank=69,
    return_estimators=True)

TypeError: rank must be an instance of None, dict, or str, got <class
'int'> instead

defining it as rank=dict( rank=69) or rank=dict(69)
does not help either.

Elena

On Tue, 31 Mar 2020 at 23:38, Denis A. Engemann <
denis-alexander.engemann at inria.fr> wrote:

>         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
>
>
> _______________________________________________
> Mne_analysis mailing list
> Mne_analysis at nmr.mgh.harvard.edu
> https://mail.nmr.mgh.harvard.edu/mailman/listinfo/mne_analysis

-- 
Best regards,
Elena V. Orekhova
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.nmr.mgh.harvard.edu/pipermail/mne_analysis/attachments/20200401/b0e6517a/attachment-0001.html 