[FieldTrip] Set regularization parameter after SSS

[Vladimir Litvak]

Hi Luca and Matti,

We have recently looked into beamforming on Elekta with several
methods people including Robert and Alex Gramfort. One thing we
discovered is that beamforming can completely fail when small
regularisation (like 1% or 5% which most people use) is applied to
data after SSS. The theoretical underpinning for this has not been
figured out yet, but in practice, the more robust thing to do is to
reduce the data dimensionality and set lambda to zero. I attach an
example developed for the phantom data we got from Aston MEG that
shows how this can be done in FT. If you want, I can share the phantom
example with you so that you can verify that this approach works.

Regarding sensor type fusion, I would suggest that you only limit your
analysis to one sensor type. Most Neuromag people use the planar
gradiometers. The reason for this suggestion is that after SSS,
channels of both types are linear combinations of the same basis
vectors and contain redundant information so trying to fuse them is
more a headache than anything else. There is a paper saying this in
more detail http://www.mdpi.com/1424-8220/17/12/2926/html

Best,

Vladimir

On Mon, Aug 20, 2018 at 2:35 PM Luca Kaiser <luca.kaiser at web.de> wrote:
>
> Hi Matti,
>
> sure-sorry and thanks for your quick reply. So here is what I am doing using the average covariance matrix (so data covariance).
>
> cfg=[];
> cfg.covariance='yes';
> cfg.channel=data.label;
> avg_data=ft_timelockanalysis(cfg, data);
>
> cfg=[];
> cfg.method='lcmv';
> cfg.grid=lf;
> cfg.vol=hdm;
> cfg.lcmv.keepfilter='yes';
> cfg.lcmv.fixedori= 'yes';
> cfg.lcmv.lambda= '10%'; %0% rank deficient data-use stronger regularization??
> lcmv_avg=ft_sourceanalysis(cfg, avg_data);
>
>
> Best,
>
> Luca
>
>
> Gesendet: Montag, 20. August 2018 um 15:10 Uhr
> Von: "Matti Stenroos" <matti.stenroos at aalto.fi>
> An: "Luca Kaiser" <luca.kaiser at web.de>
> Betreff: Re: [FieldTrip] Set regularization parameter after SSS
> Dear Luca,
>
> I think you'd need to tell a bit more, as the meaning of "lambda"
> depends on the algorithm you are using. Also "the covariance matrix" is
> not uniquely defined --- are you talking about noise covariance or
> measurement covariance?
>
> In general, the eigenvalue idea you had read does not make sense --- the
> eigenvalue text deals with condition numbers, while lambda does, in
> general, not relate directly to that. The advice might have been that
> "set lambda so that the ratio of largest and smallest eigenvalues is
> 1000"...
>
> Cheers,
> Matti
> also working with low-rankingMax-filtered MEG and coincidentally
> having a lot of covariance metrics on screen at the moment
>
>
>
> On 2018-08-20 15:39, Luca Kaiser wrote:
> > Dear FieldTrip community,
> > I am using ft_sourceanalysis on SSS preprocessed (neuromag) data. I
> > wonder if there are any suggestions on how to set lambda in this case?
> > My covariance matrix is rank deficient (rank 60, 306 sensors). I read in
> > the mailing list that I would need to look at the eigenvalues of XX' and
> > set lambda to be 1/1000 of the largest eigenvalue. However, I do not
> > really understand why to divide by 1/1000?
> > Every help would be very much appreciated,
> > Thanks!
> > Luca
> >
> >
> >
> >
> >
> > _______________________________________________
> > fieldtrip mailing list
> > https://mailman.science.ru.nl/mailman/listinfo/fieldtrip
> > https://doi.org/10.1371/journal.pcbi.1002202
> >
> _______________________________________________
> fieldtrip mailing list
> https://mailman.science.ru.nl/mailman/listinfo/fieldtrip
> https://doi.org/10.1371/journal.pcbi.1002202
-------------- next part --------------
cd C:\home\Data\Aston\newphantom

ft_defaults;
%%

amp = 200;
dip = 5;
%%

headmodel = [];
headmodel.o = [0 0 0];
headmodel.r = 0.10;

%%
suffix = '';%'_tsss';% '';%******* Compare to ;
%%
%%
dataset = sprintf('.\\%dnAm\\Amp%d_Dip%d_IASoff%s.fif',amp, amp, dip, suffix);

hdr = ft_read_header(dataset);

hdr.chantype{strmatch('SYS201', hdr.label)} = 'other trigger';

ev = ft_read_event(dataset, 'header', hdr);

cfg = [];
cfg.dataset = dataset;
cfg.event = ev;
cfg.channel = {'megplanar', '-MEG1133', '-MEG1323'};
cfg.trialfun = 'ft_trialfun_general';
cfg.trialdef.prestim = 0.1;
cfg.trialdef.poststim = 0.1;
cfg.trialdef.eventtype = 'SYS201';
cfg.lpfilter = 'yes';
cfg.lpfreq = 40;
cfg.demean = 'yes';
cfg = ft_definetrial(cfg);
data = ft_preprocessing(cfg);
%%
cfg = [];
cfg.method = 'pca';
cfg.updatesens = 'no';
cfg.channel = 'megplanar';
comp = ft_componentanalysis(cfg, data);
 
cfg = [];
cfg.updatesens = 'no';
cfg.component = comp.label(51:end);
data = ft_rejectcomponent(cfg, comp);
%%
cfg = [];
cfg.toilim = [0 0.1];
active = ft_redefinetrial(cfg, data);

cfg = [];
cfg.toilim = [-0.1 0];
baseline = ft_redefinetrial(cfg, data);

%%
cfg = [];
cfg.covariance = 'yes';
cfg.covariancewindow = 'all';
timelock_active = ft_timelockanalysis(cfg, active);
timelock_baseline = ft_timelockanalysis(cfg, baseline);
timelock_all  =  ft_timelockanalysis(cfg, data);
%%
cfg = [];
cfg.headmodel = headmodel;
cfg.grid.resolution = 0.01;
sourcemodel = ft_prepare_sourcemodel(cfg, timelock_active);
%%
cfg  = [];
cfg.headmodel = headmodel;
cfg.grid = sourcemodel;
cfg.method = 'lcmv';
cfg.channel = 'megplanar';
cfg.lcmv.lambda = '0%';
cfg.lcmv.normalize = 'no';
cfg.lcmv.keepfilter = 'yes';
source_lcmv_filters   = ft_sourceanalysis(cfg, timelock_all);

%%
cfg.grid.filter = source_lcmv_filters.avg.filter;

source_lcmv_active   = ft_sourceanalysis(cfg, timelock_active);

source_lcmv_baseline = ft_sourceanalysis(cfg, timelock_baseline);


source_lcmv_active.time   = 0;
source_lcmv_baseline.time = 0;

cfg = [];
cfg.parameter = 'pow';
cfg.operation = 'log10(x1/x2)';
source_lcmv_relative = ft_math(cfg, source_lcmv_active, source_lcmv_baseline);


%%
cfg = [];
cfg.funparameter = 'pow';
cfg.location = 'max';
%%
ft_sourceplot(cfg, source_lcmv_relative);
%%