[FieldTrip] Granger Causality Group Analysis

[Harrison Fisher]

Hello,

I have 30+ subjects, each with 2 sessions, and 16 conditions for the data trials. I ultimately want to look at the difference between connectivity in one condition compared to another. 

For selected regions of interest, I created two virtual channels from the time locked dataset of all conditions concatenated and performed a mtmconvol analysis to get a wavelet (fr_frequencyanalysis, mtmconvol) for each of the virtual channels. Then I computed granger over the wavelet (ft_connectivityanalysis) and use ft_singleplotTFR to visualize. 

My main question is at what point do I need to separate out by condition in order to do a comparison of granger causality on the difference between two conditions? 

I am working off of a previous student’s scripts, and it appears that she was running the beamformer creation of the virtual channels on a trial by trial basis (see script excerpt below). 

Or is there a way to pull the conditions out of the wavelet or granger matrices created for the entire appended dataset of all the conditions? 

Thanks for the help!

Harrison Fisher
Bowdoin College 
Class of 2017

time  = timelock.All %select the appended dataset

       % apply LCMV spatial filter to location 01 of interest
       cfg = [];
       cfg.method = 'lcmv';
       cfg.lcmv.keepfilter = 'yes';
       cfg.headmodel = vol;
       cfg.elec = egi;
       cfg.grid.pos = roi1;
       source01 = ft_sourceanalysis(cfg, time);

       % construct 3-D virtual channel at location 01
       beamformer01 = source01.avg.filter{1};

       chansel = ft_channelselection('all' , data.label); % find the names
       chansel = match_str(data.label, chansel); % find the indices

       chan01_3D = [];
       chan01_3D.label = {'x', 'y', 'z'};
       chan01_3D.time = data.time;
       for i=1:length(data.trial)
           chan01_3D.trial{i} = beamformer01 * data.trial{i}(chansel,:);
       end

       % construct a single virtual channel in the maximum power orientation
       timeseries = cat(2, chan01_3D.trial{:});
       [u, s, v] = svd(timeseries, 'econ'); 
       timeseriesmaxproj = u(:, 1)' * timeseries; 

       chan01 = [];
       chan01.label = {'source01'};
       chan01.time = data.time;
       for i = 1:length(data.trial)
           chan01.trial{i} = u(:, 1)' * beamformer01 * data.trial{i}(chansel, :);
       end