[FieldTrip] Computing virtual channel ERFs within a source-level contrast's mask

[Barnett, Benjy]

Hey Guys,

I’m trying to make sense of some source-space analyses I’ve done and I wondered if I could get some expert input. I have computed a source-level contrast between condition A and B across subjects. Each subject performed both conditions. Specifically, I ran a one tailed test to reveal clusters where condition A showed significantly greater activation than condition B. My code to run this contrast is here:

%Group Stats
cfg=[];
cfg.dim         = grandavgA{1}.dim;
cfg.method      = 'montecarlo';
cfg.statistic   = 'ft_statfun_depsamplesT';
cfg.parameter   = 'pow';
cfg.correctm    = 'cluster';
cfg.numrandomization = 1000;
cfg.clusteralpha = 0.05;
cfg.alpha       = 0.05;
cfg.tail        = 1; %One tailed, see where A>B (I hope??)


nsubj=numel(grandavgZero);
cfg.design(1,:) = [1:nsubj 1:nsubj];
cfg.design(2,:) = [ones(1,nsubj)*1 ones(1,nsubj)*2];
cfg.uvar        = 1; % row of design matrix that contains unit variable (in this case: subjects)
cfg.ivar        = 2; % row of design matrix that contains independent variable (the conditions)


stat = ft_sourcestatistics(cfg,grandavgA{:}, grandavgB{:}); %where grandavgA etc. are cell arrays containing subject-level source reconstructions of condition A and B separately

This contrast reveals a large, significant occipito/parietal/temporal cluster where condition A has greater activity than B. I am then using this positive cluster mask as an ROI to create virtual channels and visualise the source-level ERFs for both condition A and B. I do this for each subject, warping the MNI coordinates of the cluster to individual subject-specific space, calculating virtual channels for each of the N pos within the cluster mask, giving me N*3 (multiplied by 3 for each XYZ direction of the dipole) virtual channels for each trial for each condition. The relevant code to do this per subject is here:

    gridpos= ft_warp_apply(M,stat.pos(stat.mask,:));%Calculate spatial filter for each voxel over all data

    %Calculate spatial filter for each voxel over all data
    cfg = [];
    cfg.method = 'lcmv';
    cfg.sourcemodel = sourcemodel;  %subject specific source model
    cfg.sourcemodel.pos = gridpos;
    cfg.sourcemodel.inside = true(length(cfg.sourcemodel.pos),1);
    cfg.unit = sourcemodel.unit;
    cfg.headmodel = headmodel;
    cfg.lcmv.keepfilter = 'yes';
    cfg.lcmv.lambda = '5%';
    cfg.channel = {'MEG'};
    cfg.senstype = 'MEG';
    sourceavg = ft_sourceanalysis(cfg, avg);
    dot_filter = cell2mat(sourceavg.avg.filter);

    chansel  = ft_channelselection('MEG', meg_data.label); % find MEG sensor names
    chanindx = match_str(meg_data.label, chansel);

    for c = 1:length(conditions)
        cond = conditions(c);
        if strcmp(cond,’A’)

            cfg = [];
            cfg.trials = meg_data.trialinfo(:,5) == c;
            data = ft_selectdata(cfg,dot_trials);
        elseif strcmp(cond,’B’)

            cfg = [];
            cfg.trials = meg_data.trialinfo(:,5) == c;
            data = ft_selectdata(cfg,dot_trials);
        end

        condvc = [];
        condvc.label = {'x' 'y' 'z'};
        condvc.time = data.time;
        for i = 1:length(data.trial)
            condvc.trial{i} = dot_filter * data.trial{i}(chanindx,:); %compute virtual channel
        end

        cfg = [];
        avgcond = ft_timelockanalysis(cfg,condvc);%compute subject level ERF for each condition

        VCs{subj,cond} = avgcond;


    end

The issue is that once I’ve taken the grand average of the virtual channel ERFs per each condition across subjects, the ERF for condition B (blue) is statistically significantly greater than condition A (red) across the timepoints the original source contrast was computed over. This seems completely backwards. Surely my process is essentially double-dipping, and I should see a greater ERF for condition A over condition B, given that I'm using virtual channels computed within a cluster where A was significantly greater than B at the whole brain level. I’ve tried it this way but also just using the peak voxel within the cluster based on the t-stat, but I get the same counterintuitive pattern of ERF results.

%% Compute Grand Avg over Subjects
cfg = [];
cfg.channel = 'all';
cfg.latency = 'all';
cfg.parameter = 'avg';
grandavg_A = ft_timelockgrandaverage(cfg,VCs{:,1});
grandavg_B = ft_timelockgrandaverage(cfg,VCs{:,2});

%% Plot Average Source-Level ERF For Zero vs non-zero
cfg = [];
cfg.channel = 'all';
%cfg.showlegend = 'yes';
cfg.linewidth = 1.75;
cfg.linecolor = [1,0,0; 0,0,1];
ft_singleplotER(cfg, grandavg_A, grandavg_B);


Am I thinking about this in the correct way, or have I made some error in my code? I’ve not copied my entire scripts so as not to clutter this question, but I hope I have supplied enough for you to follow my method.

All the best,
Benjy
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20231212/193de425/attachment.htm>