Estimating the regularization parameter for HCP MEG resting state source reconstruction

Hi all,

I’m trying to figure out how to set the regularization parameter lambda2 for computing the inverse solution using eLORETA. I have processed the raw resting state MEG data from the Human Connectome Project according to the HCP pipeline and mne-hcp. In the inverse solution tutorial of mne-hcp a lambda2 of 1.0/9.0^2 is chosen without justification. In the MNE documentation the formula 1/SNR**2 is given as a reasonable estimate for lambda2, where SNR is the signal-to-noise ratio of the whitenend data. Since, I am interested in spontaneous resting state data, I think the signal-to-noise ratio should simply be the mean squared of the preprocessed data divided by the mean squared of the preprocessed empty room data. How is the whitened SNR computed? Is it the signal and noise whitened using the noise covariance of the sensor level empty room data? Are there any convenience functions in MNE that calculate this for you?

I would appreciate any help with this!

1 Like

Hello @Dominik and welcome to the forum!

I”m a bit in a hurry right now so I’ll have to keep it short, but what you can do is check out the estimate_snr() function. This will give you the SNR estimate for each time point. Please see the Notes section for the mathematical background.

As for how to collapse this into a single SNR measure, and how using different SNRs in the regularization parameter for source estimation, I don’t know.

I hope this helps a little bit,
Richard

Hi @richard,

I did check the estimate_snr() function and the math background, which made me wonder why MNE makes you choose a lambda2 if there seems to be a decent way to estimate the SNR?

Simply calculating the mean across timepoints of the SNR returned by estimate_snr() (which gave a warning that it did not converge) does give a reasonable value of ~4. I will simply use those estimates for calculating the source reconstruction.

thanks for your help,
Dominik

In fact, I discussed this with two other devs (@drammock and @britta-wstnr) yesterday, and none of us knew why :smiley: But I’m sure @agramfort or @larsoner can provide some insight here!

Best wishes,
Richard

1 Like

One thing that comes to mind is the time dimension: across which time points to take the average for a good estimate can be very different depending on your data/paradigm. There might not be a “one-fits-all” solution for this, i.e. one would have to query this in an extra parameter - which makes stuff convoluted.

1 Like

I agree with @britta-wstnr the time dimension makes it a bit tricky

1 Like

But now that @Dominik estimated an SNR of 4, should they use this value? Or resort to the more conservative 3 we commonly use in our examples?

Hello everyone,

I guess a lot of information regarding the important time episodes of a specific paradigm should already be contained in the definition of epochs. Another question is whether any reasonable computational estimate is better than human guesswork. If not making the SNR estimate the default for computing lambda2, it might be a great idea to make this an alternative option in the inverse solvers or at least mention this option in the docstring.

Concerning my specific experiment with empty room and resting-state data, I think it makes sense to use this estimate. However, it appears computationally demanding to do this for every subject and thus, I might use the default estimate of 3 for my study. I’ve seen this value of 3 recommended in some places but I couldn’t trace back to where this default originated from, can you tell me where this estimate came from?

Thanks for all your help and the great tools you provide!

1 Like

I believe we simply “assume” this in several inverse modeling tutorials that work with evoked data. There is at least one exception, I believe, where we’re working with single epochs (or even raw data??) and assume an SNR of 1. But I haven’t looked into the docs again before writing this response, so please take it with a grain of salt. :salt:

Have a great weekend!
Richard

If you ever come across the rationale behind using this default please let me know, I’d be curious where it comes from. Since I’m working with spontaneous resting state data, i.e. single-trial activity, I might check out this lower “default” of 1, which to me sounds much more realistic. Thanks for the discussion.

Marginally related, but are you a aware of this tutorial where an SNR “spectrum” is being derived, providing you with SNR estimates across frequency bins?

https://mne.tools/dev/auto_tutorials/time-freq/50_ssvep.html#plot-psd-and-snr-spectra

I found an older posting by @agramfort where he states when he defaults to an SNR of 1:

1 Like

This is very interesting and brings me to a related problem: If the SNR depends on the frequency band of interest, at what stage should one apply the respective narrowband filters? Let’s use the example of SSVEP. If I want to investigate the corresponding sources at the two distinct frequency bands, is it better to narrowband filter the activity around 12/15 Hz and use the respective frequency-specific SNR for the inverse operator but risk inducing temporal correlations in the noise-covariance violating its assumption of uncorrelated noise (as suggested in this old discussion). Or would you rather choose a common SNR of 1, reconstruct the sources, and narrowband filter the signals in source space? At least to me, it makes more sense to use the complete information of the broadband for the source reconstruction with a common SNR and then narrowband filter in the source space (which comes with its own problems of dipole orientations etc.). Any thoughts?

I honestly don’t know. Maybe @Denis has an informed opinion?

@Dominik

I just asked @agramfort about this while we were fighting with the coffee machine (spoiler: we won! :muscle: :coffee:)

He said this parameter (SNR of 3 for averaged, 1 for non-averaged data) “was established by Matti Hamalainen 20 years ago”, and that “no one touches this parameter”.

Hence … I would say, case closed. :wink:

Best wishes,
Richard

Haha, congrats for winning that struggle!

I already chose to blindly accept their authority as everyone else seemed to do and used SNR = 1. I would appreciate it if @agramfort could also point to the respective publication(s) for those estimates and then I’m done nagging you. :slight_smile: