code

Contains the code to do the deconvolution and plotting

code/ITW_unfold_calcERPs.m

@@ -0,0 +1,65 @@
+function [n2n,n2h,h2n,h2h,n2n_nodc,n2h_nodc,h2n_nodc,h2h_nodc,sb]=ITW_unfold_calcERPs(d2nd,codingscheme,chan,phase)
+
+subjects=d2nd.subject;
+if  strcmp(codingscheme,'mashup')
+    if strcmp(phase,'WLFO')
+        [~,parampos]= ismember({'(Intercept)','curr_HF','prev_HF','humanface:prevhumanface','samebox_1'},{d2nd.param.name});
+    elseif strcmp(phase,'Lab')
+        [~,parampos]= ismember({'2_(Intercept)','curr_HF','prev_HF','humanface:prevhumanface'},{d2nd.param.name});
+        
+    end
+    
+end
+
+n2n=[]; n2h=[]; h2n=[]; h2h=[];
+n2n_nodc=[]; n2h_nodc=[]; h2n_nodc=[]; h2h_nodc=[];
+sb=[];
+
+for sub_num=1:length(subjects)
+    
+    
+    fprintf('calculating subject %i \n',subjects(sub_num))
+    
+    unfold2=d2nd;
+    unfold2.unfold = unfold2.unfold(sub_num);
+    unfold2.beta=squeeze(d2nd.beta(:,:,:,sub_num));
+    
+    unfold2.beta_nodc=d2nd.beta_nodc(:,:,:,sub_num);
+    
+    if strcmp(phase,'WLFO')
+        uf2am = uf_predictContinuous(unfold2,'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',1}});
+        uf2am = uf_addmarginal(uf2am);
+    elseif strcmp(phase,'Lab')
+        
+        uf2am=unfold2;
+    end
+    
+    n2n_data = uf2am.beta(chan,:,parampos(1))-unfold2.beta(chan,:,parampos(3))-unfold2.beta(chan,:,parampos(2))+unfold2.beta(chan,:,parampos(4));
+    n2h_data = uf2am.beta(chan,:,parampos(1))-unfold2.beta(chan,:,parampos(3))+unfold2.beta(chan,:,parampos(2))-unfold2.beta(chan,:,parampos(4));
+    h2n_data = uf2am.beta(chan,:,parampos(1))+unfold2.beta(chan,:,parampos(3))-unfold2.beta(chan,:,parampos(2))-unfold2.beta(chan,:,parampos(4));
+    h2h_data = uf2am.beta(chan,:,parampos(1))+unfold2.beta(chan,:,parampos(3))+unfold2.beta(chan,:,parampos(2))+unfold2.beta(chan,:,parampos(4));
+    
+    n2n(sub_num,:,:)=n2n_data;
+    n2h(sub_num,:,:)=n2h_data ;
+    h2n(sub_num,:,:)=h2n_data;
+    h2h(sub_num,:,:)=h2h_data;
+    
+    n2n_nodc_data = uf2am.beta_nodc(chan,:,parampos(1))-unfold2.beta_nodc(chan,:,parampos(3))-unfold2.beta_nodc(chan,:,parampos(2))+unfold2.beta_nodc(chan,:,parampos(4));
+    n2h_nodc_data = uf2am.beta_nodc(chan,:,parampos(1))-unfold2.beta_nodc(chan,:,parampos(3))+unfold2.beta_nodc(chan,:,parampos(2))-unfold2.beta_nodc(chan,:,parampos(4));
+    h2n_nodc_data = uf2am.beta_nodc(chan,:,parampos(1))+unfold2.beta_nodc(chan,:,parampos(3))-unfold2.beta_nodc(chan,:,parampos(2))-unfold2.beta_nodc(chan,:,parampos(4));
+    h2h_nodc_data = uf2am.beta_nodc(chan,:,parampos(1))+unfold2.beta_nodc(chan,:,parampos(3))+unfold2.beta_nodc(chan,:,parampos(2))+unfold2.beta_nodc(chan,:,parampos(4));
+    
+    n2n_nodc(sub_num,:,:)=n2n_nodc_data;
+    n2h_nodc(sub_num,:,:)=n2h_nodc_data ;
+    h2n_nodc(sub_num,:,:)=h2n_nodc_data;
+    h2h_nodc(sub_num,:,:)=h2h_nodc_data;
+    
+    %samebox data
+    if strcmp(phase, 'WLFO')
+        sb_data=uf2am.beta(chan,:,parampos(1))+unfold2.beta(chan,:,parampos(3))+unfold2.beta(chan,:,parampos(2))+unfold2.beta(chan,:,parampos(4))+unfold2.beta(chan,:,parampos(5));
+        sb(sub_num,:,:)=sb_data;
+    end
+ 
+end
+
+end

code/collectUnfold_calculateTFCE.m

@@ -0,0 +1,162 @@
+%% 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

code/collectUnfold_plotERPs.m

@@ -0,0 +1,378 @@
+%%
+% This script loads all the deconvolved data and plots them as ERPs
+
+
+cd '/net/store/nbp/users/agert/itw'
+init_itw
+
+projectFolder=['/net/store/nbp/projects/IntoTheWild'];
+
+codingscheme = 'mashup';
+phase='WLFO';
+filtType='causal';
+
+timewin=[-0.5 1.0];
+
+
+% LOADING
+% needed for chanlocs
+if strcmp(phase,'Lab')
+    Data=load(['/net/store/nbp/projects/IntoTheWild/Daten/EEG/' phase '/collectedERPs_-5001000_causal_noCar.mat']);
+elseif strcmp(phase,'WLFO')
+    Data=load(['/net/store/nbp/projects/IntoTheWild/Daten/EEG/' phase '/collectedERPs_-5001000_causal.mat']);
+    
+end
+
+
+fns=fieldnames(Data);
+collData = Data.(fns{1});
+
+paramNames=[];
+%
+
+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 '_-500-1000_regularized0.mat']);
+    try
+        % load subject results
+        if ~exist(fPath,'file')
+            continue
+        end
+        fprintf('\n \n loading subject %i \n',subj)
+        
+        d = load(fPath);
+        
+        ufresult = d.ufresult;
+        
+        tmpBeta=[];
+        tmpBeta_nodc=[];
+        if strcmp(phase,'Lab') && subj==51 % sub 51 does not have predictore sacAmp -> adjust structure accordingly
+            tmpBeta(:,:,1:9)=ufresult.beta(:,:,1:9);
+            tmpBeta(:,:,10:13)=nan(length(d2nd.chanlocs),length(ufresult.times),4);
+            tmpBeta(:,:,14:20)=ufresult.beta(:,:,10:16);
+            tmpBeta_nodc(:,:,1:9)=ufresult.beta_nodc(:,:,1:9);
+            tmpBeta_nodc(:,:,10:13)=nan(length(d2nd.chanlocs),length(ufresult.times),4);
+            tmpBeta_nodc(:,:,14:20)=ufresult.beta_nodc(:,:,10:16);
+            
+            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;
+        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
+
+
+%% Calculate ERPs
+
+elecfind = @(x)find(strcmp(x,{d2nd.chanlocs.labels}));
+chan = 1:length(d2nd.chanlocs);%[elecfind('PO8')];
+
+if strcmp(codingscheme,'effects')
+    [x2n,x2h,n2x,h2x,x2n_nodc,x2h_nodc,n2x_nodc,h2x_nodc,SB]=ITW_unfold_calcERPs(d2nd,codingscheme,chan,phase);
+else
+    [n2n,n2h,h2n,h2h,n2n_nodc,n2h_nodc,h2n_nodc,h2h_nodc,SB]=ITW_unfold_calcERPs(d2nd,codingscheme,chan,phase);
+end
+
+
+%% Save in ERP structure
+
+ufERP.times=d2nd.times;
+ufERP.chanlocs=d2nd.chanlocs;
+ufERP.subject=d2nd.subject;
+ufERP.data=[];
+ufERP.data(:,:,:,1)=n2n;
+ufERP.data(:,:,:,2)=n2h;
+ufERP.data(:,:,:,3)=h2n;
+ufERP.data(:,:,:,4)=h2h;
+ufERP.cond={'n2n','n2h','h2n','h2h'};
+
+
+save(['/net/store/nbp/projects/IntoTheWild/Daten/EEG/' phase '/unfold_' phase '/' filtType '/unfoldERPs_-5001000.mat'],'ufERP')
+
+
+
+%% PLOT ERPs
+%  load(['/net/store/nbp/projects/IntoTheWild/Daten/EEG/WLFO/unfold_WLFO/' filtType '/unfoldERPs_-5001000.mat'],'ufERP')
+
+elecfind = @(x)find(strcmp(x,{d2nd.chanlocs.labels}));
+chan =[elecfind('PO8'),elecfind('P8'),elecfind('PO7'),elecfind('P7')];
+
+% plot 2x2 results of deconvolution
+    figure('units','normalized','outerposition',[0 0 1 1])
+    hold on,  plot(d2nd.times,squeeze(mean(mean(n2n(:,chan,:)),2)),'b--','LineWidth',2);
+    hold on, plot(d2nd.times,squeeze(mean(mean(n2h(:,chan,:)),2)),'r--','LineWidth',2);
+    hold on, plot(d2nd.times,squeeze(mean(mean(h2n(:,chan,:)),2)),'b','LineWidth',2);
+    hold on, plot(d2nd.times,squeeze(mean(mean(h2h(:,chan,:)),2)),'r','LineWidth',2);
+    
+    legend({'BG -> BG','BG -> HF','HF -> BG','HF -> HF'},'FontSize',20)
+    xlim([-0.5 1])
+    xt = get(gca, 'XTick');
+    set(gca, 'FontSize', 18)
+    ylim([-6 6]);
+    vline(0,'k')
+    hline(0,'k');
+    box off
+    title([phase ' Deconvolution'],'FontSize',20)
+    xlabel('Time [sec]', 'FontSize',20)
+    ylabel('Amplitude [microV]', 'FontSize',20);
+    
+%     export_fig(['/home/student/a/agert/Desktop/Into the Wild/LabvsWLFO/autoplot/2x2_ERP_Deconv_Unfold_' filtType '.eps'])
+    
+    % plot 2x2 design of NO deconvolution
+    figure('units','normalized','outerposition',[0 0 1 1])
+    hold on, plot(d2nd.times,squeeze(mean(mean(n2n_nodc(:,chan,:)),2)),'b--','LineWidth',2);
+    hold on, plot(d2nd.times,squeeze(mean(mean(n2h_nodc(:,chan,:)),2)),'r--','LineWidth',2);
+    hold on, plot(d2nd.times,squeeze(mean(mean(h2n_nodc(:,chan,:)),2)),'b','LineWidth',2);
+    hold on, plot(d2nd.times,squeeze(mean(mean(h2h_nodc(:,chan,:)),2)),'r','LineWidth',2);
+    
+    legend({'BG -> BG','BG -> HF','HF -> BG','HF -> HF'},'FontSize',20)
+    xlim([-0.5 1])
+    xt = get(gca, 'XTick');
+    set(gca, 'FontSize', 18)
+    ylim([-6 6])
+    vline(0,'k')
+    hline(0,'k')
+    box off
+    title([phase ' No Deconvolution'],'FontSize',20)
+    xlabel('Time [sec]', 'FontSize',20)
+    ylabel('Amplitude [microV]', 'FontSize',20);
+
+    
+    % within face fixation
+    if strcmp(phase,'WLFO')
+        figure('units','normalized','outerposition',[0 0 1 1])
+        hold on,plot(d2nd.times,squeeze(mean(mean(n2h(:,chan,:)),2)),'b','LineWidth',2);
+        hold on, plot(d2nd.times,squeeze(mean(mean(h2h(:,chan,:)),2)),'r--','LineWidth',2);
+        hold on, plot(d2nd.times,squeeze(mean(mean(SB (:,chan,:)),2)),'r','LineWidth',2);
+        
+        
+        legend({'BG -> HF','HF -> HF','Within Box'},'FontSize',20)
+        xlim([-0.5 1])
+        xt = get(gca, 'XTick');
+        set(gca, 'FontSize', 18)
+        ylim([-6 6])
+        vline(0,'k')
+        % vline(0.17,'--k')
+        hline(0,'k')
+        box off
+        title([phase ' Samebox Effect'],'FontSize',20)
+        xlabel('Time [sec]', 'FontSize',20)
+        ylabel('Amplitude [microV]', 'FontSize',20);
+        %     export_fig(['/home/student/a/agert/Desktop/Into the Wild/LabvsWLFO/autoplot/SB_ERP_' filtType '.png'])
+        %     export_fig(['/home/student/a/agert/Desktop/Into the Wild/LabvsWLFO/autoplot/SB_ERP.eps'])
+    end
+    %% 1x2 ERP plots
+    n2x=(n2h+n2n)./2;
+    h2x=(h2h+h2n)./2;
+    x2n=(h2n+n2n)./2;
+    x2h=(h2h+n2h)./2;
+    
+    % plot 1x2 for type of previous fixation
+    figure('units','normalized','outerposition',[0 0 1 1])
+    hold on,  plot(d2nd.times,squeeze(mean(mean(n2x(:,chan,:)),2)),'k--','LineWidth',2);
+    hold on, plot(d2nd.times,squeeze(mean(mean(h2x(:,chan,:)),2)),'k','LineWidth',2);
+    
+    
+    legend({'BG -> X','HF -> X'},'FontSize',20)
+    xlim([-0.5 1])
+    xt = get(gca, 'XTick');
+    set(gca, 'FontSize', 18)
+    ylim([-6 6])
+    vline(0,'k')
+    % vline(0.17,'--k')
+    hline(0,'k')
+    box off
+    title([phase ' Main Effect Prev'],'FontSize',20)
+    xlabel('Time [sec]', 'FontSize',20)
+    ylabel('Amplitude [microV]', 'FontSize',20);
+    %export_fig(['/home/student/a/agert/Desktop/Into the Wild/LabvsWLFO/autoplot/1x2_ERP_MainPrev_Unfold_' filtType '.eps'])
+    
+    % plot 1x2 for type of current fixation
+    x2nCI=zeros(length(squeeze(mean(x2n(:,chan,:),2))),2); 
+    x2hCI=zeros(length(squeeze(mean(x2h(:,chan,:),2))),2); 
+    for timep=1:length(squeeze(mean(x2n(:,chan,:),2)))
+        x2nCI(timep,:)= bootci(10000,@mean,mean(x2n(:,chan,timep),2));
+        x2hCI(timep,:)= bootci(10000,@mean,mean(x2h(:,chan,timep),2));
+    end
+    
+    figure('units','normalized','outerposition',[0 0 1 1])
+%     hold on, plot(d2nd.times,squeeze(mean(mean(x2h(:,chan,:)),2)),'r','LineWidth',2);
+%     
+    x1 = [squeeze(d2nd.times), fliplr(squeeze(d2nd.times))];
+    inBetween = [x2nCI(:,1)', fliplr(x2nCI(:,2)')];
+    fill(x1, inBetween, [0,0.5,0.5]);
+    hold on;
+    plot(d2nd.times,squeeze(mean(mean(x2n(:,chan,:)),2)),'b','LineWidth',2);
+%     hold on, plot(d2nd.times,x2nCI,'b--','LineWidth',2);
+%     hold on, plot(d2nd.times,x2hCI,'r--','LineWidth',2);
+
+    inBetween2 = [x2hCI(:,1)', fliplr(x2hCI(:,2)')];
+    fill(x1, inBetween2, [1,0.5,0]);
+    hold on;
+    plot(d2nd.times,squeeze(mean(mean(x2h(:,chan,:)),2)),'r','LineWidth',2);
+
+    legend({'X -> BG','X -> HF'},'FontSize',20)
+    xlim([-0.5 1])
+    xt = get(gca, 'XTick');
+    set(gca, 'FontSize', 18)
+    ylim([-7 7])
+    vline(0,'k')
+    vline(0.12,'k--') %N170 start
+    vline(0.2,'k--') %N170 stop
+    vline(-0.3,'k--')
+    % vline(0.17,'--k')
+    hline(0,'k')
+    
+    box off
+    title([phase ' Main Effect Curr'],'FontSize',20)
+    xlabel('Time [sec]', 'FontSize',20)
+    ylabel('Amplitude [microV]', 'FontSize',20);
+    export_fig(['/home/student/a/agert/Desktop/Into the Wild/LabvsWLFO/autoplot/1x2_ERP_MainCurr_Unfold_' filtType '_' phase '_withCI.eps'])
+   
+
+
+%% Plot effect topos (timerange specified for N170)
+
+effectData=squeeze(mean(x2h(:,:,:)))-squeeze(mean(x2n(:,:,:)));
+tp1=find(d2nd.times>0.120, 1, 'first');
+tp2=find(d2nd.times<0.200, 1, 'last');
+topoData=squeeze(mean(effectData(:,tp1:tp2),2));
+figure, topoplot(topoData',d2nd.chanlocs,'conv','on');
+title(['N170 topo 120:200 ' phase])
+
+%% viewing behavior
+d2nd2 = subjectwise_getparam(d2nd,{{'sac_amplitude',[0.5 1 2.5 5 7.5 10]}},1);
+
+uf_plotParam(d2nd2,'channel','Oz','plotSeparate','all','plotParam',{'fix_avgpos_y','fix_avgpos_x'})
+uf_plotParam(d2nd2,'channel','Oz', 'plotSeparate','all','plotParam','sac_amplitude','deconv',1)
+set(findobj(gca,'-property','LineWidth'),'LineWidth',2)
+
+%% Component amplitudes
+%P100: 80-130
+%N170: 120-200
+
+
+p100Start=find(d2nd.times>=0.080,1,'first');
+p100End=find(d2nd.times<=0.130,1,'last');
+n170Start=find(d2nd.times>=0.120,1,'first');
+n170End=find(d2nd.times<=0.200,1,'last');
+
+meanN2N.p100=squeeze(mean(ufERP.data(:,chan,p100Start:p100End,1),2));
+meanN2H.p100=squeeze(mean(ufERP.data(:,chan,p100Start:p100End,2),2));
+meanH2N.p100=squeeze(mean(ufERP.data(:,chan,p100Start:p100End,3),2));
+meanH2H.p100=squeeze(mean(ufERP.data(:,chan,p100Start:p100End,4),2));
+
+meanN2N.n170=squeeze(mean(ufERP.data(:,chan,n170Start:n170End,1),2));
+meanN2H.n170=squeeze(mean(ufERP.data(:,chan,n170Start:n170End,2),2));
+meanH2N.n170=squeeze(mean(ufERP.data(:,chan,n170Start:n170End,3),2));
+meanH2H.n170=squeeze(mean(ufERP.data(:,chan,n170Start:n170End,4),2));
+
+
+meanX2N.p100=squeeze(mean(mean(ufERP.data(:,chan,p100Start:p100End,[1,3]),2),4));
+meanX2H.p100=squeeze(mean(mean(ufERP.data(:,chan,p100Start:p100End,[2,4]),2),4));
+
+meanX2N.n170=squeeze(mean(mean(ufERP.data(:,chan,n170Start:n170End,[1,3]),2),4));
+meanX2H.n170=squeeze(mean(mean(ufERP.data(:,chan,n170Start:n170End,[2,4]),2),4));
+
+
+meanN2X.p100=squeeze(mean(mean(ufERP.data(:,chan,p100Start:p100End,[1,2]),2),4));
+meanH2X.p100=squeeze(mean(mean(ufERP.data(:,chan,p100Start:p100End,[3,4]),2),4));
+
+meanN2X.n170=squeeze(mean(mean(ufERP.data(:,chan,n170Start:n170End,[1,2]),2),4));
+meanH2X.n170=squeeze(mean(mean(ufERP.data(:,chan,n170Start:n170End,[3,4]),2),4));
+
+
+
+%%
+p100=meanH2H.p100;
+n170=meanH2H.n170;
+[p100Amp,p100Time]=max(p100,[],2);
+[n170Amp,n170Time]=min(n170,[],2);
+p2pEff=p100Amp-n170Amp;
+p100PeakTime=ufERP.times(p100Start+p100Time-1);
+n170PeakTime=ufERP.times(n170Start+n170Time-1);
+p100AmpCI = bootci(10000,@mean,p100Amp);
+n170AmpCI = bootci(10000,@mean,n170Amp);
+p2pAmpCI  = bootci(10000,@mean,p2pEff);
+p100TimeCI = bootci(10000,@mean,p100PeakTime);
+n170TimeCI = bootci(10000,@mean,n170PeakTime);
+
+t=0;
+
+if t==100
+    fprintf(['P100 for ' phase ' H2H on average ' num2str(mean(p100Amp)) 'uV 95CI [' num2str(p100AmpCI(1)) ';' num2str(p100AmpCI(2))  '] at ' num2str(mean(p100PeakTime)) 'ms [' num2str(p100TimeCI(1)) ';' num2str(p100TimeCI(2))  '].\n'])
+elseif t==170
+    fprintf(['N170 for ' phase ' H2H on average ' num2str(mean(n170Amp)) 'uV 95CI [' num2str(n170AmpCI(1)) ';' num2str(n170AmpCI(2))  '] at ' num2str(mean(n170PeakTime)) 'ms [' num2str(n170TimeCI(1)) ';' num2str(n170TimeCI(2))  '].\n'])
+elseif t==0
+    fprintf(['P2P for ' phase ' H2H on average ' num2str(mean(p2pEff)) 'uV 95CI [' num2str(p2pAmpCI(1)) ';' num2str(p2pAmpCI(2))  '].\n'])
+    
+end
+
+
+%% conditionwise differences
+
+cond1=meanX2N;
+cond2=meanX2H;
+
+[cond1p100Amp,cond1p100Time]=max(cond1.p100,[],2);
+[cond1n170Amp,cond1n170Time]=min(cond1.n170,[],2);
+[cond2p100Amp,cond2p100Time]=max(cond2.p100,[],2);
+[cond2n170Amp,cond2n170Time]=min(cond2.n170,[],2);
+condDiffp100=cond2p100Amp-cond1p100Amp;
+condDiffn170=cond2n170Amp-cond1n170Amp;
+p100AmpCI = bootci(10000,@mean,condDiffp100);
+n170AmpCI = bootci(10000,@mean,condDiffn170);
+
+fprintf(['Diff P100 for ' phase ' on average ' num2str(mean(condDiffp100)) 'uV 95CI [' num2str(p100AmpCI(1)) ';' num2str(p100AmpCI(2))  '].\n'])
+
+fprintf(['Diff N170 for ' phase ' on average ' num2str(mean(condDiffn170)) 'uV 95CI [' num2str(n170AmpCI(1)) ';' num2str(n170AmpCI(2))  '].\n'])
+
+
+% [N2NAmp,N2NTime]=min(meanN2N,[],2);
+% [N2HAmp,N2HTime]=min(meanN2H,[],2);
+% [H2NAmp,H2NTime]=min(meanH2N,[],2);
+% [H2HAmp,H2HTime]=min(meanH2H,[],2);
+%
+%
+% [X2HAmp,X2HTime]=max(meanX2H,[],2);
+% [X2NAmp,X2NTime]=max(meanX2N,[],2);
+%
+% X2NPeakTime=ufERP.times(n170Start+X2NTime);
+% X2HPeakTime=ufERP.times(n170Start+X2HTime);
+
+% X2NAmpCI = bootci(10000,@mean,X2NAmp);
+% X2HAmpCI = bootci(10000,@mean,X2HAmp);
+% X2NTimeCI = bootci(10000,@mean,X2NPeakTime);
+% X2HTimeCI = bootci(10000,@mean,X2HPeakTime);
+
+% fprintf(['N170 for ' phase ' N2N on average ' num2str(mean(N2NRange)) 'uV.\n'])
+% fprintf(['N170 for ' phase ' N2H on average ' num2str(mean(N2HRange)) 'uV.\n'])
+% fprintf(['N170 for ' phase ' H2N on average ' num2str(mean(H2NRange)) 'uV.\n'])
+% fprintf(['N170 for ' phase ' H2H on average ' num2str(mean(H2HRange)) 'uV.\n'])
+
+% fprintf(['N170 for ' phase ' X2N on average ' num2str(mean(X2NAmp)) 'uV 95CI [' num2str(X2NCI(1)) ' ' num2str(X2NAmpCI(2))  '] at ' num2str(mean(X2NPeakTime)) 'ms.\n'])
+% fprintf(['N170 for ' phase ' X2H on average ' num2str(mean(X2HAmp)) 'uV 95CI [' num2str(X2HCI(1)) ' ' num2str(X2HAmpCI(2))  '] at ' num2str(mean(X2HPeakTime)) 'ms.\n'])
+
+

code/unfold/unfold_calculateSE.m

@@ -0,0 +1,157 @@
+% loads unfold results, adds ufresult.data and calculates the
+% SE over trials for each subject
+
+function se = unfold_calculateSE(subj,phase)
+cd '/net/store/nbp/users/agert/itw'
+init_itw
+
+addpath('/net/store/nbp/users/agert/unfold')
+init_unfold
+
+projectFolder=['/net/store/nbp/projects/IntoTheWild'];
+
+filtType='causal';
+% time ranges for P100 and N170
+t_p100=[0.07 0.150];
+t_n170=[0.120 0.200];
+
+
+%% LOADING
+codingscheme = 'mashup';
+
+
+%% load ufresults
+ufPath = fullfile(projectFolder,['Daten/EEG/' phase '/unfold_' phase '/' filtType '/']);
+ufPath = fullfile(ufPath,['ufresult_subj' num2str(subj) ,'_allElec_' codingscheme '_-500-1000_regularized0.mat']);
+
+fprintf('\n \n loading subject %i \n',subj)
+
+d = load(ufPath);
+ufresult = d.ufresult;
+clear d; % make some space
+
+if strcmp(phase,'WLFO')
+    ufam = 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',1}});
+    ufam = uf_addmarginal(ufam);
+elseif strcmp(phase,'Lab')
+    ufam=ufresult;
+end
+
+
+%% load EEG data
+cfg = struct();
+cfg.subject = num2str(subj);
+cfg.filtering='causal';
+cfg.mainfolder = ['/net/store/nbp/projects/IntoTheWild/Daten/EEG/' phase '/VP' cfg.subject '/preprocessed/' cfg.filtering];
+cfg.file = ['3_ITW_' phase '_subj' cfg.subject '_channelrejTriggersXensor.set'];
+tmp = pop_loadset('filename',['2_ITW_' phase '_subj' cfg.subject '_bandpass_resample_deblank.set'],'filepath',cfg.mainfolder,'loadmode','info');
+tmp.chanlocs(find(strcmp({tmp.chanlocs.labels},'CzAmp2'))).labels = 'Iz';
+tmp=pop_chanedit(tmp, 'lookup','/net/store/nbp/projects/IntoTheWild/EEG_analysis/eeglab13_5_4b/plugins/dipfit2.3/standard_BESA/standard-10-5-cap385.elp');
+cfg.urchanlocs= tmp.chanlocs;
+cfg.avgref=1;
+cfg.saccade=1;
+cfg.regularize=0;
+tmp = load(fullfile(cfg.mainfolder,['6_ITW_' phase '_subj' cfg.subject '_ICAcleancont.mat']));
+cfg.badcomponents = tmp.comps_to_rej;
+cfg.timewin=[-0.5 1];
+
+%load EEG
+EEG = pop_loadset('filename',cfg.file,'filepath',cfg.mainfolder);
+EEG =  pop_select(EEG,'nochannel',{'AUX1','AUX2'});
+
+% Load ICA and clean
+mod = loadmodout12(fullfile(cfg.mainfolder,'amica'));
+EEG.icaweights = mod.W;
+EEG.icasphere = mod.S;
+EEG.icawinv = [];EEG.icaact = [];EEG.icachansind = [];
+EEG = eeg_checkset(EEG);
+EEG = pop_subcomp(EEG,cfg.badcomponents);
+
+% interpolate missing channels
+% ET channels (Gaze etc.) not interpolated
+alldel = {'CzAmp2' 'BIP1' 'BIP2' 'BIP3' 'BIP4' 'BIP5' 'BIP6' 'BIP7' 'BIP8' 'AUX1' 'AUX2' 'AUX3' 'AUX4' 'AUX5' 'AUX6' 'AUX7' 'AUX8','TIME', 'L_GAZE_X', 'L_GAZE_Y', 'L_AREA','R_GAZE_X', 'R_GAZE_Y', 'R_AREA', 'INPUT', 'L-GAZE-X', 'L-GAZE-Y', 'L-AREA','R-GAZE-X', 'R-GAZE-Y', 'R-AREA', 'INPUT'};
+%remove VEOG from EEG, as there is no corresponding location
+EEG = pop_select(EEG, 'nochannel', {'VEOG'});
+
+if cfg.avgref
+    %calculate avg ref
+    EEG = pop_reref( EEG, []); %Participants’ averages were then re-referenced to a common average reference. (Rossion & Caharel, 2011)
+end
+
+%interpolate missing channels
+idxs = ~ismember({cfg.urchanlocs.labels},alldel);
+EEG= pop_interp(EEG,cfg.urchanlocs(idxs),'spherical');
+
+
+%% add the EEG data to the unfold data
+ufresult.data=EEG.data;
+
+clear EEG mod; % make some space
+
+
+%% find P100 and N170 peaks
+
+elecfind = @(x)find(strcmp(x,{ufresult.chanlocs.labels}));
+chan = [elecfind('PO8'),elecfind('P8'),elecfind('P7'),elecfind('PO7')];
+
+P100Start=find(ufresult.times>t_p100(1),1,'first');
+P100Stop=find(ufresult.times<t_p100(2),1,'last');
+N170Start=find(ufresult.times>t_n170(1),1,'first');
+N170Stop=find(ufresult.times<t_n170(2),1,'last');
+
+paramNames=ufresult.unfold.colnames;
+if strcmp(phase,'WLFO')
+    [~,intPos]   =find(ismember(paramNames,'(Intercept)'));
+elseif strcmp(phase,'Lab')
+    [~,intPos]   =find(ismember(paramNames,'2_(Intercept)'));
+end
+[~,currHFPos]=find(ismember(paramNames,'curr_HF'));
+
+ufFace=mean(ufam.beta(chan,:,intPos)+ufresult.beta(chan,:,currHFPos));
+ufObj =mean(ufam.beta(chan,:,intPos)-ufresult.beta(chan,:,currHFPos));
+clear ufam;% make some space
+
+[~,P100Face]=max(ufFace(P100Start:P100Stop));
+P100Face=P100Start+P100Face-1;
+[~,N170Face]=min(ufFace(N170Start:N170Stop));
+N170Face=N170Start-1+N170Face; %-1 as correction for shift by one position
+
+[~,P100Obj]=max(ufObj(P100Start:P100Stop));
+P100Obj=P100Start+P100Obj-1;
+[~,N170Obj]=min(ufObj(N170Start:N170Stop));
+N170Obj=N170Start-1+N170Obj; %-1 as correction for shift by one position
+
+
+%% calculate SE
+contrast = zeros(7,size(ufresult.unfold.beta_dc,2),size(ufresult.unfold.beta_dc,3));
+
+%y = b0==N1 + bHF==N1 - (b0==P1 + bHF==P1)
+
+contrast(1, N170Face,intPos)    = 1;
+contrast(1, N170Face,currHFPos) = 1;
+
+contrast(2, P100Face,intPos)   = (1);
+contrast(2, P100Face,currHFPos)= (1);
+
+%y = b0==N1 - bHF==N1 - (y = b0==P1 - bHF==P1)
+contrast(3, N170Obj,intPos)    = 1;
+contrast(3, N170Obj,currHFPos) = -1;
+
+contrast(4, P100Obj,intPos)   = (1);
+contrast(4, P100Obj,currHFPos)= (-1);
+
+contrast=uf_contrast_addmarginal(ufresult,contrast);
+%y =  bHF==N1+ bHF==N1  - bHF==P1  - bHF==P1   - b0==P1 + b0==P1 + b0==N1 - b0==N1
+contrast(5, :,:) =   contrast(1, :,:)-contrast(2, :,:); %N1Face - P1Face
+contrast(6, :,:) =   contrast(3, :,:)-contrast(4, :,:); %N1Obje - P1Obje
+contrast(7, :,:) =   contrast(5, :,:)-contrast(6, :,:); %P2PFace - P2Pobj
+
+tic
+[se] = uf_se(ufresult,'channels', chan,'contrast',contrast);
+
+toc
+
+save(['/net/store/nbp/projects/IntoTheWild/Daten/EEG/' phase '/unfold_' phase '/' filtType 'unfoldSE_subj' num2str(subj) '.mat'],'se')
+end
+
+%end

code/unfold/unfold_freeviewing.m

@@ -0,0 +1,358 @@
+function [ufresult]=unfold_freeviewing(sub_num)
+%%
+cfg = struct();
+cfg.subject = num2str(sub_num);
+cfg.filtering='causal';
+cfg.mainfolder = ['/net/store/nbp/projects/IntoTheWild/Daten/EEG/WLFO/VP' cfg.subject '/preprocessed/' cfg.filtering];
+tmp = load(fullfile(cfg.mainfolder,['4_ITW_WLFO_subj' cfg.subject '_cleaningTimes.mat']));
+cfg.cleantimes = tmp.rej(:,1:2);
+
+cfg.file = ['3_ITW_WLFO_subj' cfg.subject '_channelrejTriggersXensor.set'];
+tmp = pop_loadset('filename',['2_ITW_WLFO_subj' cfg.subject '_bandpass_resample_deblank.set'],'filepath',cfg.mainfolder,'loadmode','info');
+tmp=pop_chanedit(tmp, 'lookup','/net/store/nbp/projects/IntoTheWild/EEG_analysis/eeglab13_5_4b/plugins/dipfit2.3/standard_BESA/standard-10-5-cap385.elp');
+tmp.chanlocs(find(strcmp({tmp.chanlocs.labels},'CzAmp2'))).labels = 'Iz';
+cfg.urchanlocs= tmp.chanlocs;
+cfg.avgref=1;
+cfg.saccade=1;
+cfg.regularize=0;
+
+tmp = load(fullfile(cfg.mainfolder,['6_ITW_WLFO_subj' cfg.subject '_ICAcleancont.mat']));
+cfg.badcomponents = tmp.comps_to_rej;
+
+cfg.timewin=[-0.5 1];
+
+
+%% load EEG data
+EEG = pop_loadset('filename',cfg.file,'filepath',cfg.mainfolder);
+EEG =  pop_select(EEG,'nochannel',{'AUX1','AUX2'});
+
+
+%% Load ICA and clean
+addpath('/net/store/nbp/projects/EEG/blind_spot/amica')
+mod = loadmodout12(fullfile(cfg.mainfolder,'amica'));
+EEG.icaweights = mod.W;
+EEG.icasphere = mod.S;
+EEG.icawinv = [];EEG.icaact = [];EEG.icachansind = [];
+EEG = eeg_checkset(EEG);
+EEG = pop_subcomp(EEG,cfg.badcomponents);
+
+
+%% interpolate missing channels
+% ET channels (Gaze etc.) not interpolated
+alldel = {'CzAmp2' 'BIP1' 'BIP2' 'BIP3' 'BIP4' 'BIP5' 'BIP6' 'BIP7' 'BIP8' 'AUX1' 'AUX2' 'AUX3' 'AUX4' 'AUX5' 'AUX6' 'AUX7' 'AUX8','TIME', 'L_GAZE_X', 'L_GAZE_Y', 'L_AREA','R_GAZE_X', 'R_GAZE_Y', 'R_AREA', 'INPUT', 'L-GAZE-X', 'L-GAZE-Y', 'L-AREA','R-GAZE-X', 'R-GAZE-Y', 'R-AREA', 'INPUT'};
+
+%remove VEOG from EEG, as there is no corresponding location
+EEG = pop_select(EEG, 'nochannel', {'VEOG'});
+
+if cfg.avgref
+    %calculate avg ref
+    EEG = pop_reref( EEG, []); %Participants’ averages were then re-referenced to a common average reference. (Rossion & Caharel, 2011)
+end
+
+idxs = ~ismember({cfg.urchanlocs.labels},alldel);
+
+%interpolate missing channels
+EEG= pop_interp(EEG,cfg.urchanlocs(idxs),'spherical');
+EEG = eeg_checkset(EEG);
+
+%% adjust cleaning times to be between 1 and length of data
+[a,b]=find(cfg.cleantimes(:,1)<1);
+if ~isempty(a)
+    cfg.cleantimes(a,1)=1;
+end
+
+[a,b]=find(cfg.cleantimes(:,1)>length(EEG.data));
+if ~isempty(a)
+    cfg.cleantimes(a,:)=[];
+end
+
+
+[a,b]=find(cfg.cleantimes(:,2)>length(EEG.data));
+if ~isempty(a)
+    cfg.cleantimes(a,2)=length(EEG.data);
+end
+
+%% for each fixation we need to save ALL the necessary information in the EEG structure
+% for e = 1:length(EEG.event)
+for e = length(EEG.event):-1:1
+    
+    
+    % if the trigger was the beginning of a trial
+    if str2num(EEG.event(e).type) <= 171
+        EEG.event(e).urtype = EEG.event(e).type;
+        EEG.event(e).type = 'stimulus';
+        continue
+    end
+    
+    % 1-171 stimulus
+    % 180 stim offset
+    % 200 exp start
+    % 255 exp stop
+    
+    % if the trigger is a fixation
+    if ismember(EEG.event(e).type, ...
+            {'HFtoHF'    'HFtoHH'    'HFtoNH'    'HFtoOL'    'HFtoNone'  'HFtoOS' ...
+            'HHtoHF'    'HHtoHH'    'HHtoNH'    'HHtoOL'    'HHtoNone'   'HHtoOS' ...
+            'NHtoHF'    'NHtoHH'    'NHtoNH'    'NHtoOL'    'NHtoNone'   'NHtoOS' ...
+            'OLtoHF'    'OLtoHH'    'OLtoNH'    'OLtoOL'    'OLtoNone'   'OLtoOS' ...
+            'NonetoHF'  'NonetoHH'  'NonetoNH'  'NonetoOL'  'NonetoNone' 'NonetoOS' ...
+            'OStoHF'    'OStoHH'    'OStoNH'    'OStoOL'    'OStoNone'   'OStoOS' ...
+            })
+        
+        
+        EEG.event(e).urtype = EEG.event(e).type;
+        
+        % if it is a meaningful fixation
+        toIDX=strfind(EEG.event(e).type,'to');
+        
+        if toIDX>0
+            
+            EEG.event(e).curr=EEG.event(e).type((toIDX+2):end);
+            EEG.event(e).prev=EEG.event(e).type(1:toIDX-1);
+            
+        else error(['You forgot to code a fixation type:' EEG.event(e).type])
+        end
+        
+        EEG.event(e).type = 'fixation';
+        
+        %         continue      end
+        
+        
+    end
+    
+    % recode all fix that are not Face or None to "OtherHeads"
+    if strcmp(EEG.event(e).type, 'fixation')
+        isFix=1;
+    else
+        isFix=0;
+    end
+    
+    if isFix==1
+        if ismember(EEG.event(e).curr,{'OS'})==1
+            EEG.event(e).curr='outside';
+        elseif ismember(EEG.event(e).curr,{'HF','None'})==0
+            EEG.event(e).curr='OtherHeads';
+            
+        end
+        
+        if ismember(EEG.event(e).prev,{'OS'})==1
+            EEG.event(e).prev='outside';
+        elseif ismember(EEG.event(e).prev,{'HF','None'})==0
+            EEG.event(e).prev='OtherHeads';
+        end
+        
+        if EEG.event(e).prevprevhumanface==1
+            EEG.event(e).prevprev='HF';
+        elseif EEG.event(e).prevprevnone==1
+            EEG.event(e).prevprev='None';
+        elseif EEG.event(e).prevprevoutside==1
+            EEG.event(e).prevprev='outside';
+        elseif EEG.event(e).prevprevhumanhead==1
+            EEG.event(e).prevprev='OtherHeads';
+        elseif EEG.event(e).prevprevnonhuman==1
+            EEG.event(e).prevprev='OtherHeads';
+        else
+            %first in a trial has no prevprev, therefore we define it as
+            %'outside' as outside will not be analysed (easiest solution)
+            EEG.event(e).prevprev='outside';
+        end
+        
+    end
+    
+end
+
+
+%% save information about saccaded and fixations
+if sum(strcmp({EEG.event.type},'L_saccade'))>0
+    ixSac = find(strcmp({EEG.event.type},'L_saccade'));
+else
+    ixSac = find(strcmp({EEG.event.type},'R_saccade'));
+end
+
+ixFix = find(strcmp({EEG.event.type},'fixation'));
+
+%find blinks
+ixBlk=find(cellfun(@(x)strcmp(x,'R_blink'),{EEG.event.type}));
+if isempty(ixBlk)
+    ixBlk=find(cellfun(@(x)strcmp(x,'L_blink'),{EEG.event.type}));
+end
+
+latSac = [EEG.event(ixSac).latency];
+latFix = [EEG.event(ixFix).latency];
+latBlk = [EEG.event(ixBlk).latency];
+
+endSac = latSac+[EEG.event(ixSac).duration];
+endBlk = latBlk+[EEG.event(ixBlk).duration];
+
+% put the velocity for blink sac/data loss sections to NaN
+for ixS=1:length(ixSac)
+    
+    blkStarAftertSacStart=find(latBlk>=latSac(ixS),1,'first');
+    blkStopBeforeSacStop=find(endBlk<=endSac(ixS),1,'last');
+    if blkStarAftertSacStart== blkStopBeforeSacStop
+        EEG.event(ixSac(ixS)).sac_vmax=NaN;
+    end
+end
+
+for e = 1:length(latFix)
+    % for each fixation, find the closest saccade before the event and copy
+    % the saccade information
+    
+    % when  Sac --> Fix, then lat(Fix) - lat(Sac) should be positive
+    ix = find((latFix(e) - latSac) > 0,1,'last');
+    
+    for fn = {'sac_endpos_x','sac_endpos_y','sac_amplitude','sac_startpos_x','sac_startpos_y','sac_vmax'}
+        EEG.event(ixFix(e)).(fn{1})= EEG.event(ixSac(ix)).(fn{1});
+    end
+    
+end
+
+% remove all fixation with saccade velocity of 0 (= blink saccades)
+for ixF=1:length(ixFix)
+    if EEG.event(ixFix(ixF)).sac_vmax==0
+        EEG.event(ixFix(ixF)).type='blinkFix';
+    end
+end
+
+ixFix = find(strcmp({EEG.event.type},'fixation'));
+latFix =  [EEG.event(ixFix).latency];
+
+% remove all fixations close to a blink/data loss section
+if ~isempty(ixBlk)
+    ixBlk=fliplr(ixBlk);
+    for iB=ixBlk
+        blkStop=EEG.event(iB).latency+EEG.event(iB).duration;
+        delFix=find(abs(latFix-blkStop)<50);
+        
+        if ~isempty(delFix)
+            for dF=1:length(delFix)
+                EEG.event(ixFix(delFix(dF))).type='blinkFix';
+            end
+        end
+    end
+end
+
+for e = 1:length(latSac)
+    % for each fixation, find the closest saccade before the event and copy
+    % the saccade information
+    
+    % when  Sac --> Fix, then lat(Fix) - lat(Sac) should be positive
+    ix = find((latSac(e) - latFix) > 0,1,'last');
+    if ~isempty(ix) && ~isempty( EEG.event(ixFix(ix)).prev)
+        for fn = {'prev','curr','samebox'}
+            EEG.event(ixSac(e)).(fn{1})= EEG.event(ixFix(ix)).(fn{1});
+            
+        end
+        EEG.event(ixSac(e)).type='saccade';
+    end
+    
+end
+
+
+%% Remove large saccades?!
+fprintf('\n\n %u large saccades (>40) will be removed. \n\n\n', length(find([EEG.event.sac_amplitude]>40)))
+EEG.event(cellfun(@(x)(x>40),{EEG.event.sac_amplitude})) = [];
+
+
+%% Config
+EEG.unfold = [];
+cfgDesign = [];
+
+cfgDesign.splinespacing = 'quantile';
+
+cfgDesign.eventtypes = {
+    {'fixation'}   % fixation onset
+    {'stimulus'}}; % stimulus onset
+
+cfgDesign.formula = {
+    'y~1 + cat(curr) + cat(prev)+ cat(samebox) + humanface:prevhumanface + spl(fix_avgpos_x,5) + spl(fix_avgpos_y,5) +  spl(sac_amplitude,5)',
+    'y~1'};
+
+cfgDesign.codingschema = 'effects';
+cfgDesign.mashup=1;
+
+% if codingschema=effects -> reference cat is the last one
+% if codingschema=reference -> reference cat is the first one
+if strcmp(cfgDesign.codingschema, 'effects')
+    cfgDesign.categorical = {
+        'curr',{'outside','HF','OtherHeads','None'};
+        'prev',{'outside','HF','OtherHeads','None'};
+        'samebox',{1,0};
+        };
+elseif strcmp(cfgDesign.codingschema, 'reference')
+    cfgDesign.categorical = {
+        'curr',{'None','HF','OtherHeads','outside'};
+        'prev',{'None','HF','OtherHeads','outside'};
+        'samebox',{0,1};
+        };
+end
+
+cfgFit = [];
+cfgFit.precondition = 1;
+cfgFit.lsmriterations = 1500;
+% if subset of channels is wanted:
+% elecfind = @(x)find(strcmp(x,{EEG.chanlocs.labels}));
+% cfgFit.channel = [elecfind('PO8')];
+cfgFit.channel  = 1:length(EEG.chanlocs);
+if cfg.regularize
+    cfgFit.method='glmnet';
+    cfgFit.glmnetalpha = 0;
+end
+
+EEG = uf_designmat(EEG,cfgDesign);
+
+% if we want a mashup design, exchange everything <1 with 0, in all
+% uninteresting parameters
+
+if cfgDesign.mashup
+    parampos=find(ismember(EEG.unfold.colnames,{'curr_OtherHeads','prev_OtherHeads','curr_outside','prev_outside'}));
+    
+    for p=1:length(parampos)
+        pos=find(EEG.unfold.X(:,parampos(p))<0);
+        EEG.unfold.X(pos,parampos(p))=0;
+    end
+    
+    sameboxpos=find(strncmp(EEG.unfold.colnames,'same',4));
+    sbpos0=find(EEG.unfold.X(:,sameboxpos)<0);
+    EEG.unfold.X(sbpos0,sameboxpos)=0;
+    
+    sameboxnpos=find(strncmp(EEG.unfold.colnames,'sameboxn',8));
+    psbpos0=find(EEG.unfold.X(:,sameboxnpos)<0);
+    EEG.unfold.X(psbpos0,sameboxnpos)=0;
+    
+    %the interaction is basically only the multiplicatiom of curr and prev
+    interactpos=find(ismember(EEG.unfold.colnames,{'humanface:prevhumanface'}));
+    mainpos =find(ismember(EEG.unfold.colnames,{'curr_HF','prev_HF'}));
+    EEG.unfold.X(:,interactpos)=EEG.unfold.X(:,mainpos(1)).*EEG.unfold.X(:,mainpos(2));
+end
+
+
+%% Timeshift the design matrix
+for t=1:size(cfg.timewin,1)
+    cfgTimeshift = [];
+    cfgTimeshift.timelimits = cfg.timewin(t,:);
+    
+    EEG = uf_timeexpandDesignmat(EEG,cfgTimeshift);
+    
+    EEG = uf_continuousArtifactExclude(EEG,'winrej',cfg.cleantimes);
+    
+    
+    %% Fit it iteratively
+    EEG= uf_glmfit(EEG,cfgFit);
+    
+    
+    %% Make a massive uni-variate fit without de-convolution
+    EEGepoch = uf_epoch(EEG,'winrej',cfg.cleantimes,'timelimits',cfgTimeshift.timelimits);
+    EEGepoch= uf_glmfit_nodc(EEGepoch);
+    
+    
+    %% Save the data
+    
+    ufresult = uf_condense(EEGepoch);
+    
+    save(['/net/store/nbp/projects/IntoTheWild/Daten/EEG/WLFO/unfold_WLFO/' cfg.filtering '/ufresult_subj' num2str(sub_num) '_allElec_mashup_' num2str(cfg.timewin(t,1)*1000) '-' num2str(cfg.timewin(t,2)*1000) '_regularized' num2str(cfg.regularize) '.mat'],'ufresult')
+    
+    disp('successfully saved')
+end
+
+

code/unfold/unfold_lab.m

@@ -0,0 +1,378 @@
+function [ufresult]=unfold_lab(sub_num)
+%%
+cfg = struct();
+cfg.subject = num2str(sub_num);
+cfg.filtering='causal';
+cfg.mainfolder = ['/net/store/nbp/projects/IntoTheWild/Daten/EEG/Lab/VP' cfg.subject '/preprocessed/' cfg.filtering];
+tmp = load(fullfile(cfg.mainfolder,['4_ITW_Lab_subj' cfg.subject '_cleaningTimes.mat']));
+cfg.cleantimes = tmp.rej(:,1:2);
+
+cfg.file = ['3_ITW_Lab_subj' cfg.subject '_channelrejTriggersXensor.set'];
+tmp = pop_loadset('filename',['2_ITW_Lab_subj' cfg.subject '_bandpass_resample_deblank.set'],'filepath',cfg.mainfolder,'loadmode','info');
+tmp.chanlocs(find(strcmp({tmp.chanlocs.labels},'CzAmp2'))).labels = 'Iz';
+tmp=pop_chanedit(tmp, 'lookup','/net/store/nbp/projects/IntoTheWild/EEG_analysis/eeglab13_5_4b/plugins/dipfit2.3/standard_BESA/standard-10-5-cap385.elp');
+cfg.urchanlocs= tmp.chanlocs;
+cfg.avgref=1;
+cfg.saccade=1;
+cfg.regularize=0;
+tmp = load(fullfile(cfg.mainfolder,['6_ITW_Lab_subj' cfg.subject '_ICAcleancont.mat']));
+cfg.badcomponents = tmp.comps_to_rej;
+
+cfg.timewin=[-0.5 1];
+
+
+%% load EEG data
+EEG = pop_loadset('filename',cfg.file,'filepath',cfg.mainfolder);
+EEG =  pop_select(EEG,'nochannel',{'AUX1','AUX2'});
+
+
+%% Load ICA and clean
+addpath('/net/store/nbp/projects/EEG/blind_spot/amica')
+mod = loadmodout12(fullfile(cfg.mainfolder,'amica'));
+EEG.icaweights = mod.W;
+EEG.icasphere = mod.S;
+EEG.icawinv = [];EEG.icaact = [];EEG.icachansind = [];
+EEG = eeg_checkset(EEG);
+EEG = pop_subcomp(EEG,cfg.badcomponents);
+
+
+%% interpolate missing channels
+% ET channels (Gaze etc.) not interpolated
+alldel = {'CzAmp2' 'BIP1' 'BIP2' 'BIP3' 'BIP4' 'BIP5' 'BIP6' 'BIP7' 'BIP8' 'AUX1' 'AUX2' 'AUX3' 'AUX4' 'AUX5' 'AUX6' 'AUX7' 'AUX8','TIME', 'L_GAZE_X', 'L_GAZE_Y', 'L_AREA','R_GAZE_X', 'R_GAZE_Y', 'R_AREA', 'INPUT', 'L-GAZE-X', 'L-GAZE-Y', 'L-AREA','R-GAZE-X', 'R-GAZE-Y', 'R-AREA', 'INPUT'};
+
+%remove VEOG from EEG, as there is no corresponding location
+EEG = pop_select(EEG, 'nochannel', {'VEOG'});
+
+if cfg.avgref
+    %calculate avg ref
+    EEG = pop_reref( EEG, []); %Participants’ averages were then re-referenced to a common average reference. (Rossion & Caharel, 2011)
+end
+
+idxs = ~ismember({cfg.urchanlocs.labels},alldel);
+
+%interpolate missing channels
+EEG= pop_interp(EEG,cfg.urchanlocs(idxs),'spherical');
+
+
+%% adjust for trigger delay
+% as calculated from photosensor
+if sub_num<36
+    trigger_delay=0.0102; %10.2ms
+else
+    trigger_delay=0.0348; %34.8ms
+end
+
+for evt=1:length(EEG.event)
+    if str2num(EEG.event(evt).type)<=5 %only for stim onset [1 2 3 4 5]
+        assert(str2num(EEG.event(evt).type)>0)
+        EEG.event(evt).latency = EEG.event(evt).latency + (trigger_delay*EEG.srate);
+    end
+end
+EEG = eeg_checkset(EEG);
+
+%% adjust cleaning times to be between 1 and length of data
+[a,b]=find(cfg.cleantimes(:,1)<1);
+if ~isempty(a)
+    cfg.cleantimes(a,1)=1;
+end
+
+[a,b]=find(cfg.cleantimes(:,1)>length(EEG.data));
+if ~isempty(a)
+    cfg.cleantimes(a,:)=[];
+end
+
+
+[a,b]=find(cfg.cleantimes(:,2)>length(EEG.data));
+if ~isempty(a)
+    cfg.cleantimes(a,2)=length(EEG.data);
+end
+numstim=0;
+%% for each fixation we need to save ALL the necessary information in the EEG structure
+% for e = 1:length(EEG.event)
+for e = length(EEG.event):-1:1
+    
+    % 1-5 stimulus
+    % 180 stim offset
+    % 200 exp start
+    % 255 exp stop
+    
+    % if the trigger is a fixation
+    if ismember(EEG.event(e).type, ...
+            {'L_fixation' 'R_fixation'})
+
+        EEG.event(e).urtype = EEG.event(e).type;
+        
+        % if it is a meaningful fixation
+        
+        EEG.event(e).type = 'fixation';
+        
+        continue
+
+    end
+
+    % if the trigger was the beginning of a trial
+    if str2num(EEG.event(e).type) <= 5
+        EEG.event(e).urtype = EEG.event(e).type;
+        EEG.event(e).type = 'stimulus';
+    end
+    
+    
+    
+    % recode all fix that are not Face or None to "OtherHeads"
+    if strcmp(EEG.event(e).type, 'stimulus')
+        isStim=1;
+    else
+        isStim=0;
+    end
+    
+    
+    
+    if isStim==1
+        
+        numstim=numstim+1;
+        
+        if EEG.event(e).humanface==1
+            EEG.event(e).curr='HF';
+            
+        elseif EEG.event(e).car==1
+            EEG.event(e).curr='Car';
+        elseif EEG.event(e).obj==1
+            EEG.event(e).curr='Obj';
+        elseif EEG.event(e).none==1
+            EEG.event(e).curr='None';
+        else
+            %first in a trial has no prevprev, therefore we define it as
+            %'outside' as outside will not be analysed (easiest solution)
+            EEG.event(e).curr='outside';
+        end
+        
+        if EEG.event(e).prevhumanface==1
+            EEG.event(e).prev='HF';
+        elseif EEG.event(e).prevcar==1
+            EEG.event(e).prev='Car';
+        elseif EEG.event(e).prevobj==1
+            EEG.event(e).prev='Obj';
+        elseif EEG.event(e).prevnone==1
+            EEG.event(e).prev='BlockBeginning';
+            
+        end
+        
+        
+        if EEG.event(e).prevprevhumanface==1
+            EEG.event(e).prevprev='HF';
+        elseif EEG.event(e).prevprevnone==1
+            EEG.event(e).prevprev='None';
+        else
+            %first in a trial has no prevprev, therefore we define it as
+            %'outside' as outside will not be analysed (easiest solution)
+            EEG.event(e).prevprev='BlockBeginning';
+        end
+        
+    end
+  
+end
+
+
+%% save information about saccaded and fixations
+if sum(strcmp({EEG.event.type},'L_saccade'))>0
+    ixSac = find(strcmp({EEG.event.type},'L_saccade'));
+else
+    ixSac = find(strcmp({EEG.event.type},'R_saccade'));
+end
+
+ixFix = find(strcmp({EEG.event.type},'fixation'));
+
+%find blinks
+ixBlk=find(cellfun(@(x)strcmp(x,'R_blink'),{EEG.event.type}));
+if isempty(ixBlk)
+    ixBlk=find(cellfun(@(x)strcmp(x,'L_blink'),{EEG.event.type}));
+end
+
+
+latSac = [EEG.event(ixSac).latency];
+latFix = [EEG.event(ixFix).latency];
+latBlk = [EEG.event(ixBlk).latency];
+
+endSac = latSac+[EEG.event(ixSac).duration];
+endBlk = latBlk+[EEG.event(ixBlk).duration];
+
+% put the velocity for blink sac/data loss sections to NaN
+for ixS=1:length(ixSac)
+    
+    blkStarAftertSacStart=find(latBlk>=latSac(ixS),1,'first');
+    blkStopBeforeSacStop=find(endBlk<=endSac(ixS),1,'last');
+    if blkStarAftertSacStart== blkStopBeforeSacStop
+        EEG.event(ixSac(ixS)).sac_vmax=NaN;
+    end
+end
+
+
+
+for e = 1:length(latFix)
+    % for each fixation, find the closest saccade before the event and copy
+    % the saccade information
+    
+    % when  Sac --> Fix, then lat(Fix) - lat(Sac) should be positive
+    ix = find((latFix(e) - latSac) > 0,1,'last');
+    
+    if ~isempty(ix)
+        for fn = {'sac_endpos_x','sac_endpos_y','sac_amplitude','sac_startpos_x','sac_startpos_y','sac_vmax'}
+            EEG.event(ixFix(e)).(fn{1})= EEG.event(ixSac(ix)).(fn{1});
+        end
+    end
+    
+end
+
+% remove all fixation with saccade velocity of 0 (= blink saccades)
+for ixF=1:length(ixFix)
+    if EEG.event(ixFix(ixF)).sac_vmax==0
+        EEG.event(ixFix(ixF)).type='blinkFix';
+    end
+end
+
+ixFix = find(strcmp({EEG.event.type},'fixation'));
+latFix =  [EEG.event(ixFix).latency];
+
+% remove all fixations close to a blink/data loss section
+if ~isempty(ixBlk)
+    ixBlk=fliplr(ixBlk);
+    for iB=ixBlk
+        blkStop=EEG.event(iB).latency+EEG.event(iB).duration;
+        delFix=find(abs(latFix-blkStop)<50);
+        
+        if ~isempty(delFix)
+            for dF=1:length(delFix)
+                EEG.event(ixFix(delFix(dF))).type='blinkFix';
+            end
+        end
+    end
+end
+
+for e = 1:length(latSac)
+    % for each fixation, find the closest saccade before the event and copy
+    % the saccade information
+    
+    % when  Sac --> Fix, then lat(Fix) - lat(Sac) should be positive
+    ix = find((latSac(e) - latFix) > 0,1,'last');
+    if ~isempty(ix) && ~isempty( EEG.event(ixFix(ix)).prev)
+        for fn = {'prev','curr'}
+            EEG.event(ixSac(e)).(fn{1})= EEG.event(ixFix(ix)).(fn{1});
+            
+        end
+        EEG.event(ixSac(e)).type='saccade';
+    end
+    
+end
+
+
+%% Remove large saccades
+fprintf('\n\n %u large saccades (>40) will be removed. \n\n\n', length(find([EEG.event.sac_amplitude]>40)))
+EEG.event(cellfun(@(x)(x>40),{EEG.event.sac_amplitude})) = [];
+
+
+
+%% Config
+EEG.unfold = [];
+cfgDesign = [];
+
+cfgDesign.splinespacing = 'quantile';
+
+cfgDesign.eventtypes = {
+    {'fixation'}   % fixation onset
+    {'stimulus'}}; % stimulus onset
+
+if sub_num==51 %51 has only 2 relevant sac before fix - maybe due to stopping recording between trials
+    cfgDesign.formula = {
+        'y~1 + spl(fix_avgpos_x,5) + spl(fix_avgpos_y,5)'
+        'y~1 + cat(curr) + cat(prev)+ humanface:prevhumanface'};
+else
+    cfgDesign.formula = {
+        'y~1 + spl(fix_avgpos_x,5) + spl(fix_avgpos_y,5) +  spl(sac_amplitude,5)'
+        'y~1 + cat(curr) + cat(prev)+ humanface:prevhumanface'};
+end
+
+
+cfgDesign.codingschema = 'effects';
+cfgDesign.mashup=1;
+
+% if codingschema=effects -> reference cat is the last one
+% if codingschema=reference -> reference cat is the first one
+if strcmp(cfgDesign.codingschema, 'effects')
+    cfgDesign.categorical = {
+        'curr',{'Car','HF','Obj'};
+        'prev',{'BlockBeginning','Car','HF','Obj'};
+        };
+elseif strcmp(cfgDesign.codingschema, 'reference')
+    cfgDesign.categorical = {
+        'curr',{'Obj','HF','Car','BlockBeginning'};
+        'prev',{'Obj','HF','Car','BlockBeginning'};
+        };
+end
+
+cfgFit = [];
+cfgFit.precondition = 1;
+cfgFit.lsmriterations = 1500;
+% if subset of channels is wanted:
+% elecfind = @(x)find(strcmp(x,{EEG.chanlocs.labels}));
+% cfgFit.channel = [elecfind('PO8')];
+cfgFit.channel  = 1:length(EEG.chanlocs);
+if cfg.regularize
+    cfgFit.method='glmnet';
+    cfgFit.glmnetalpha = 0;
+end
+
+EEG = uf_designmat(EEG,cfgDesign);
+
+% if we want a mashup design, exchange everything <1 with 0, in all
+% uninteresting parameters
+if cfgDesign.mashup
+    parampos=find(ismember(EEG.unfold.colnames,{'prev_Car','curr_Car','prev_BlockBeginning'}));
+    
+    for p=1:length(parampos)
+        pos=find(EEG.unfold.X(:,parampos(p))<0);
+        EEG.unfold.X(pos,parampos(p))=0;
+    end
+    
+    %the interaction is basically only the multiplicatiom of curr and prev
+    interactpos=find(ismember(EEG.unfold.colnames,{'humanface:prevhumanface'}));
+    mainpos =find(ismember(EEG.unfold.colnames,{'curr_HF','prev_HF'}));
+    EEG.unfold.X(:,interactpos)=EEG.unfold.X(:,mainpos(1)).*EEG.unfold.X(:,mainpos(2));
+end
+
+
+%% Add cyclical regressor
+% if cfgDesign.cyclicalAngle
+%     EEG = scenes_add_cyclical(EEG);
+% end
+
+
+%% Timeshift the design matrix
+for t=1:size(cfg.timewin,1)
+    cfgTimeshift = [];
+    cfgTimeshift.timelimits = cfg.timewin(t,:);
+    
+    EEG = uf_timeexpandDesignmat(EEG,cfgTimeshift);
+    
+    EEG = uf_continuousArtifactExclude(EEG,'winrej',cfg.cleantimes);
+    
+    
+    %% Fit it iteratively
+    EEG= uf_glmfit(EEG,cfgFit);
+ 
+    
+    %% Make a massive uni-variate fit without de-convolution
+    EEGepoch = uf_epoch(EEG,'winrej',cfg.cleantimes,'timelimits',cfgTimeshift.timelimits);
+    EEGepoch= uf_glmfit_nodc(EEGepoch);
+
+    
+    %% Save the data
+    
+    ufresult = uf_condense(EEGepoch);
+    
+    save(['/net/store/nbp/projects/IntoTheWild/Daten/EEG/Lab/unfold_Lab/' cfg.filtering '/ufresult_subj' num2str(sub_num) '_allElec_mashup_' num2str(cfg.timewin(t,1)*1000) '-' num2str(cfg.timewin(t,2)*1000) '_regularized' num2str(cfg.regularize) '.mat'],'ufresult')
+    
+    
+    disp('successfully saved')
+end
+
+