function plot_fwd_vs_bwd_GLME_visNet_deltaDI(sessDirectTransf,sel_time_window,params,computedROIsGroupsLab,sort_idxs,paths)
% Session-level stat pooling links

allGroups = unique(computedROIsGroupsLab(sort_idxs),'stable');

figBox=figure('Position',[360,198,348,420]);

xpos_increm = 0;
% maxYFIT = max(max([nanmean(sessDirectTransf.all.FIT_S.(params.band_label).allROIs.forward);nanmean(sessDirectTransf.all.FIT_S.(params.band_label).allROIs.feedback);...
%     nanmean(sessDirectTransf.all.FIT_C.(params.band_label).allROIs.forward);nanmean(sessDirectTransf.all.FIT_C.(params.band_label).allROIs.feedback)]));
maxYFIT = 0;

for infoIdx = 1:numel(params.info_type)
    infoLab = params.info_type{infoIdx};
    
    % Pick feedforward and feedback links
    measure1 = sessDirectTransf.all.(infoLab).(params.band_label).allROIs.forward(:)';
    measure2 = sessDirectTransf.all.(infoLab).(params.band_label).allROIs.feedback(:)';

%     measure1 = nanmean(sessDirectTransf.all.(infoLab).(params.band_label).allROIs.forward,1);
%     measure2 = nanmean(sessDirectTransf.all.(infoLab).(params.band_label).allROIs.feedback,1);
    
    % Compute statistics
    if params.doGLMEstat
        subjMeas1 = sessDirectTransf.all.(infoLab).(params.band_label).allROIs.forward(:); % (:) sorts the elements by colums first (i.e. we have that first values are all from subj 1)
        subjMeas2 = sessDirectTransf.all.(infoLab).(params.band_label).allROIs.feedback(:);
        fwd_fbk_idx = [zeros(1,numel(subjMeas1)),ones(1,numel(subjMeas2))];
        subj_groups = 1:params.nSubj;
        %subj_groups = repelem(subj_groups,1,params.nSess);
        subj_groups = repelem(subj_groups,1,params.nSess*size(sessDirectTransf.all.(infoLab).(params.band_label).allROIs.forward,2));
        subj_groups = [subj_groups,subj_groups]';

        measures=[subjMeas1;subjMeas2];
        tbl = table(measures,subj_groups,fwd_fbk_idx','VariableNames',{'measures','Subj','Cond'});
        model = fitglme(tbl,'Cond~measures + (1|Subj)', 'Distribution', 'Binomial');
        p = model.Coefficients.pValue(2);
    else
        [a,p] = ttest(measure1,measure2);
    end
    
    % Boxplot
    figure(figBox)
    hold on
    mean1 = mean(measure1);
    sem1 = std(measure1)/sqrt(numel(measure1));
    mean2 = mean(measure2);
    sem2 = std(measure2)/sqrt(numel(measure2));

    scatter(1,mean1,12,'r','filled')
    scatter(2,mean2,12,'b','filled')
    errorbar(1,mean1,sem1,'r')
    errorbar(2,mean2,sem2,'b')
    
    maxY = max([mean1+sem1,mean2+sem2]);
    if infoIdx > 1
        maxYFIT = max([maxYFIT,maxY]);
        maxY = maxYFIT;
    end
    
    ylim([0,(3/2)*maxY])
    xlim([0.5,2.5])
    set(gca, 'XTick', [1 2], 'XTickLabel', {'fwd','fbk'}, 'fontsize', 12)
    pvalues_plot(p,1.5,(5/4)*maxY,maxY/10,maxY/10,14,0.5,'k',0)
    
end
title(['Fwd vs Fbk deltaDI T = [',num2str(sel_time_window(1)),',',num2str(sel_time_window(2)),'] ', params.bands_labels{params.bandIdx}])
if params.save_figures
    fnamePNG = [paths.figuresPath,'\allROIs_fwd_vs_fbk_deltaDI_ttest_',date,'.png'];
    fnameFIG = [paths.figuresPath,'\allROIs_fwd_vs_fbk_deltaDI_ttest_',date,'.fig'];
    fnameSVG = [paths.figuresPath,'\allROIs_fwd_vs_fbk_deltaDI_ttest_',date,'.svg'];
    saveas(figBox,fnamePNG)
    saveas(figBox,fnameFIG)
    figBox.Renderer = 'painters';
    saveas(figBox,fnameSVG)
end

end

