function stat_ds = cosmo_stat(ds, stat_name, output_stat_name)
% compute t-test or F-test (ANOVA) statistic
%
% stat_ds=cosmo_stats(ds, stat_name[, output_stat_name])
%
% Inputs:
% ds dataset struct with
% .samples PxQ, for P observations on Q features
% .sa.targets Px1 observation conditions (classes)
% .sa.chunks Px1 observation chunks (e.g. subjects).
% stat_name One of [*]:
% 't' : one-sample t-test against zero (nclasses==1),
% or paired t-test (nclasses==2)
% 't2': two-sample t-test with equal variance,
% (nclasses==2)
% 'F' : one-way ANOVA or repeated measures ANOVA
% (nclasses>=2)
% [*] nclasses is the number of unique values in
% ds.sa.targets
% [**] if nclasses==2, the contrast of unq(1) minus
% unq(2) is returned, where
% unq=unique(ds.sa.targets);
% Because unique returns the results in sorted
% order from small to large, this means that
% the contrast is computed of the smaller value
% of targets minus the larger. This may be somewhat
% confusing, in the sense that the output statistics
% (t- or z-score) gives a positive value if samples
% associated with the lower .sa.targets value have a
% higher mean than those associated with the
% higher .sa.targets value.
% output_stat_name (optional) 'left', 'right', 'both', 'z', 'p', or
% the empty string '' (default).
% - 'left', 'right', and 'both' return a p-value with
% the specified tail.
% - 'p' returns a p-value, with tail='right' if
% stat_name='F' and tail='both' otherwise.
% - 'z' returns a z-score corresponding to the p-value
% - '' (empty) returns the t or F statistic.
% Missing values can be indicated by NaNs; if these are
% present and a p-value or z-score is returned, then
% these values are computed with a possible variable
% number of degrees of freedom across features.
%
% Returns:
% stat_ds dataset struct with fields:
% .samples 1xQ statistic value, or (if output_stat_name is
% non-empty) z-score or p-value. See the Notes below
% for interpreting p-values.
% .sa struct with field:
% .stats One of 'Ftest(df1,df2)', 'Ttest(df)', 'Zscore()', or
% 'Pval()', where df* are the degrees of freedom
% .[f]a identical to the input ds.[f]a, if present.
%
% Examples:
% % one-sample t-test
% % make a simple dataset
% ds=struct();
% ds.samples=reshape(mod(1:7:(12*3*7),13)',[],3)-3;
% ds.sa.targets=ones(12,1);
% ds.sa.chunks=(1:12)';
% cosmo_disp(ds.samples);
% %|| [ -2 4 -3
% %|| 5 -2 4
% %|| -1 5 -2
% %|| : : :
% %|| 9 2 8
% %|| 3 9 2
% %|| -3 3 9 ]@12x3
% %
% % run one-sample t-test
% s=cosmo_stat(ds,'t');
% cosmo_disp(s);
% %|| .samples
% %|| [ 2.49 3.36 2.55 ]
% %|| .sa
% %|| .stats
% %|| { 'Ttest(11)' }
% %
% % compute z-score of t-test
% s=cosmo_stat(ds,'t','z');
% cosmo_disp(s);
% %|| .samples
% %|| [ 2.17 2.73 2.21 ]
% %|| .sa
% %|| .stats
% %|| { 'Zscore()' }
% %
% % compute (two-tailed) p-value of t-test
% s=cosmo_stat(ds,'t','p');
% cosmo_disp(s);
% %|| .samples
% %|| [ 0.03 0.00633 0.0268 ]
% %|| .sa
% %|| .stats
% %|| { 'Pval()' }
% %
% % compute left-tailed p-value of t-test
% s=cosmo_stat(ds,'t','left');
% cosmo_disp(s);
% %|| .samples
% %|| [ 0.985 0.997 0.987 ]
% %|| .sa
% %|| .stats
% %|| { 'Pval()' }
%
% % one-way ANOVA
% % each observation is independent and thus each chunk is unique;
% % there are three conditions with four observations per condition
% ds=struct();
% ds.samples=reshape(mod(1:7:(12*3*7),13)',[],3)-3;
% ds.sa.targets=repmat(1:3,1,4)';
% ds.sa.chunks=(1:12)';
% s=cosmo_stat(ds,'F');
% cosmo_disp(s);
% %|| .samples
% %|| [ 0.472 0.0638 0.05 ]
% %|| .sa
% %|| .stats
% %|| { 'Ftest(2,9)' }
% % compute z-score
% s=cosmo_stat(ds,'F','z'); % convert to z-score
% cosmo_disp(s);
% %|| .samples
% %|| [ -0.354 -1.54 -1.66 ]
% %|| .sa
% %|| .stats
% %|| { 'Zscore()' }
%
%
% % two-sample t-test
% % each observation is independent and thus each chunk is unique;
% % there are two conditions with four observations per condition
% ds=struct();
% ds.samples=reshape(mod(1:7:(12*3*7),13)',[],3)-3;
% ds.sa.targets=repmat(1:2,1,6)';
% ds.sa.chunks=(1:12)';
% s=cosmo_stat(ds,'t2','p'); % return p-value
% cosmo_disp(s);
% %|| .samples
% %|| [ 0.0307 0.000242 7.07e-05 ]
% %|| .sa
% %|| .stats
% %|| { 'Pval()' }
% %
% % for illustration, this test gives the same p-values as a
% % repeated measures ANOVA
% s=cosmo_stat(ds,'F','p');
% cosmo_disp(s);
% %|| .samples
% %|| [ 0.0307 0.000242 7.07e-05 ]
% %|| .sa
% %|| .stats
% %|| { 'Pval()' }
%
% Notes:
% - If output_stat_name is not provided or empty, then this function runs
% considerably faster than the builtin matlab functions
% (ttest, ttest2, or anova1).
% - In the case of two classes, the sign of the result may seem counter-
% intuitive (see the [**]) comment above).
% - When output_stat_name=='p' then the p-values returned are the same as
% the builtin matlab functions anova1, ttest, and ttest2 with the
% default tails.
% - To run a one-sample t-tests against x (if x~=0), one has to
% subtract x from ds.samples before using ds as input to this function
% - The .sa.chunks and .sa.targets determine which test is performed:
% * statname=='t': all chunks are unique => one-sample t-test
% : each chunk present twice => paired-sample t-test
% * statname=='F': all chunks are unique => one-way ANOVA
% : each chunk present N times => repeated measures ANOVA
% See cosmo_montecarlo_cluster_stat for examples on how .sa.targets and
% .sa.chunks should be set for different statistics.
% - Missing values can be indicated by NaNs, and if the output is a
% p-value (or a z-score based on the p-value), then the output is
% computed for different features using varying degrees of freedom.
% For example, if the dataset has 10 samples and a one-sample t-test is
% used, z-scores for samples with no NaN values is based on the
% t-statistic with df=9, but those with two missing values (NaN values)
% are based on the t-statistic with df=7. A use case is computing fMRI
% group statistics where overlap across brains is not perfect at the
% voxel-by-voxel level, in as similar approach as AFNI's 3dttest++ with
% the '-toz' option.
% - This function computes feature-wise statistics that are not corrected
% for multiple comparisons. To correct for multiple comparisons, see
% cosmo_montecarlo_cluster_stat.
%
% See also: anova1, ttest, ttest2, cosmo_montecarlo_cluster_stat
%
% # For CoSMoMVPA's copyright information and license terms, #
% # see the COPYING file distributed with CoSMoMVPA. #
if nargin < 3
output_stat_name = '';
end
[output_stat_name, tail] = get_stat_definition(stat_name, ...
output_stat_name);
[samples, targets, chunks, type] = get_descriptors(ds);
% run specified helper function
if isfield(ds.sa, 'contrast')
contrast = ds.sa.contrast;
else
contrast = [];
end
% ensure stat_name matches the number of classes
verify_class_count(stat_name, targets);
% get ttest1, ttest2 or Ftest function handle
stat_func = get_stat_func(stat_name);
% apply statistic function handle
[stat, df_struct, stat_label] = apply_stat_func(stat_func, ...
samples, targets, chunks, type, contrast);
% convert to p-value or z-score, if necessary
[stat, stat_label] = compute_output_stat(stat, df_struct, stat_label, ...
tail, output_stat_name);
% store output
stat_ds = struct();
if isfield(ds, 'a')
stat_ds.a = ds.a;
end
if isfield(ds, 'fa')
stat_ds.fa = ds.fa;
end
stat_ds.samples = stat;
stat_ds.sa.stats = {stat_label};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% helper functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [stat, df_struct, stat_label] = apply_stat_func(stat_func, ...
samples, targets, chunks, ...
type, contrast)
% applies the stat_func to samples, taking into account missing values
% (indicated by NaN) in samples
samples_nan_msk = isnan(samples);
has_missing = any(any(samples_nan_msk, 1) & ~all(samples_nan_msk, 1));
if ~has_missing
[stat, df_struct, stat_label] = apply_stat_func_no_nans(stat_func, ...
samples, targets, chunks, ...
type, contrast);
return
end
[nsamples, nfeatures] = size(samples);
stat = zeros(1, nfeatures);
zeroed_samples = samples;
zeroed_samples(~samples_nan_msk) = 0;
nan_counts = sum(samples_nan_msk, 1);
unq_counts = unique(nan_counts);
% build indicator matrix, so that a value in zeroed_samples(i,j)=v
% - v=0 means that samples(i,j) is not NaN
% - v=k, k>0 means that samples(i,j) is the k-th NaN value in that
% column
for c = 1:numel(unq_counts)
count = unq_counts(c);
col_msk = nan_counts == count;
msk = bsxfun(@and, samples_nan_msk, col_msk);
zeroed_samples(msk) = repmat((1:count)', 1, sum(col_msk));
end
% using the indicator matrix, split up the samples in blocks so that
% in each block, the same rows in samples are NaN. Each block is
% processed separately in the 'for' loop below.
[indices, zeroed_unq] = cosmo_index_unique(zeroed_samples');
nclasses = max(targets);
for k = 1:numel(indices)
zeroed = zeroed_unq(k, :)';
cols = indices{k};
nan_msk = zeroed > 0;
% keep samples where there are non-NaN values for the corresponding
% chunk
keep_msk = any(bsxfun(@eq, chunks(:), chunks(~nan_msk)'), 2);
keep_samples = samples(keep_msk, cols);
% ensure chunks are numbered 1:max(chunks)
[unused, unused, keep_chunks] = unique(chunks(keep_msk));
% select targets
keep_targets = targets(keep_msk);
% compute t or F statistic
[stat_k, df, stat_label_k] = stat_func(keep_samples, ...
keep_targets, ...
keep_chunks, ...
type, contrast);
% if for some samples in a chunk the value was NaN and for others
% it was not NaN, then the output is NaN
has_mixed_nan_and_non_nan = any(nan_msk(keep_msk));
% if a target value was completely non-present, then
% the output is NaN
keep_targets_msk = false(1, nclasses);
keep_targets_msk(targets(keep_msk)) = true;
if has_mixed_nan_and_non_nan || any(~keep_targets_msk)
stat_k(:) = NaN;
end
% store stat val
stat(cols) = stat_k;
if k == 1
% allocate space
df_count = numel(df);
df_matrix = zeros(df_count, nfeatures);
stat_label = stat_label_k;
else
% stat label must always be the same
assert(isequal(stat_label, stat_label_k));
end
ncols = numel(cols);
df_matrix(:, cols) = repmat(df(:), 1, ncols);
end
single_sample = samples(:, 1);
single_sample(:) = 0;
[unused, max_df] = stat_func(single_sample, targets, chunks, ...
type, contrast);
df_struct = struct();
df_struct.max_df = max_df(:);
df_struct.feature_wise_df = df_matrix;
function [stat, df_struct, stat_label] = apply_stat_func_no_nans(stat_func, ...
samples, targets, chunks, ...
type, contrast)
% much faster computation of statistics, if there are no NaNs
[stat, max_df, stat_label] = stat_func(samples, targets, chunks, ...
type, contrast);
nfeatures = size(samples, 2);
df_struct = struct();
df_struct.max_df = max_df(:);
df_struct.feature_wise_df = repmat(max_df(:), 1, nfeatures);
function f = get_stat_func(stat_name)
stat_name2func = struct();
stat_name2func.t = @ttest1_wrapper;
stat_name2func.t2 = @ttest2_wrapper;
stat_name2func.F = @ftest_wrapper;
assert(isfield(stat_name2func, stat_name));
f = stat_name2func.(stat_name);
function verify_class_count(stat_name, targets)
stat_name2interval = struct();
stat_name2interval.t = [1, 2];
stat_name2interval.t2 = [2, 2];
stat_name2interval.F = [2, Inf];
if ~isfield(stat_name2interval, stat_name)
error('illegal statname ''%s'', supported are: %s', stat_name, ...
cosmo_strjoin(fieldnames(stat_name2interval), ', '));
end
valid_interval = stat_name2interval.(stat_name);
count = max(targets);
assert(count == numel(unique(targets)));
min_count = valid_interval(1);
if count < min_count
error('statname ''%s'' requires at least %d unique targets', ...
stat_name, min_count);
end
max_count = valid_interval(2);
if count > max_count
error('statname ''%s'' requires at most %d unique targets', ...
stat_name, max_count);
end
function [stat, stat_label] = compute_output_stat(stat, df_struct, stat_name, ...
tail, output_stat_name)
% transform output is required
if isempty(output_stat_name)
% raw t or F value
stat_fa_name = stat_name;
% features with missing values (as indicated by df values that
% are different from the maximum df possible) are set to NaN
df_is_max_msk = bsxfun(@eq, df_struct.max_df, ...
df_struct.feature_wise_df);
has_missing = ~all(df_is_max_msk, 1);
stat(has_missing) = NaN;
% set degrees of freedom
df_str = arrayfun(@(x) sprintf('%d', x), df_struct.max_df, ...
'UniformOutput', false);
stat_label = sprintf('%s(%s)', stat_fa_name, ...
cosmo_strjoin(df_str, ','));
else
% transform to left-tailed p-value for each unique combination
% of degrees of freedom
stat = apply_cdf_wrapper_different_dfs(stat_name, stat, df_struct);
switch output_stat_name
case 'z'
% transform to z-score
stat = cosmo_norminv(stat);
stat_fa_name = 'Zscore';
case 'p'
switch tail
case 'left'
% do nothing
case 'right'
% invert p-value
stat = 1 - stat;
case 'both'
% take whichever tail is more extreme
stat = (.5 - abs(stat - .5)) * 2;
otherwise
assert(false, 'this should not happen');
end
stat_fa_name = 'Pval';
otherwise
error('illegal output type %s', output_stat_name);
end
stat_label = sprintf('%s()', stat_fa_name);
end
function cdf_stat = apply_cdf_wrapper_different_dfs(stat_name, stat, df_struct)
% applies the cdf wrapper to each unique combination of degrees of
% freedom stored in df_struct.feature_wise_df
[idx, unique_dfs] = cosmo_index_unique(df_struct.feature_wise_df');
% allocate space for output
cdf_stat = zeros(size(stat)) + NaN;
% compute result
for k = 1:numel(idx)
cols = idx{k};
df_cell = num2cell(unique_dfs(k, :));
col_stat = stat(cols);
msk = ~isnan(col_stat);
ps = cdf_wrapper(stat_name, col_stat(msk), df_cell{:});
% In Octave, extreme statistics can lead to negative p values
ps(ps < 0) = 0;
ps(ps > 1) = 1;
cdf_stat(cols(msk)) = ps;
end
function y = cdf_wrapper(name, x, df1, df2)
ensure_has_stats_functions();
if ~(df1 > 0)
y = zeros(size(x)) + NaN;
return
end
switch name
case 'Ttest'
assert(nargin == 3);
y = tcdf(x, df1);
case 'Ftest'
assert(nargin == 4);
if ~(df2 > 0)
y = zeros(size(x)) + NaN;
return
end
y = fcdf(x, df1, df2);
otherwise
assert(false);
end
function [stat, df, stat_label] = ttest1_wrapper(samples, targets, chunks, ...
type, contrast)
if ~isempty(contrast)
error('contrast is not supported for t-stat');
end
nclasses = max(targets);
if nclasses == 2
samples = pairwise_differences(samples, targets, chunks);
nclasses = 1;
end
if nclasses ~= 1
error('t-stat: expected 1 or 2 classes, found %d', ...
nclasses);
end
[stat, df] = quick_ttest(samples);
stat_label = 'Ttest';
function [stat, df, stat_label] = ttest2_wrapper(samples, targets, chunks, ...
type, contrast)
if ~isempty(contrast)
error('contrast is not supported for t-stat');
end
if strcmp(type, 'within')
error(['ttest2 stat: the values in chunks and targets suggest '...
'a within-subject design. If you want to '...
'run a paired-test, use ''t'' (not ''t2'') ', ...
'as the second argument']);
end
m1 = targets == 1;
m2 = targets == 2;
[stat, df] = quick_ttest2(samples(m1, :), ...
samples(m2, :));
nclasses = max(targets);
if nclasses ~= 2
% missing targets, set all to NaN
assert(nclasses < 2);
stat(:) = NaN;
end
stat_label = 'Ttest';
function [stat, df, stat_label] = ftest_wrapper(samples, targets, chunks, ...
type, contrast)
nclasses = max(targets);
if nclasses < 2
error('F stat: expected >=2 classes, found %d', nclasses);
end
switch type
case 'between'
[stat, df] = quick_ftest_between(samples, targets, ...
chunks, contrast);
case 'within'
[stat, df] = quick_ftest_within(samples, targets, ...
chunks, contrast);
end
stat_label = 'Ftest';
function [f, df] = quick_ftest_between(samples, targets, chunks, contrast)
% one-way ANOVA
has_contrast = ~isempty(contrast);
contrast_sum = 0;
nclasses = max(targets);
[ns, nf] = size(samples);
mu = sum(samples, 1) / ns; % grand mean
b = zeros(nclasses, nf); % between-class sum of squares
nsc = zeros(nclasses, 1);
wss = 0; % within-class sum of squares
for k = 1:nclasses
msk = k == targets;
nsc(k) = sum(msk); % number of samples in this class
sample = samples(msk, :);
muc = sum(sample, 1) / nsc(k); % class mean
% between- and within-class sum of squares
if has_contrast
cmsk = contrast(msk);
if ~all(cmsk(1) == cmsk)
error('Contrast has differerent values in level %d', k);
end
contrast_sum = contrast_sum + cmsk(1);
b(k, :) = sum(bsxfun(@times, contrast(msk), mu - muc), 1);
else
b(k, :) = (mu - muc);
end
wss = wss + sum(bsxfun(@minus, muc, sample).^2, 1);
end
if has_contrast
if contrast_sum ~= 0
error('contrast has sum %d, should be 0', contrast_sum);
end
bss = sum(b, 1).^2 / sum(contrast.^2);
df1 = 1;
else
bss = sum(bsxfun(@times, nsc, b.^2), 1);
df1 = nclasses - 1;
end
df = [df1, ns - nclasses];
mbss = bss / df(1);
mwss = wss / df(2);
f = zeros(1, nf) + NaN;
msk = mbss > 0;
f(msk) = mbss(msk) ./ mwss(msk);
function [f, df] = quick_ftest_within(samples, targets, chunks, contrast)
% repeated measures anova
if ~isempty(contrast)
error('contrast is not supported for within-subject design');
end
nchunks = max(chunks);
nclasses = max(targets);
nfeatures = size(samples, 2);
gm = mean(samples, 1); % grand mean
sst = zeros(1, nfeatures);
ssw = zeros(1, nfeatures);
for k = 1:nclasses
xk = samples(k == targets, :);
n = size(xk, 1);
if n == 0
ssw(:) = NaN;
break
end
mu = sum(xk, 1) / n;
sst = sst + n * (gm - mu).^2;
ssw = ssw + sum(bsxfun(@minus, mu, xk).^2);
end
sss = zeros(1, nfeatures);
for k = 1:nchunks
xk = samples(k == chunks, :);
n = size(xk, 1);
mu = mean(xk, 1);
sss = sss + n * (gm - mu).^2;
end
df1 = (nclasses - 1);
mst = sst / df1;
df2 = df1 * (nchunks - 1);
sse = ssw - sss;
mse = sse / df2;
msk = mse > 0;
f = zeros(1, nfeatures) + NaN;
f(msk) = mst(msk) ./ mse(msk);
df = [df1 df2];
function [t, df] = quick_ttest(x)
% one-sample t-test against zero
[ns, nf] = size(x);
mu = sum(x, 1) / ns; % grand mean
df = ns - 1;
scaling = ns * df;
% sum of squares
ss = sum(bsxfun(@minus, x, mu).^2, 1);
t = zeros(1, nf) + NaN;
msk = ss > 0;
t(msk) = mu(msk) .* sqrt(scaling ./ ss(msk));
function [t, df] = quick_ttest2(x, y)
% two-sample t-test with equal variance assumption
[nx, nf] = size(x);
ny = size(y, 1);
df = nx + ny - 2;
if nx == 0 || ny == 0
t = zeros(1, nf) + NaN;
return
end
mux = sum(x, 1) / nx; % mean of class x
muy = sum(y, 1) / ny; % " " y
scaling = (nx * ny) * df / (nx + ny);
% sum of squares
ss = sum([bsxfun(@minus, x, mux); bsxfun(@minus, y, muy)].^2, 1);
t = zeros(1, nf) + NaN;
msk = ss > 0;
t(msk) = (mux(msk) - muy(msk)) .* sqrt(scaling ./ ss(msk));
function ensure_has_stats_functions()
% - Octave has the required functionality in the octave-forge
% statistics toolbox and will raise an error if it is not installed.
% - Matlab will raise an error message that the statistics toolbox is
% required
persistent cached_has_stat_funcs
if isequal(cached_has_stat_funcs, [])
if cosmo_wtf('is_matlab') ...
&& isempty(which('tinv')) ...
&& isempty(which('fcdf'))
raise_exception_if_absent = true;
cached_has_stat_funcs = cosmo_check_external('@stats', ...
raise_exception_if_absent);
else
cached_has_stat_funcs = false;
end
end
function [samples, targets, chunks, type] = get_descriptors(ds)
cosmo_isfield(ds, {'samples', 'sa.chunks', 'sa.targets'}, true);
samples = ds.samples;
% unique targets
[unused, unused, targets] = unique(ds.sa.targets);
% unique chunks
[unq_chunks, unused, chunks] = unique(ds.sa.chunks);
nc = max(chunks);
if isequal(sort(ds.sa.chunks), unq_chunks)
type = 'between';
else
combis = (targets - 1) * nc + chunks;
if isequal(sort(combis), (1:numel(combis))')
type = 'within';
else
error(['Either all chunks must be unique, or each chunk '...
'must contain the same targets']);
end
end
function delta = pairwise_differences(samples, targets, chunks)
% for one-sample t-test: compute differences between first and second
% target
n = size(samples, 1);
assert(numel(targets) == n);
assert(numel(chunks) == n);
assert(mod(n, 2) == 0);
n2 = n / 2;
idxs = zeros(n2, 2);
for k = 1:n
assert(idxs(chunks(k), targets(k)) == 0);
idxs(chunks(k), targets(k)) = k;
end
assert(isequal(size(idxs), [n2, 2]));
assert(all(idxs(:) > 0));
delta = samples(idxs(:, 1), :) - samples(idxs(:, 2), :);
function [output_stat_name, tail] = get_stat_definition(stat_name, ...
output_stat_name)
if any(cosmo_match({'left', 'right', 'both'}, output_stat_name))
tail = output_stat_name;
output_stat_name = 'p';
elseif strcmp(output_stat_name, 'p')
switch stat_name
case 'F'
tail = 'right'; % show anova1 behaviour w.r.t. p-values
otherwise
tail = 'both'; % show ttest[2] " "
end
else
tail = 'both';
end
