[FieldTrip] EEG MNE with standard_bem and statistics sanity

[hm ham]

Hello fellow fieldtrippers,

I sent these questions a few weeks back, but I think they got buried in the
weekend emails, so I'm asking the same things again. Sorry for the spam.

I'm using fieldtrip to do a MNE for EEG data with no subject MRIs and I
have a few questions to which I can't get a straight answer from tutorials
or from mailing list history. I guess this is sort of a sanity check.

First, I'm wondering about the correct sourcemodel, as loading the
standard_bem and setting the third mesh as the source points in
ft_prepare_leadfield gives an error: Warning: dipole lies on boundary of
volume model.

Code:
% Electrodes, vol, sourcemodel all set in same scale with ft_convert_units

vol = load('standard_bem');
cfg = [];
cfg.elec = elec_aligned;
   % electrodes aligned to skin surface
cfg.grid.pos = vol.bnd(3).pnt;              % source points
cfg.grid.inside = 1:size(vol.bnd(3).pnt,1);
cfg.vol = vol;
cfg.reducerank = 3;
leadfield = ft_prepare_leadfield(cfg);

Using the code below works:

vol = load('standard_bem');
sourcemodel = ft_read_headshape('cortex_8196.surf.gii');
cfg = [];
cfg.elec = elec_aligned;                 % electrodes aligned to skin
surface
cfg.grid.pos = sourcemodel;              % source points
cfg.grid.inside = 1:size(sourcemodel,1);
cfg.vol = vol;
cfg.reducerank = 3;
leadfield = ft_prepare_leadfield(cfg);

But as the http://www.fieldtriptoolbox.org/tutorial/headmodel_eeg seems to
suggest the third mesh in vol should be the brain? Or is the third
compartment there just for conduction calculations?

As using the sourcemodel 'cortex_8196.surf.gii' at least gave an output I
went forward with those for now, to see what other problems I'd encounter.
Running the MNE:

cfg        = [];
cfg.elec = elec_aligned;
cfg.method = 'mne';
cfg.grid   = leadfield;
cfg.vol    = vol;
cfg.mne.prewhiten = 'yes';
cfg.mne.lambda    = 60; %Hmmm?
cfg.mne.scalesourcecov = 'yes';
cfg.mne.normalize = 'yes';
cfg.channel = 'all';
cond_12_source{subject_id}  = ft_sourceanalysis(cfg,cond_12{subject_id});


So secondly, I'm concerned with the lambda value, as different tutorials
give very different advice ranging from 3 to 1e8. Running the MNE with 3
seems to give very unstable results when comparing the results between
subjects or conditions. The data is somewhat contaminated with leftover
eye-movement artefacts (removed with ICA). And that 3 is from the tutorial
where the data is from MEG. So a bigger lambda is in order, but how big? Or
just let the minimumnormestimate.m calculate it from the data? I have to
point out that the eye-movements are somewhat spread around the trial
durations and not present in the calculation window of the .cov field from
ft_timelockanalysis. Should I maybe calculate the .cov field from the whole
trial?

Third, there has been talk of MNE statistics in the mailing list, and it
seems that the ft_timelockstatistics will do the job if a neighbours
structure is provided. I build the neighbours structure from the
sourcemodel mesh like this:

sourcemodel = ft_read_headshape('cortex_8196.surf.gii');
nsources = length(sourcemodel.pnt);
neighbours=struct;
for i=1:nsources
    neighbours(i).label     = num2str(i);
    neighb_nodes            = sourcemodel.tri((sourcemodel.tri(:,1)==i |
sourcemodel.tri(:,2)==i | sourcemodel.tri(:,3)==i),:);
    neighb_nodes            = unique(neighb_nodes)';
    idx                     = (neighb_nodes ~= i);
    neighb_nodes            = neighb_nodes(idx);
    for k=1:length(neighb_nodes)
        neighbours(i).neighblabel{k} = num2str(neighb_nodes(k));
    end
end

This seemed to give reasonable neighbours and just going by distance gave
some sources that were on the other side of a gyrus. Although now with the
above method the distances between the sources are not uniform.

I also had to tweak the data structures a little bit, labels and such.
Then I ran the tests with:

cfg=[];
cfg.parameter           = 'avg.pow';
cfg.method              = 'montecarlo';
cfg.statistic           = 'depsamplesT';
cfg.parameter           = 'avg';
cfg.correctm            = 'cluster';
cfg.clusterstatistic    = 'maxsum';
cfg.minnbchan           = 2;
cfg.numrandomization    = 1000;
cfg.tail                = 1;
cfg.correcttail         ='alpha';
cfg.alpha               = 0.05;
cfg.latency             = [0.4 0.5];
cfg.avgovertime         = 'yes';

cfg.neighbours  = neighbours;

cfg.design = [1:15 1:15;ones(1,15), ones(1,15)*2];
cfg.uvar        = 1;
cfg.ivar        = 2;
cond1_stat = ft_timelockstatistics(cfg, cond_11_source{:},
cond_12_source{:});

Looking at the resulting clusters from ft_timelockstatistics there were
clusters that when plotted showed disconnected sources. For example,
plotted like this for the first cluster:

paint = ones(8196,3);
indx = find(cond1_stat.posclusterslabelmat == 1);
paint(indx,2) = 0;
ft_plot_mesh(sourcemodel,'vertexcolor',paint);


Gave the output in the attachment, in that case the third source is atleast
close, but still not a neighbour to the two next to it. Is there maybe
something I'm missing with using the ft_timelockstatistics like this?

I hope someone can clarify some of these issues.

Thank you already in advance!

Tatu
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20150415/89ea4ec1/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: disconnected_cluster.png
Type: image/png
Size: 83016 bytes
Desc: not available
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20150415/89ea4ec1/attachment-0001.png>