Fclust_corr.m

%Calculate cluster correction for a F-observed and the permutation
%distribution
%
%EXAMPLE USAGE
% >> [h, p, clust_info, est_alpha] = Fclust_corr(F_obs, F_dist, 0.05, chan_hood, 3.84)
%
%REQUIRED INPUTS
% F_obs          - An electrode x time points array of observed F-values
% F-dist         - A permutation x electrode x time point array of the
%                  permutation F-distribution
% alpha          - Alpha level to use for hypothesis test
% chan_hood      - A symmetric 2d matrix indicating which channels are
%                  neighbors with other channels.  If chan_hood(a,b)=1,
%                  then Channel A and B are neighbors. This is produced by the
%                  function spatial_neighbors.m.
% thresh_F       - F-statistic threshold for inclusion in cluster
%
%OUTPUT
% h              - electrode x time point array indicating which locations
%                  are part of a statistically significant cluster
% p              - electrode x time point array of p-values
% clust_info     - struct containing information about each cluster
% est_alpha      - estimated achieved alpha level of the test; may not be
%                  accurate if F_dist was generated by an approximate
%                  permutation method
%
%
%VERSION DATE: 4 March 2020
%AUTHOR: Eric Fields
%
%NOTE: This function is provided "as is" and any express or implied warranties 
%are disclaimed. 

%Copyright (c) 2017-2019, Eric Fields
%All rights reserved.
%This code is free and open source software made available under the 3-clause BSD license.

function [h, p, clust_info, est_alpha] = Fclust_corr(F_obs, F_dist, alpha, chan_hood, thresh_F)

    global VERBLEVEL
    if VERBLEVEL
        fprintf('Calculating clusters . . . ')
    end
    
    if ~isequal(chan_hood, chan_hood')
        error('Your chan_hood matrix is not symmetric')
    end
    
    %Some useful numbers
    [n_perm, n_electrodes, n_time_pts] = size(F_dist);
    
    %Find clusters and cluster mass distribution
    clust_mass_dist = zeros(n_perm, 1);
    for i = 1:n_perm
        F = reshape(F_dist(i, :, :), [n_electrodes, n_time_pts]);
        %Find clusters in F
        [clust_ids, n_clust] = find_clusters(F, thresh_F, chan_hood, 1);
        %Save observed data
        if i == 1
            clust_ids_obs = clust_ids;
            n_clust_obs = n_clust;
        end
        %Find largest cluster mass for this permutation
        for c = 1:n_clust
            use_mass = sum(F(clust_ids == c));
            if use_mass > clust_mass_dist(i)
                clust_mass_dist(i) = use_mass;
            end
        end 
    end
    
    %Get cluster mass critical value
    clust_mass_dist = sort(clust_mass_dist);
    clust_mass_crit = clust_mass_dist(ceil((1-alpha) * n_perm));
    
    %Find cluster mass for each observed cluster and compare to critical
    %value
    clust_mass_obs = NaN(1, n_clust_obs);
    clust_pval = NaN(1, n_clust_obs);
    null_test = false(1, n_clust_obs);
    p = ones(n_electrodes, n_time_pts);
    for c = 1:n_clust_obs
        clust_mass_obs(c) = sum(F_obs(clust_ids_obs == c));
        if clust_mass_obs(c) > clust_mass_crit
            null_test(c) = true;
        end
        clust_pval(c) = mean(clust_mass_dist >= clust_mass_obs(c));
        p(clust_ids_obs == c) = clust_pval(c);
    end
    h = (p <= alpha);
    est_alpha = mean(clust_mass_dist>clust_mass_crit);
    if VERBLEVEL
        fprintf('DONE\n');
        fprintf('Estimated alpha level is %f\n', est_alpha);
    end
    
    clust_info = struct('null_test', null_test, ...
                        'pval', clust_pval, ...
                        'clust_mass', clust_mass_obs, ...
                        'clust_ids', clust_ids_obs);

end