[FieldTrip] Cluster-based permutation statistics using ft_timelockstatistics assistance

[Merlin Kelly]

Hello everyone
I have understood the correct fix to my issue, and I wanted to document
what I discovered just in case someone else encounters the same confusion
due to difficulty finding answers. The struct variable stat.stat
demonstrates the t-values of each element of the timelocked data.
I haven't been able to figure out how to run a three-way analysis using
ft_timelockstatistics, but instead, I use a for loop to run each
combination I desire to use.
The variable stat.prob demonstrates the location of each cluster with their
probability at that specific element, and the elements with value 1
demonstrate that no cluster was identified in that specific location.
Additionally, the error about unique permutations in a nutshell is I didn't
have enough subjects for 1000 unique permutations to occur.

There is still some confusion on my part on what other stat variable means,
but I hoped these brief descriptions above can help
Best
Merlin

On Tue, Sep 5, 2023 at 5:31 PM Merlin Kelly <merlinkelly75 at gmail.com> wrote:

> Hello
> I am currently trying to use the command ft_timelockstatistics to
> replicate the method of data analysis as demonstrated in the paper: "EEG
> alpha–theta dynamics during mind wandering in the context of breath focus
> meditation: An experience sampling approach with novice meditation
> practitioners". They also have used fieldtrip for permutation data
> analysis, but the code they provide doesn't seem to function for the
> alpha:theta ratio data.
>
> I have recorded my data in a within-subject study with a baseline and two
> independent variables. I have my data separated between these three classes
> in a 4D array with dimensions, subject x trial/epoch x EEG channel x ratio
> data. Below is the code on how I convert this data to be fieldtrip friendly:
> function [DATA, layout]=permutation(data, Participant_index, index)
> load('electrode19.mat')
> load('chanlocs.mat')
> DATA = cell(1, max(unique(Participant_index)));
> %for each participant, load data in struct form and prepare for
> %fieldtrip
> for i=1:max(unique(Participant_index))
> clear Data
> Data=struct();
> Data.label=electrode19;
> Data.hdr.nChans=numel(electrode19);
> Data.hdr.label=electrode19';
> Data.hdr.Fs=500;
> Data.fsample=500;
> data_sel = data(Participant_index(index) == i, :,:,:);
> if size(data_sel, 1) ~= 1
> data_sel = reshape(data_sel, [], size(data_sel, 3), size(data_sel, 4));
> else
> data_sel = squeeze(data_sel);
> end
> Data.trial = data_sel;
> Data.time=1:0.1:3.5;
> Data.hdr.nSamples=1;
> Data.dimord = 'chan_time';
> try
> cfg=[];
> Data=ft_preprocessing(cfg, Data);
> Data = ft_timelockanalysis(cfg, Data);
> Data = ft_datatype_timelock(Data);
> % Data.trialinfo=labels;
> cfg.layout = 'biosemi32.lay';
> Data.layout = ft_prepare_layout(cfg);
> layout = ft_prepare_layout(cfg);
> Data.layout=Data;
> DATA{i} = Data;
> catch
> DATA{i} = NaN;
> end
> end
> end
> and below is the code for which I use this converted data and apply it to
> the command ft_timelockstatistics:
> clear cls
> %Loads the directory names and the processed data
> dirs = dir('C:/Users/merli/OneDrive - University College
> London/Desktop/Study_!/P*/*.mat');
> load('C:\Users\merli\OneDrive - University College
> London\Desktop\Study_!\Data.mat');
> %For loop and unique command seperates each class and participant so we
> can extract an index
> for i=1:size(dirs,1)
> class = split(dirs(i).name, '_');
> Par = split(dirs(i).folder, 'P');
> cls{i} = class{1};
> P{i} = Par{2};
> end
> [uniqueStrings, ~, index] = unique(cls);
> [PuniqueStrings, ~, Pindex] = unique(P);
> %Seperate data to each class based on the previous index
> %percentage ratio is a 4D array with dimensions subject x trial/epoch x
> %channel x frequency data
> %Testing error by removing incomplete sessions
> % x = percentage_ratio(3:58,:,:,:);
> % size(x)
> % index = index(3:58);
> % Pindex = Pindex(3:58);
> BaselineData = percentage_ratio(index == 1, :, :, :);
> FrameData = percentage_ratio(index == 2, :, :, :);
> ShamData = percentage_ratio(index == 3, :, :, :);
> % Pindex(index == 1) == 1
> %Converts 4D array to timelock data for fieldtrip analysis
> [BDataP, BlayoutP] = permutation(BaselineData, Pindex, index == 1);
> [FDataP, FlayoutP] = permutation(FrameData, Pindex, index == 2);
> [SDataP, SlayoutP] = permutation(ShamData, Pindex, index == 3);
> % figure,plot(1:0.1:3.5, BDataP{14}.avg)
> % figure,plot(1:0.1:3.5, FDataP{14}.avg)
> % figure,plot(1:0.1:3.5, SDataP{14}.avg)
> %configuration definition for permutation statistics
> cfg=[];
> cfg.method = 'triangulation';
> cfg.layout = BlayoutP;
> neighbours = ft_prepare_neighbours(cfg);
> cfg.method = 'montecarlo'; % use the Monte Carlo Method to calculate the
> significance probability
> cfg.statistic = 'ft_statfun_depsamplesT';% use dependent samples t
> statistic
> cfg.correctm = 'cluster';
> cfg.clusteralpha = 0.05; % alpha level of the sample-specific test
> statistic that will be used for thresholding
> cfg.clusterstatistic = 'maxsum'; % test statistic that will be evaluated
> under the permutation distribution.
> cfg.neighbours = neighbours; % definition of neighbours
> cfg.tail = 0; %two-sided test
> cfg.clustertail = 0;
> cfg.alpha = 0.025; % alpha level of the permutation test
> cfg.numrandomization = 1000; % number of draws from the permutation
> distribution
> %define design matrix
> Bsubj=size(BDataP,2);
> Ssubj=size(SDataP,2);
> Fsubj=size(FDataP,2);
> %Currently instead of using all participants, am just doing it for some
> %as it'll cause errors when using everyones data
> cfg.design(1,:) = ([2:8 2:8]); %Each number represents the time it
> happened
> cfg.design(2, :) = ([ones(1,7)*1 ones(1,7)*2]); %Each number is the
> different independant variable
> 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_timelockstatistics(cfg, SDataP{2:8}, FDataP{2:8});
>
>
> This code is currently functioning, but not as intended. The variable
> stat.stat is giving a channel x len(ratio data) array
> which I assume is supposed to be the probability between both conditions,
> but all values are 1, except two random
> values that equal 0.23. I don't understand why I'm reciving these results
> as I assume it's impossible to recieve a
> probability of 1 unless the data is exactly the same for each observation.
> Additionally, when I run the code it states:
> Warning: the number of randomizations (1000) is larger than the maximum
> number of unique
> permutations (128), better use cfg.numrandomization='all'
> Even though I have 601 different trials/epochs for each class per subject.
> Can anyone identify if I have made any mistakes or point me in the
> direction required to analyse my data correctly. Any help would be
> appreciated.
> Best
> Merlin
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20230908/16e137ab/attachment.htm>