WildLab/code/collectUnfold_calculateTFCE.m

%% collectUnfold_calculateTFCE.m
% This script collects all unfold data and calculates the TFCE

cd '/net/store/nbp/users/agert/itw'
init_itw

projectFolder=['/net/store/nbp/projects/IntoTheWild'];

codingscheme = 'mashup';
filtType='causal';
phase='WLFO';

timewin=[-500 1000];


%% Collecting unfold

d2nd = [];
for subj = [8 9 29:36 38:41 44 45 48 51:53]
    
    fPath = fullfile(projectFolder,['/Daten/EEG/' phase '/unfold_' phase '/' filtType '/']);
    fPath = fullfile(fPath,['ufresult_subj' num2str(subj) ,'_allElec_' codingscheme '_' num2str(timewin(1)) '-' num2str(timewin(2)) '_regularized0.mat']);
    
    
    try
      
        fprintf('loading subject %i \n',subj)
        
        d = load(fPath);    
        ufresult=d.ufresult;

        
        ufresult = uf_predictContinuous(ufresult,'predictAt',{{'sac_amplitude',[0.5 1 2.5 5 7.5 10]},{'fix_avgpos_y',linspace(500,1500,5)},{'fix_avgpos_x',linspace(900,3000,5)},{'samebox',1}});
        % I don't have to add marginals, as I am interested in the pure
        % betas = differences between conditions
        
        if strcmp(phase,'Lab') && subj==51 % sub 51 does not have predictore sacAmp -> adjust structure accordingly
            tmpBeta(:,:,1:11)=ufresult.beta(:,:,1:11);
            tmpBeta(:,:,12:17)=nan(128,length(ufresult.times),6);
            tmpBeta(:,:,18:24)=ufresult.beta(:,:,12:18);
            tmpBeta_nodc(:,:,1:11)=ufresult.beta_nodc(:,:,1:11);
            tmpBeta_nodc(:,:,12:17)=nan(128,length(ufresult.times),6);
            tmpBeta_nodc(:,:,18:24)=ufresult.beta_nodc(:,:,12:18);
            
            ufresult.beta=tmpBeta;
            ufresult.beta_nodc=tmpBeta_nodc;
            
            if ~isempty(d2nd)
                ufresult.param=d2nd.param;
            end
        end
        
 % first subject: initialize d2nd (2nd-level analysis)
        if isempty(d2nd)
            d2nd = ufresult;
            d2nd.unfold(1)    = d2nd.unfold;
            d2nd.subject      = subj;
            d2nd.chanlocs     = ufresult.chanlocs;
            d2nd.codingscheme = codingscheme;
        else
            d2nd.unfold(end+1)          = ufresult.unfold;
            d2nd.beta(:,:,:,end+1)      = ufresult.beta;
            d2nd.beta_nodc(:,:,:,end+1) = ufresult.beta_nodc;
            d2nd.subject(end+1)         = subj;
        end
    catch e
        display(['problem with subject ',num2str(subj), ':   ' ,e.message])
    end
end


%% TFCE

elecfind = @(x)find(strcmp(x,{d2nd.chanlocs.labels}));

cfg = [];
cfg.chan = 1:length(d2nd.chanlocs); % the other channels are eye-movement channels from the Eye-EEG toolbox
cfg.blcorrect = 0;

cfg.save = 1;
ept_tfce_nb = ept_ChN2(d2nd.chanlocs,0);
cfg.nperm=10000;

cfgPlot = [];
cfgPlot.highlighted_channel= [elecfind('PO8') elecfind('P8') elecfind('PO7') elecfind('P7')];
cfgPlot.colormap = {{'div','RdYlBu'},'seq'};
cfgPlot.topoalpha = 0.05; % where to put the significance dots?
cfgPlot.individualcolorscale = 'row'; % different rows have very different interpretation
cfgPlot.time = [-0.2 0.5]; % we will zoom in
cfgPlot.n_topos = 14; %every 50 ms
cfgPlot.quality = 100; % put higher for better quality

paramNames={d2nd.param.name};

cfg.effecttype = 'prev'; %'prev' or 'curr' - only relevant for HF
paramlist = {'HF'};%HF,'interaction','samebox'};
    
data=0;
%%
for param = paramlist
    
    
    switch param{1}
        case 'HF'
            [~,paramPos]=find(ismember(paramNames,[cfg.effecttype '_HF']));
            assert(~isempty(paramPos))
            data = squeeze(d2nd.beta(cfg.chan,:,paramPos,:));
            data=data*2; %because the data is effect coded
            
        case 'interaction'
            [~,paramPos]=find(ismember(paramNames,'humanface:prevhumanface'));
            assert(~isempty(paramPos))
            data = squeeze(d2nd.beta(cfg.chan,:,paramPos,:));
            data=data*2;
        case 'samebox'
            [~,paramPos]=find(ismember(paramNames,'samebox_1'));
            assert(~isempty(paramPos))
            data = squeeze(d2nd.beta(cfg.chan,:,paramPos,:));
        otherwise
            error('wrong condition')
  
    end
    
    
    if cfg.blcorrect
        data = bsxfun(@minus,data,mean(data(:,d2nd.times<0,:),2));
    end
    
    [res,info] = be_ept_tfce_diff(struct('nperm',cfg.nperm,'neighbours',ept_tfce_nb(cfg.chan,:)),permute(data,[3 1 2]));
    cfgPlot.pvalues = res.P_Values;
    
    TFCEresults.res=res;
    TFCEresults.info=info;
    TFCEresults.channels=d2nd.chanlocs;
    TFCEresults.times=d2nd.times;
   
    
    %% actual plotting

    tmpPVal = log10(res.P_Values); % change p-values to log-scale
    
    hA =  plot_topobutter(cat(3,mean(data,3),tmpPVal),d2nd.times,d2nd.chanlocs,cfgPlot);
    hA.topo.colorbar{2}.XTick = log10([0.001 0.005 0.05]);
    hA.topo.colorbar{2}.TickLabels = [0.001 0.005 0.05];
    

   
        if strcmp(param,'interaction')
            figurename = [codingscheme ' ' param{1} ', iter:' num2str(cfg.nperm) ', time:' num2str(timewin(1)) ' to ' num2str(timewin(2))];
        elseif strcmp(param,'samebox') || strcmp(param,'boxn2')
            figurename = [codingscheme ' ' param{1} ', iter:' num2str(cfg.nperm) ', time:' num2str(timewin(1)) ' to ' num2str(timewin(2))];
        else
            
            figurename = [codingscheme ' main effect: ' cfg.effecttype ', iter:' num2str(cfg.nperm) ', time:' num2str(timewin(1)) ' to ' num2str(timewin(2))];
        end

    title(figurename)
    drawnow
    set(gcf,'Position',[5  515  2550  571])
    drawnow
   
end