%% Second level statistics - across subjects
%specify: ROI name
% testType : 'memory','category','complexity','memorability'


function second_level_stats(subjects,testType,nrep,rsa_tag,roi_tag,stats_tag,dir_tag)
    
    fs = '\';
    baseDir = '..\Data';
    
    %load folder names of subjects
    load('..\subject_list.mat'); %subNames for easier folder accessing
    
    
    %% load ROI first-level stats for each subject and store empirical and permutation t-values in one struct
    groupstats = {};
    groupstats.subnum = [];
    groupstats.channum = [];
    for subject = 1:length(subjects)
        first_lvl_stats_dir = [baseDir fs subNames{subjects(subject)} fs 'stats' fs 'ERS' fs rsa_tag fs roi_tag fs stats_tag];
        
        load([first_lvl_stats_dir fs 'first_level_stats.mat'])
            
        %add subnum and channels to groupstats file
        groupstats.subnum(end+1,1) = subjects(subject);
        groupstats.channum(end+1,1) = length(stats.channels);
        
        if strcmp(testType,'memory')
            groupstats.sub.emp(subject,:,:) = stats.emp.memory(:,:,5); %idx 5 is t-statistics (mean(rem), mean(forg), h, p, t)
            groupstats.sub.perm(subject,:,:,:) = squeeze(stats.perm.memory(:,:,:,5));
        elseif strcmp(testType,'category')
            groupstats.sub.emp(subject,:,:) = stats.emp.category(:,:,5);
            groupstats.sub.perm(subject,:,:,:) = squeeze(stats.perm.category(:,:,:,5));
        elseif strcmp(testType,'complexity')
            groupstats.sub.emp(subject,:,:) = stats.emp.complexity(:,:,5);
            groupstats.sub.perm(subject,:,:,:) = squeeze(stats.perm.complexity(:,:,:,5));
        elseif strcmp(testType,'memorability')
            groupstats.sub.emp(subject,:,:) = stats.emp.memorability(:,:,5);
            groupstats.sub.perm(subject,:,:,:) = squeeze(stats.perm.memorability(:,:,:,5));
        end
    end
        
    %% perform second-level ttests across subjects  
    for win1 = 1:size(groupstats.sub.emp,2) %encoding time
        for win2 = 1:size(groupstats.sub.emp,3) %recognition time
            [~,P,~,STATS] = ttest(squeeze(groupstats.sub.emp(:,win1,win2)));
            
            groupstats.emp.t(win1,win2) = STATS.tstat;
            groupstats.emp.p(win1,win2) = P;
            
            for rep = 1:nrep
                [~,P,~,STATS] = ttest(squeeze(groupstats.sub.perm(:,rep,win1,win2)));
            
                groupstats.perm.t(rep,win1,win2) = STATS.tstat;
                groupstats.perm.p(rep,win1,win2) = P;
            end
        end
    end
            
    %% Cluster-statistics
    %empirical clusters
    tmap = groupstats.emp.t;
    pmap = groupstats.emp.p;
    groupstats.emp_cluster = extract_clusterstats(tmap,pmap);
    
    %permutation clusters
    for rep = 1:nrep
        tmap = squeeze(groupstats.perm.t(rep,:,:));
        pmap = squeeze(groupstats.perm.p(rep,:,:));
        
        groupstats.perm_cluster(rep) = extract_clusterstats(tmap,pmap);
        
        perm_pos_cluster(rep) = groupstats.perm_cluster(rep).pos.sum;
        perm_neg_cluster(rep) = groupstats.perm_cluster(rep).neg.sum;
    end
    
    % permutation cluster size distribution and ranking of empirical in
    % surrogate data
    perm_pos_dist = sort(perm_pos_cluster,'ascend');
    perm_neg_dist = sort(perm_neg_cluster,'descend');
    
    groupstats.emp_cluster.pos.rank = nearest(perm_pos_dist, groupstats.emp_cluster.pos.sum);
    groupstats.emp_cluster.pos.p = 1 - (groupstats.emp_cluster.pos.rank / nrep);
    
    groupstats.emp_cluster.neg.rank = nearest(perm_neg_dist, groupstats.emp_cluster.neg.sum);
    groupstats.emp_cluster.neg.p = 1 - (groupstats.emp_cluster.neg.rank / nrep);

%     %% plot clusterstats results
%     f = plot_clusterstats(groupstats,rsawidth,rsastep,maxT,ROI,testType);
    
    %% save stats 
    second_level_stats_dir = ['..\STATS' fs dir_tag fs rsa_tag fs stats_tag fs testType];
    status = mkdir(second_level_stats_dir);

    save([second_level_stats_dir fs roi_tag '_stats.mat'],'groupstats');
end

    
    
            
            
            
            
        
        