function [stats, sdat] = multi_group_compare(data,alpha,group,forcenonpara)
% STATS = MULTI_GROUP_COMPARE(DATA,ALPHA,GROUP) Performs a multiple
% comparison test on the data matrix or cell array at significance level
% ALPHA. If DATA is a matrix, each column contains data from a particular
% group. If DATA is a cell array, each cell contains data for a particular
% group. GROUP is a char array that contains the names for each group (one
% per cell/column of DATA).
%
% [STATS, DMAT] = MULTI_GROUP_COMPARE(DATA,ALPHA,GROUP) Returns the matrix
% version of the input data with missing values replaced with nans.

% Reformat data as matrix, even if there are different numbers of samples
% in each group
if iscell(data)
    n = cellfun(@numel,data);
    sdat = nan(max(n),numel(data));
    for i = 1:numel(n)
        sdat(1:n(i),i) = data{i};
    end
else
    sdat = data;
end

ngroup = size(sdat,2);

% Default parameters
if nargin < 4 || isempty(forcenonpara)
    forcenonpara = 0;
end

if nargin < 3 || isempty(group)
    group = (1:ngroup)';
    group = num2str(group);
end
stats.groups = group;

if nargin < 2 || isempty(alpha)
    alpha = 0.05;
end

%% Test for Normality
[isnorm, h_norm, p_norm] = norm_test(sdat,alpha);

% Try to normalize data with a sqrt transform
if ~isnorm && sum(sum(sdat < 0)) == 0
    sdat_nrm = log(sdat);
    %     sdat_nrm = sqrt(sdat);
    [isnorm, h_norm, p_norm] = norm_test(sdat_nrm,alpha);
    if isnorm && ~forcenonpara
        normalized = 1;
    else
        normalized = 0;
    end
else
    sdat_nrm = sdat;
    normalized = 0;
end

% Gather statistical results
stats.norm.forcenp = forcenonpara;
stats.norm.isnorm = isnorm;
stats.normalized = normalized;
stats.norm.h_norm = h_norm;
stats.norm.p_norm = p_norm;

%% Multiple Comparison Testing
if isnorm && ~forcenonpara
    % Proceed with parametric testing
    [p,table,stat] = anova1(sdat_nrm,group,'off');
    stats.test.type = 'anova1';
    stats.test.p = p;
    stats.test.table = table;
    stats.test.stats = stat;
    
    % Multiple comparisons
    tests = combnk(1:ngroup,2);
    if tests(1) ~= 1
        tests = flipud(tests);
    end
    p = nan(size(tests,1),1);
    if stats.test.p < alpha
        for i = 1:size(tests,1)
            [~, p(i)] = ttest(sdat_nrm(:,tests(i,1)),sdat_nrm(:,tests(i,2)),alpha);
        end
    end
    
    % Bonferroni correction
    p = p*size(tests,1);
    p(p>1) = 1;
    h = p < alpha;
    stats.multcomp.pairs = tests;
    stats.multcomp.h = logical(h);
    stats.multcomp.p = p;
    
    [c,m] = multcompare(stat,'ctype','bonferroni','alpha',alpha,'display','off');
    c = [c ~(c(:,3) < 0 & c(:,5) > 0)]; % Are the groups significantly different?
    stats.multcomp.c = c;
    stats.multcomp.m = m;
    stats.multcomp.usem = [nanmean(sdat)' nanstd(sdat)'/sqrt(size(sdat,1))];
    
else
    % Proceed with non-parametric testing
    [p,table,stat] = kruskalwallis(sdat,group,'off');
    stats.test.type = 'kruskalwallis';
    stats.test.p = p;
    stats.test.table = table;
    stats.test.stats = stat;
    
    % Multiple comparisons
    tests = combnk(1:ngroup,2);
    if tests(1) ~= 1
        tests = flipud(tests);
    end
    
    p = nan(size(tests,1),1);
    if stats.test.p < alpha
        for i = 1:size(tests,1)
            p(i) = ranksum(sdat(:,tests(i,1)),sdat(:,tests(i,2)),'alpha', alpha);
        end
    end
    
    % Bonferroni correction
    p = p*size(tests,1);
    p(p>1) = 1;
    h = p < alpha;
    stats.multcomp.pairs = tests;
    stats.multcomp.h = logical(h);
    stats.multcomp.p = p;
    
    % Multiple comparison with Bonferroni correction
    [c,m] = multcompare(stat,'ctype','bonferroni','alpha',alpha,'display','off');
    c = [c ~(c(:,3) <= 0 & c(:,5) >= 0)]; % Are the groups significantly different?
    
    stats.multcomp.c = c;
    stats.multcomp.m = m;
    stats.multcomp.usem = [nanmean(sdat)' nanstd(sdat)'/sqrt(size(sdat,1))];
    
end

end


function [isnorm, h_norm, p_norm] = norm_test(data,alpha)


%% Test for Normality
h_norm = []; % reject normality hypothesis for each group?
p_norm = []; % p-value for each group
if alpha < 0.001
    alpha_n = 0.001;
else
    alpha_n = alpha;
end
for i = 1:size(data,2)
    [h, p] = lillietest(data(:,i),alpha_n);
    h_norm = [h_norm h];
    p_norm = [p_norm p];
end

% Do any of the normality null hypotheses get rejected?
isnorm = ~sum(h_norm);
end


