function ds_sa = cosmo_correlation_measure(ds, varargin)
% Computes a split-half correlation measure
%
% d=cosmo_correlation_measure(ds[, args])
%
% Inputs:
% ds dataset structure with fields .samples, .sa.targets and
% .sa.chunks
% args optional struct with the following optional fields:
% .partitions struct with fields .train_indices and .test_indices.
% Both should be a Px1 cell for P partitions. If omitted,
% it is set to cosmo_nchoosek_partitioner(ds,'half').
% .template QxQ matrix for Q classes in each chunk. This matrix
% weights the correlations across the two halves.
% If ds.sa.targets has only one unique value, it must be
% set to the scalar value 1; otherwise it should
% have a mean of zero. If omitted, it has positive values
% of (1/Q) on the diagonal and (-1/(Q*(Q-1)) off the
% diagonal.
% (Note: this can be used to test for representational
% similarity matching)
% .merge_func A function handle used to merge data from matching
% targets in the same chunk. Default is @(x) mean(x,1),
% meaning that values are averaged over the same samples.
% It is assumed that isequal(args.merge_func(y),y) if y
% is a row vector.
% .corr_type Type of correlation: 'Pearson','Spearman','Kendall'.
% The default is 'Pearson'.
% .post_corr_func Operation performed after correlation. (default:
% @atanh)
% .output 'mean' (default): correlations weighted by template
% 'raw' or 'correlation': correlations between all classes
% 'one_minus_correlation': 1 minus correlations
% 'mean_by_fold': provide weighted correlations for each
% fold in the partitions.
%
%
% Output:
% ds_sa Struct with fields:
% .samples Scalar indicating how well the template matrix
% correlates with the correlation matrix from the two
% halves (averaged over partitions). By default:
% - this value is based on Fisher-transformed correlation
% values, not raw correlation values
% - this is the average of the (Fisher-transformed)
% on-diagonal minus the average of the
% (Fisher-transformed) off-diagonal elements of the
% correlation matrix based on the two halves of the data.
% .sa Struct with field:
% .labels if output=='corr'
% .half1 } if output=='raw': (N^2)x1 vectors with indices of data
% .half2 } from two halves, with N the number of unique targets.
%
% Example:
% ds=cosmo_synthetic_dataset();
% %
% % compute on-minus-off diagonal correlations
% c=cosmo_correlation_measure(ds);
% cosmo_disp(c)
% %|| .samples
% %|| [ 1.23 ]
% %|| .sa
% %|| .labels
% %|| { 'corr' }
% %
% % Spearman correlation requires the Matlab statistics toolbox
% cosmo_skip_test_if_no_external('@stats');
% %
% c=cosmo_correlation_measure(ds,'corr_type','Spearman');
% cosmo_disp(c)
% %|| .samples
% %|| [ 1.28 ]
% %|| .sa
% %|| .labels
% %|| { 'corr' }
%
% ds=cosmo_synthetic_dataset();
% % get raw correlations without fisher transformation
% c_raw=cosmo_correlation_measure(ds,'output','correlation',...
% 'post_corr_func',[]);
% cosmo_disp(c_raw)
% %|| .samples
% %|| [ 0.363
% %|| -0.404
% %|| -0.447
% %|| 0.606 ]
% %|| .sa
% %|| .half1
% %|| [ 1
% %|| 2
% %|| 1
% %|| 2 ]
% %|| .half2
% %|| [ 1
% %|| 1
% %|| 2
% %|| 2 ]
% %|| .a
% %|| .sdim
% %|| .labels
% %|| { 'half1' 'half2' }
% %|| .values
% %|| { [ 1 [ 1
% %|| 2 ] 2 ] }
% %
% % convert to matrix form (N x N x P, with N the number of classes and
% % P=1)
% matrices=cosmo_unflatten(c_raw,1);
% cosmo_disp(matrices)
% %|| [ 0.363 -0.447
% %|| -0.404 0.606 ]
%
% % compute for each fold separately, using a custom take-one-chunk
% % out partitioning scheme and without Fisher transformation of
% % the correlations.
% % Note that c.sa.chunks in the output reflects the test chunk
% % in each partition
% ds=cosmo_synthetic_dataset('type','fmri','nchunks',4);
% partitions=cosmo_nfold_partitioner(ds);
% c=cosmo_correlation_measure(ds,'output','mean_by_fold',...
% 'partitions',partitions,...
% 'post_corr_func',[]);
% cosmo_disp(c.samples);
% %|| [ 1.32
% %|| 0.512
% %|| 1.05
% %|| 1.23 ]
% cosmo_disp(c.sa);
% %|| .partition
% %|| [ 1
% %|| 2
% %|| 3
% %|| 4 ]
%
% % minimal searchlight example
% ds=cosmo_synthetic_dataset('type','fmri');
% % use searchlight with radius 1 voxel (radius=3 is more typical)
% nbrhood=cosmo_spherical_neighborhood(ds,'radius',1,'progress',false);
% % run searchlight
% res=cosmo_searchlight(ds,nbrhood,@cosmo_correlation_measure,...
% 'progress',false);
% cosmo_disp(res.samples)
% %|| [ 1.87 1.25 1.51 1.68 1.71 0.879 ]
% cosmo_disp(res.sa)
% %|| .labels
% %|| { 'corr' }
%
% Notes:
% - by default the post_corr_func is set to @atanh. This is equivalent to
% a Fisher transformation from r (correlation) to z (standard z-score).
% The underlying math is z=atanh(r)=.5*log((1+r)./log(1-r)).
% The rationale is to make data more normally distributed under the
% null hypothesis.
% Fisher-transformed correlations can be transformed back to
% their original correlation values using 'tanh', which is the inverse
% of 'atanh'.
% - To disable the (by default used) Fisher-transformation, set the
% 'post_corr_func' option to [].
% - if multiple samples are present with the same chunk and target, they
% are averaged *prior* to computing the correlations.
% - if multiple partitions are present, then the correlations are
% computed separately for each partition, and then averaged (unless
% the 'output' option is set, and set to a different value than
% 'mean'.
% - When more than two chunks are present in the input, partitions
% consist of all possible half splits for which the number of unique
% chunks in the train and test set differ by 1 at most.
% For illustration, up to 6 chunks, the
% partitions are:
% - 2 chunks - partition #1
% chunks first half: { 1
% chunks second half: { 2
%
% - 3 chunks - partition #1 #2 #3
% chunks first half: { 3 2 1
% { 1 1 2
% chunks second half: { 2 3 3
%
% - 4 chunks - partition #1 #2 #3
% { 3 2 2
% chunks first half: { 4 4 3
% { 1 1 1
% chunks second half: { 2 3 4
%
% - 5 chunks - partition #1 #2 #3 #4 #5 #6 #7 #8 #9 #10
% { 4 3 3 2 2 2 1 1 1 1
% chunks first half: { 5 5 4 5 4 3 5 4 3 2
% { 1 1 1 1 1 1 2 2 2 3
% chunks second half: { 2 2 2 3 3 4 3 3 4 4
% { 3 4 5 4 5 5 4 5 5 5
%
% - 6 chunks - partition #1 #2 #3 #4 #5 #6 #7 #8 #9 #10
% { 4 3 3 3 2 2 2 2 2 2
% chunks first half: { 5 5 4 4 5 4 4 3 3 3
% { 6 6 6 5 6 6 5 6 5 4
% { 1 1 1 1 1 1 1 1 1 1
% chunks second half: { 2 2 2 2 3 3 3 4 4 5
% { 3 4 5 6 4 5 6 5 6 6
% Thus, with an increasing number of chunks, the number of partitions
% (and thus the time required to run this function) increases
% quadratically. To use simpler partition schemes (e.g. odd-even, as
% provided by cosmo_oddeven_partitioner), specify the 'partitions'
% argument.
%
% References
% - Haxby, J. V. et al (2001). Distributed and overlapping
% representations of faces and objects in ventral temporal cortex.
% Science 293, 2425?2430
%
% # For CoSMoMVPA's copyright information and license terms, #
% # see the COPYING file distributed with CoSMoMVPA. #
persistent cached_partitions
persistent cached_chunks
persistent cached_params
persistent cached_varargin
% optimize parameters parsing: if same arguments were used in
% a previous call, do not use cosmo_structjoin (which is relative
% expensive)
if ~isempty(cached_params) && isequal(cached_varargin, varargin)
params = cached_params;
else
defaults = struct();
defaults.partitions = [];
defaults.template = [];
defaults.merge_func = [];
defaults.corr_type = 'Pearson';
defaults.post_corr_func = @atanh;
defaults.output = 'mean';
defaults.check_partitions = true;
params = cosmo_structjoin(defaults, varargin);
show_warning_if_no_defaults(defaults, params);
cached_params = params;
cached_varargin = varargin;
end
partitions = params.partitions;
template = params.template;
merge_func = params.merge_func;
post_corr_func = params.post_corr_func;
check_partitions = params.check_partitions;
chunks = ds.sa.chunks;
if isempty(partitions)
if ~isempty(cached_chunks) && isequal(cached_chunks, chunks)
partitions = cached_partitions;
% assume that partitions were already checked
check_partitions = false;
else
partitions = cosmo_nchoosek_partitioner(ds, 'half');
cached_chunks = chunks;
cached_partitions = partitions;
end
end
if check_partitions
cosmo_check_partitions(partitions, ds, params);
end
targets = ds.sa.targets;
nsamples = size(targets, 1);
[classes, unused, class_ids] = fast_unique(targets);
nclasses = numel(classes);
switch nclasses
case 0
error('No classes found - this is not supported');
case 1
if ~isequal(template, 1)
error(['Only one unique value for .sa.targets was found; '...
'this is only allowed when the ''template'' '...
'parameter is explicitly set to 1.\n'...
'Note that this option does not compute the '...
'correlation *difference* between matching and '...
'non-matching values in .sa.targets; instead it '...
'computes the direct correlation '...
'between two halves of the data. A typical use '...
'case is when .samples already contains a '...
'difference score comparing two different '...
'conditions']);
end
otherwise
if isempty(template)
template = (eye(nclasses) - 1 / nclasses) / (nclasses - 1);
else
max_tolerance = 1e-8;
if abs(sum(template(:))) > max_tolerance
error('Template matrix does not have a sum of zero');
end
expected_size = [nclasses, nclasses];
if ~isequal(size(template), [nclasses, nclasses])
error('Template must be of size %dx%d', expected_size);
end
end
end
template_msk = isfinite(template);
npartitions = numel(partitions.train_indices);
% space for output
pdata = cell(npartitions, 1);
% keep track of how often each chunk was used in test_indices,
% and which chunk was used last
test_chunks_last = NaN(nsamples, 1);
test_chunks_count = zeros(nsamples, 1);
for k = 1:npartitions
train_indices = partitions.train_indices{k};
test_indices = partitions.test_indices{k};
% get data in each half
half1 = get_data(ds, train_indices, class_ids, merge_func);
half2 = get_data(ds, test_indices, class_ids, merge_func);
% compute raw correlations
raw_c = cosmo_corr(half1', half2', params.corr_type);
% apply post-processing (usually Fisher-transform, i.e. atanh)
c = apply_post_corr_func(post_corr_func, raw_c);
% aggregate results
pdata{k} = aggregate_correlations(c, template, template_msk, params.output);
test_chunks_last(test_indices) = chunks(test_indices);
test_chunks_count(test_indices) = test_chunks_count(test_indices) + 1;
end
ds_sa = struct();
switch params.output
case 'mean'
ds_sa.samples = mean(cat(2, pdata{:}), 2);
ds_sa.sa.labels = {'corr'};
case {'raw', 'correlation'}
ds_sa.samples = mean(cat(2, pdata{:}), 2);
nclasses = numel(classes);
ds_sa.sa.half1 = reshape(repmat((1:nclasses)', nclasses, 1), [], 1);
ds_sa.sa.half2 = reshape(repmat(1:nclasses, nclasses, 1), [], 1);
ds_sa.a.sdim = struct();
ds_sa.a.sdim.labels = {'half1', 'half2'};
ds_sa.a.sdim.values = {classes, classes};
case 'mean_by_fold'
ds_sa.sa.partition = (1:npartitions)';
ds_sa.samples = [pdata{:}]';
otherwise
assert(false, 'this should be caught by get_data');
end
function c = apply_post_corr_func(post_corr_func, c)
if ~isempty(post_corr_func)
c = post_corr_func(c);
end
function agg_c = aggregate_correlations(c, template, template_msk, output)
switch output
case {'mean', 'mean_by_fold'}
pcw = c(template_msk) .* template(template_msk);
agg_c = sum(pcw(:));
case {'raw', 'correlation'}
agg_c = c(:);
case 'one_minus_correlation'
error(['the ''output'' option ''one_minus_correlation'' '...
'has been removed. Please contact CoSMoMVPA''s '...
'authors if you really need this option']);
otherwise
error('Unsupported output %s', output);
end
function [unq, pos, ids] = fast_unique(x)
% optimized for vectors
[y, i] = sort(x, 1);
msk = [true; diff(y) > 0];
j = find(msk);
pos = i(j);
unq = y(j);
vs = cumsum(msk);
ids = vs;
ids(i) = vs;
function data = get_data(ds, sample_idxs, class_ids, merge_func)
samples = ds.samples(sample_idxs, :);
target_ids = class_ids(sample_idxs);
nclasses = max(class_ids);
merge_by_averaging = isempty(merge_func);
if isequal(target_ids', 1:nclasses) && merge_by_averaging
% optimize standard case of one sample per class and normal
% averaging over samples
data = samples;
return
end
nfeatures = size(samples, 2);
data = zeros(nclasses, nfeatures);
for k = 1:nclasses
msk = target_ids == k;
n = sum(msk);
class_samples = samples(msk, :);
if merge_by_averaging
if n == 1
data(k, :) = class_samples;
else
data(k, :) = sum(class_samples, 1) / n;
end
else
data(k, :) = merge_func(class_samples);
end
if n == 0
error('missing target class %d', class_ids(k));
end
end
function show_warning_if_no_defaults(defaults, params)
keys = {'output', 'post_corr_func'};
has_defaults = isequal(select_fields(defaults, keys), ...
select_fields(params, keys));
if ~has_defaults && isequal(params.post_corr_func, @atanh)
name = mfilename();
msg = sprintf( ...
['Please note that the ''%s'' function applies '...
'Fisher transformation after the correlations '...
'have been computed. This was a somewhat unfortunate '...
'implementation decision that will not be changed '...
'to avoid breaking behaviour with earlier versions.\n'...
'To disable using the Fisher transformation, use the '...
'''%s'' function while setting the ''post_corr_func'''...
'option to the empty array ([])'], name, name);
cosmo_warning(msg);
end
function subset = select_fields(s, keys)
n = numel(keys);
for k = 1:n
key = keys{k};
subset.(key) = s.(key);
end
