diff --git a/+rsa/+stat/FDRthreshold.m b/+rsa/+stat/FDRthreshold.m index 82ff06c..35445fc 100644 --- a/+rsa/+stat/FDRthreshold.m +++ b/+rsa/+stat/FDRthreshold.m @@ -54,13 +54,13 @@ % statistical maps in functional neuroimaging using the false % discovery rate. NeuroImage 15, 870-878. 2002. -import rsa.* -import rsa.fig.* -import rsa.fmri.* -import rsa.rdm.* -import rsa.sim.* -import rsa.spm.* -import rsa.stat.* +import rsa.* +import rsa.fig.* +import rsa.fmri.* +import rsa.rdm.* +import rsa.sim.* +import rsa.spm.* +import rsa.stat.* import rsa.util.* %% PREPARATIONS @@ -110,4 +110,8 @@ disp([num2str(nSig),' tests exceed the FDR threshold for q<',num2str(q),'.']); +if isempty(pThreshold) + pThreshold = inf; +end + end%function diff --git a/+rsa/+stat/bootConfint.m b/+rsa/+stat/bootConfint.m new file mode 100644 index 0000000..c354bdf --- /dev/null +++ b/+rsa/+stat/bootConfint.m @@ -0,0 +1,173 @@ +function [ci_lo, ci_up, p] = bootConfint(x, B, alternative, user_opts) + +% bootConfint +% Compute confidence intervals based on bootstrap distribution using +% various methods and use them to derive p-values for hypothesis testing. +% +% [ci_lo, ci_up, p] = bootConfint(x, B, alternative, user_opts) +% +% Currently two methods are available, "percentile" and "bca". P-values +% can be computed for both one-tailed and two-tailed testing. +% +% INPUT +% ----- +% x : vector, original sample of length #subjects +% B : vector, bootstrap distribution of the statistic of +% interest of length #bootstrap_iterations +% alternative : string, either 'two-tailed' or 'greater' +% user_opts : struct, user options used by RSA Toolbox. Options are ... +% user_opts.bootCImethod: +% - "percentile": bootstrap percentile method (default) +% - "bca": bias-corrected and accelerated CI +% user_opts.bootCIalpha: type 1 error probability to use +% for CI definition, default is .05. +% +% OUTPUT +% ------ +% ci_lo, ci_up: lower and upper bounds of CI, respectively +% p : p value of hypothesis test +% +% Literature +% ---------- +% "An Introduction to the Bootstrap" (Efron & Tibshirani, 1993) +% +% 2018-01-27 michael.bannert@tuebingen.mpg.de + +user_opts = rsa.util.setIfUnset(user_opts, 'bootCImethod', 'percentile'); +user_opts = rsa.util.setIfUnset(user_opts, 'bootCIalpha', .05); + +method = user_opts.bootCImethod; +ci_alpha = user_opts.bootCIalpha; + +n_B = length(B); +n_samples = length(x); +B_sorted = sort(B); + +switch lower(method) + case 'bca' + % Bias-correction, eq. 14.14 + z_hat0 = norminv(sum(B < mean(x)) / n_B); + + % Compute acceleration using jackknife + theta_i = nan(n_samples, 1); + for n = 1:n_samples + x_i = x(setdiff(1:n_samples, n)); + theta_i(n) = mean(x_i); + end + theta_pt = mean(theta_i); + + % Eq. 14.15 + a_hat = sum((theta_pt * ones(n_samples, 1) - theta_i).^3) / ... + 6 * (sum((theta_pt * ones(n_samples, 1) - theta_i).^2)) ^(3/2); + + z_alpha = norminv(ci_alpha / 2); + z_one_minus_alpha = norminv(1 - ci_alpha / 2); + + % Eq. 14.10 + alpha1 = normcdf( ... + z_hat0 + (z_hat0 + z_alpha) / ... + (1 - a_hat * (z_hat0 + z_alpha)) ... + ); + + alpha2 = normcdf( ... + z_hat0 + (z_hat0 + z_one_minus_alpha) / ... + (1 - a_hat * (z_hat0 + z_one_minus_alpha)) ... + ); + + % Find indices of bootstrap distribution corresponding to desired + % alpha + alpha1_n = max(floor(alpha1 * n_B), 1); % At least 1 + alpha2_n = ceil(alpha2 * n_B); + + % Eq. 14.9 + ci_lo = B_sorted(alpha1_n); + ci_up = B_sorted(alpha2_n); + + % Calculate p-value by solving eq. 14.10 for alpha. P-value should + % be one-tailed for the null hypothesis that the true value is <= 0. + % Find the z_alpha for which the threshold 0 equals the lower CI + % bound. + + % 1. Precompute inverse cumulative Gaussian of lower confidence + % bound that would just include 0. (Used several times in equation + % below.) + + % Index must be at least 1. Increase number of bootstrap iterations + % for more precise p-values + alpha0_n = find(B_sorted > 0, 1) - 1; + alpha1_norminv = norminv(alpha0_n / n_B); + + % 2. Find CI alpha for which the lower bound would just enclose + % 0. + % This is eq. 14.10 solved for z_alpha + z_alpha0 = ... + (alpha1_norminv * (1 - a_hat * z_hat0) + ... + z_hat0 * (z_hat0 * a_hat - 2)) / ... + (a_hat * (alpha1_norminv - z_hat0) - 1); + + % Super non-significant, i.e., no positive values in bootstrap + % distribution + if isempty(z_alpha0) + z_alpha0 = -inf; + end + + % Super significant, i.e., only positive values in bootstrap + % distribution. P-value only limited by number of bootstrap + % iterations. Increasing number of iterations would possibly help. + if isnan(z_alpha0) + z_alpha0 = norminv(1 - 1/n_B); + end + + switch lower(alternative) + case 'greater' + % Convert z_alpha to p-value + p = 1 - normcdf(z_alpha0); + + case 'two-tailed' + % For which z_alpha does the upper CI bound equal 0? See + % one-tailed case for explanation + + % Index must be at least n_B - 1. Increase number of + % bootstraps for more precise p-values + alpha0_n_up = min(n_B-1, find(B_sorted < 0, 1, 'last') + 1); + if isempty(alpha0_n_up) + alpha0_n_up = n_B-1; + end + alpha1_norminv_up = norminv(alpha0_n_up / n_B); + + z_alpha0_up = ... + (alpha1_norminv_up * (1 - a_hat * z_hat0) + ... + z_hat0 * (z_hat0 * a_hat - 2)) / ... + (a_hat * (alpha1_norminv_up - z_hat0) - 1); + + % Compare the two one-tailed p-values and return the + % smaller one multiplied by two. + p = min(1 - normcdf(z_alpha0), normcdf(z_alpha0_up)) * 2; + if isempty(p) || p > 1 + keyboard; + end + end + + case 'percentile' + % Efron & Tibshirani (1993, p. 170) + ci_lo = B_sorted(max(1, floor((ci_alpha / 2) * n_B))); + ci_up = B_sorted(ceil((1 - ci_alpha / 2) * n_B)); + + switch lower(alternative) + case 'greater' + % One-tailed p-value for null hypothesis that true value is + % below or equal to 0. + p = max(sum(B < 0) / n_B, 1/n_B); + + case 'two-tailed' + % Two-tailed p-value for null hypothesis that true value + % lies within confidence interval + p = max((min(sum(B < 0), sum(B >= 0)) / n_B) * 2, 1/n_B); + end + + otherwise + error(['Valid method ''%s'' to compute bootstrap CIs\nValid ' ... + 'methods are: ''bca'', ''percentile'''], method); +end + +assert(p <= 1 && p >= 0, '''p'' is not a probability.'); diff --git a/+rsa/+stat/bootstrapRDMs.m b/+rsa/+stat/bootstrapRDMs.m index 0e99fa9..bc850cf 100644 --- a/+rsa/+stat/bootstrapRDMs.m +++ b/+rsa/+stat/bootstrapRDMs.m @@ -1,4 +1,4 @@ -function [realRs bootstrapEs pairwisePs bootstrapRs] = bootstrapRDMs(bootstrappableReferenceRDMs, candRDMs, userOptions) +function [realRs, bootstrapEs, pairwisePs, bootstrapRs] = bootstrapRDMs(bootstrappableReferenceRDMs, candRDMs, userOptions) % [realRs bootstrapEs pairwisePs bootstrapRs] = ... % bootstrapRDMs(bootstrappableReferenceRDMs, ... % candRDMs, ... @@ -38,13 +38,13 @@ %__________________________________________________________________________ % Copyright (C) 2010 Medical Research Council -import rsa.* -import rsa.fig.* -import rsa.fmri.* -import rsa.rdm.* -import rsa.sim.* -import rsa.spm.* -import rsa.stat.* +import rsa.* +import rsa.fig.* +import rsa.fmri.* +import rsa.rdm.* +import rsa.sim.* +import rsa.spm.* +import rsa.stat.* import rsa.util.* % Sort out defaults @@ -115,31 +115,38 @@ n = 0; k=0; fprintf('\n'); + +% Need to create one candidate RDM for each reference because we don't want +% to correlate averaged RDMs. +candRDMs3rd = cell(nCandRDMs, 1); +for candRDMI = 1:nCandRDMs + candRDMs3rd{candRDMI} = repmat(candRDMs(:,:,candRDMI), ... + [1, 1, nSubjects]); +end + for candRDMI = 1:nCandRDMs for b = 1:userOptions.nBootstrap n = n + 1; localReferenceRDMs = bootstrappableReferenceRDMs(resampledConditionIs(b,:),resampledConditionIs(b,:),resampledSubjectIs(b,:)); - localTestRDM = candRDMs(resampledConditionIs(b,:), resampledConditionIs(b,:), candRDMI); - - averageBootstrappedRDM = mean(localReferenceRDMs, 3); - + localTestRDM = candRDMs3rd{candRDMI}(resampledConditionIs(b,:), resampledConditionIs(b,:), resampledSubjectIs(b,:)); + if isequal(userOptions.RDMcorrelationType,'Kendall_taua') % tic - bootstrapRs(candRDMI, b)=rankCorr_Kendall_taua(vectorizeRDMs(averageBootstrappedRDM)',vectorizeRDMs(localTestRDM)'); + bootstrapRs(candRDMI, b)=mean(diag(rankCorr_Kendall_taua(squeeze(vectorizeRDMs(localReferenceRDMs)),squeeze(vectorizeRDMs(localTestRDM))))); % toc elseif isequal(userOptions.RDMcorrelationType,'raeSpearman') % tic - % bootstrapRs(candRDMI, b)=raeSpearmanCorr(vectorizeRDMs(averageBootstrappedRDM)',vectorizeRDMs(localTestRDM)'); + % bootstrapRs(candRDMI, b)=raeSpearmanCorr(vectorizeRDMs(averageBootstrappedRDM),vectorizeRDMs(localTestRDM)); % toc % tic - bootstrapRs(candRDMI, b)=raeSpearmanCorr(vectorizeRDMs(averageBootstrappedRDM)',vectorizeRDMs(localTestRDM)'); + bootstrapRs(candRDMI, b)=mean(diag(raeSpearmanCorr(vectorizeRDMs(localReferenceRDMs),vectorizeRDMs(localTestRDM)))); % toc else % tic - bootstrapRs(candRDMI, b) = corr(vectorizeRDMs(averageBootstrappedRDM)',vectorizeRDMs(localTestRDM)','type',userOptions.distanceMeasure,'rows','pairwise'); + bootstrapRs(candRDMI, b) = mean(diag(corr(squeeze(vectorizeRDMs(localReferenceRDMs)), squeeze(vectorizeRDMs(localTestRDM)), 'type',userOptions.distanceMeasure,'rows','pairwise'))); % toc end - + if mod(n,floor(userOptions.nBootstrap*nCandRDMs/100))==0 fprintf('%d%% ',floor(100*n/(userOptions.nBootstrap*nCandRDMs))); if mod(n,floor(userOptions.nBootstrap*nCandRDMs/10))==0, fprintf('\n'); end; diff --git a/+rsa/compareRefRDM2candRDMs.m b/+rsa/compareRefRDM2candRDMs.m index 602555f..ace9773 100644 --- a/+rsa/compareRefRDM2candRDMs.m +++ b/+rsa/compareRefRDM2candRDMs.m @@ -20,8 +20,8 @@ % to assess, for each pair of candidate RDMs, which of the two RDMs is more % consistent with the reference RDM. A significant pairwise difference is % indicated by a horizontal line above the corresponding two bars. Each bar -% comes with an error bar, which indicates the standard error, estimated by -% the same procedure as is used for the pairwise candidate comparisons +% comes with an error bar, which indicates the standard error, estimated by +% the same procedure as is used for the pairwise candidate comparisons % (dependent on input data and userOptions, see below). % % Statistical inference on the correlation between the reference RDM and @@ -112,7 +112,7 @@ % correlation. (Note that multiple independent measurements of the % reference or candidate RDMs could also come from repeated % measurements within one subject. We refer to the instances as -% “subjects”, because subject random-effects inference is the most +% �subjects�, because subject random-effects inference is the most % common case.) % % 'randomisation': Test the relatedness of the reference RDM to each @@ -148,7 +148,7 @@ % subjects and conditions. This more conservative test attempts to % support inference generalising across both subjects and stimuli (to % their respective populations). Again, the usual caveats for basic -% bootstrap tests apply. +% bootstrap tests apply. % % 'none': Omit the test of RDM relatedness. % @@ -183,7 +183,7 @@ % this setting, userOptions.RDMrelatednessThreshold is interpreted as % the uncorrected p threshold. % -% userOptions.candRDMdifferencesTest +% userOptions.candRDMdifferencesTest % 'subjectRFXsignedRank' (default, data permitting): For each pair of % candidate RDMs, perform a statistical comparison to determine which % candidate RDM better explains the reference RDM by using the @@ -309,7 +309,7 @@ % pairwise candidate RDM comparisons. The first panel shows the % uncorrected-p matrix. The second panel shows the thresholded % uncorrected-p matrix. The third panel shows the FDR-thresholded p -% matrix. The fourth panel shows the Bonferroni-thresholded p matrix. +% matrix. The fourth panel shows the Bonferroni-thresholded p matrix. % % userOptions.resultsPath (default: pwd, i.e. current working directory) % This argument specifies the absolute path in which both figures are @@ -327,7 +327,7 @@ % If true (default: false), the figures are saved in Matlab .fig % format in userOptions.resultsPath. % -% userOptions.figure1filename +% userOptions.figure1filename % The filename for the bargraph figure, if chosen to be saved % (default: 'compareRefRDM2candRDMs_barGraph'). % @@ -341,11 +341,11 @@ % effect sizes and p values. % stats_p_r.candRelatedness_r: average correlations to reference RDM % stats_p_r.candRelatedness_p: corresponding uncorrected p values -% stats_p_r.SEs: standard errors of average RDM +% stats_p_r.SEs: standard errors of average RDM % correlations % stats_p_r.candDifferences_r: matrix of bar-height differences % (i.e. average RDM-correlation -% differences) +% differences) % stats_p_r.candDifferences_p: matrix of p values for all pairwise % candidate comparisons % stats_p_r.orderedCandidateRDMnames: candidate RDM names in the @@ -462,9 +462,9 @@ % end % ceilingUpperBound=mean(rs_ceilUpper); % ceilingLowerBound=mean(rs_ceilLower_LOO); -% +% % stats_p_r.ceiling = [ceilingLowerBound,ceilingUpperBound]; -% +% [ceilingUpperBound, ceilingLowerBound, bestFitRDM]=ceilingAvgRDMcorr(refRDM_stack,userOptions.RDMcorrelationType,false); stats_p_r.ceiling = [ceilingLowerBound,ceilingUpperBound]; @@ -526,7 +526,7 @@ fprintf('Averaging all instances of each candidate RDM.\n'); else fprintf('Found less than 12 instances of the reference RDM.\n'); - + if std(nsCandRDMinstances)==0 && nsCandRDMinstances(1)>=12 fprintf('Found %d instances of all candidate RDMs. Using these for random-effects inference.\n',nsCandRDMinstances(1)); fprintf('Averaging all instances of the reference RDM.\n'); @@ -547,10 +547,10 @@ % decide RDM-relatedness test if isequal(userOptions.RDMrelatednessTest,'subjectRFXsignedRank') || isequal(userOptions.RDMrelatednessTest,'subjectRFXbootstrap') || isequal(userOptions.RDMrelatednessTest,'subjectConditionRFXbootstrap') % user requested subject random-effects inference (by signed-rank or bootstrap test) - + if isequal(subjectRFXacross,'refRDMinstances'); fprintf('Using %d instances of the reference RDM for random-effects test of RDM relatedness.\n',nRefRDMinstances); - + elseif isequal(subjectRFXacross,'candRDMinstances') % fprintf('Found less than 12 instances of the reference RDM.'); % % we are already saying this earlier @@ -571,40 +571,51 @@ nSubjects = size(cand2refSims,1); es = std(cand2refSims)/sqrt(nSubjects); es = es(sortedIs); - - - + + + case 'subjectRFXbootstrap' fprintf('\nPerforming subject bootstrap test of RDM relatedness (subject as random effect).\n'); userOptions.resampleSubjects=1;userOptions.resampleConditions=0; - - [realRs bootstrapEs pairwisePs_bootSubj bootstrapRs] = bootstrapRDMs(refRDM_stack, meanCandRDMs, userOptions); + + [realRs bootstrapEs pairwisePs_bootSubj bootstrapRs] = ... + bootstrapRDMs( ... + refRDM_stack, meanCandRDMs, userOptions); + % bootstrapRs is now nCandRDMs x nBootstrap - bootstrapRs=bootstrapRs';% transpose to make it nBootstrap x nCandRDMs + bootstrapRs=bootstrapRs';% transpose to make it nBootstrap + % x nCandRDMs for candI = 1:numel(candRDMs) - [ps(candI)] = signrank_onesided(bootstrapRs(:,candI)); + [ci_lo(candI), ci_up(candI), ps(candI)] = bootConfint( ... + cand2refSims(:,candI), bootstrapRs(:,candI), 'greater', ... + userOptions); end ps = ps(sortedIs); stats_p_r.candRelatedness_p = ps; es = std(bootstrapRs); es = es(sortedIs); - - + + case 'subjectConditionRFXbootstrap' fprintf('\nPerforming subject and condition bootstrap test of RDM relatedness (subject and condition as random effects).\n'); userOptions.resampleSubjects=1;userOptions.resampleConditions=1; - - [realRs bootstrapEs pairwisePs_bootSubjCond bootstrapRs] = bootstrapRDMs(refRDM_stack, meanCandRDMs, userOptions); + + [realRs bootstrapEs pairwisePs_bootSubj bootstrapRs] = ... + bootstrapRDMs( ... + refRDM_stack, meanCandRDMs, userOptions); + % bootstrapRs is now nCandRDMs x nBootstrap bootstrapRs=bootstrapRs';% nBootstrap x nCandRDMs for candI = 1:numel(candRDMs) - [ps(candI)] = signrank_onesided(bootstrapRs(:,candI)); + [ci_lo(candI), ci_up(candI), ps(candI)] = bootConfint( ... + cand2refSims(:,candI), bootstrapRs(:,candI), 'greater', ... + userOptions); end ps = ps(sortedIs); stats_p_r.candRelatedness_p = ps; es = std(bootstrapRs);es = es(sortedIs); - - + + case 'randomisation' fprintf('\nPerforming condition-label randomisation test of RDM relatedness (fixed effects).\n'); nRandomisations = userOptions.nRandomisations; @@ -623,7 +634,7 @@ end%if % make space for null-distribution of correlations rs_null=nan(userOptions.nRandomisations,numel(candRDMs)); - + % index method would require on the order of n^2*nRandomisations % memory, so i'll go slowly for now... %tic @@ -634,7 +645,7 @@ else randomIndexSeq = randomPermutation(n); end%if - + rdmA_rand_vec=vectorizeRDM(meanRefRDM(randomIndexSeq,randomIndexSeq)); for candI = 1:numel(candRDMs) rs_null(randomisationI,candI)=rankCorr_Kendall_taua(rdmA_rand_vec',rdms(candI,:)'); @@ -652,7 +663,7 @@ else randomIndexSeq = randomPermutation(n); end%if - + rdmA_rand_vec=vectorizeRDM(meanRefRDM(randomIndexSeq,randomIndexSeq)); for candI = 1:numel(candRDMs) rs_null(randomisationI,candI)=raeSpearmanCorr(rdmA_rand_vec',rdms(candI,:)'); @@ -670,7 +681,7 @@ else randomIndexSeq = randomPermutation(n); end%if - + rdmA_rand_vec=vectorizeRDM(meanRefRDM(randomIndexSeq,randomIndexSeq)); rs_null(randomisationI,:)=corr(rdmA_rand_vec',rdms','type',userOptions.RDMcorrelationType,'rows','pairwise'); if mod(randomisationI,floor(userOptions.nRandomisations/100))==0 @@ -680,7 +691,7 @@ end % randomisationI fprintf('\n'); end - + % p-values from the randomisation test for candI = 1:numel(candRDMs) p_randCondLabels(candI) = 1 - relRankIn_includeValue_lowerBound(rs_null(:,candI),y(candI)); % conservative @@ -692,7 +703,7 @@ end p_randCondLabels_fwe = p_randCondLabels_fwe(sortedIs); stats_p_r.candRelatedness_p = [p_randCondLabels;p_randCondLabels_fwe]; - + otherwise fprintf('Not performing any test of RDM relatedness.\n'); end @@ -704,10 +715,10 @@ switch subjectRFXacross case 'refRDMinstances'; fprintf('Using %d instances of the reference RDM for random-effects tests comparing pairs of candidate RDMs.\n',nRefRDMinstances); - + case 'candRDMinstances' fprintf('Using %d instances of all candidate RDMs for random-effects tests comparing pairs of candidate RDMs.\n',nsCandRDMinstances(1)); - + otherwise fprintf('Attempting to revert to condition-bootstrap tests comparing pairs of candidate RDMs.\n'); if conditionRFX @@ -742,39 +753,75 @@ fprintf('Performing subject bootstrap test comparing candidate RDMs (subject as random effect).\n'); if ~exist('pairwisePs_bootSubj','var') [realRs bootstrapEs pairwisePs_bootSubj bootstrapRs] = bootstrapRDMs(refRDM_stack, meanCandRDMs, userOptions); + for candI = 1:numel(candRDMs) + [ci_lo(candI), ci_up(candI)] = bootConfint( ... + cand2refSims(:,candI), bootstrapRs(:,candI), 'greater', ... + userOptions); + end end hold on - errorbar(1:numel(candRDMs),Ys,bootstrapEs(sortedIs),'Color',[0 0 0],'LineWidth',... + errorbar(1:numel(candRDMs),y_sorted, y_sorted - ci_lo(sortedIs), ... + ci_up(sortedIs) - y_sorted,'Color',[0 0 0],'LineWidth',... 2,'LineStyle','none'); - - pairWisePs = pairwisePs_bootSubj; + + pairWisePs = nan(numel(candRDMs)); + for candI = 1:(numel(candRDMs)-1) + bsI = bootstrapRs(:,candI); + cand2refSimsI = cand2refSims(:,candI); + for candJ = (candI+1):numel(candRDMs) + bsJ = bootstrapRs(:,candJ); + cand2refSimsJ = cand2refSims(:,candJ); + [~, ~, pairWisePs(candI, candJ)] = bootConfint( ... + cand2refSimsI - cand2refSimsJ, bsI - bsJ, ... + 'two-tailed', userOptions); + pairWisePs(candJ, candI) = pairWisePs(candI, candJ); + end + end stats_p_r.candDifferences_p = pairWisePs(sortedIs,sortedIs); stats_p_r.SEs = bootstrapEs(sortedIs); - + case 'subjectConditionRFXbootstrap' fprintf('Performing condition and subject bootstrap test comparing candidate RDMs (condition and subject as random effects).\n'); if ~exist('pairwisePs_bootSubjCond','var') [realRs bootstrapEs pairwisePs_bootSubjCond bootstrapRs] = bootstrapRDMs(refRDM_stack, meanCandRDMs, userOptions); end - + hold on errorbar(1:numel(candRDMs),y_sorted,bootstrapEs(sortedIs),'Color',[0 0 0],'LineWidth',... 2,'LineStyle','none'); - + pairWisePs = pairwisePs_bootSubjCond; stats_p_r.candDifferences_p = pairWisePs(sortedIs,sortedIs); stats_p_r.SEs = bootstrapEs(sortedIs); - + case 'conditionRFXbootstrap' fprintf('Performing condition bootstrap test comparing candidate RDMs (subject as random effect).\n'); userOptions.resampleSubjects=0;userOptions.resampleConditions=1; [realRs bootstrapEs pairwisePs bootstrapRs] = bootstrapRDMs(refRDM_stack, meanCandRDMs, userOptions); + for candI = 1:numel(candRDMs) + [ci_lo(candI), ci_up(candI)] = bootConfint( ... + cand2refSims(:,candI), bootstrapRs(:,candI), 'greater', ... + userOptions); + end + pairWisePs = nan(numel(candRDMs)); + for candI = 1:(numel(candRDMs)-1) + bsI = bootstrapRs(:,candI); + cand2refSimsI = cand2refSims(:,candI); + for candJ = (candI+1):numel(candRDMs) + bsJ = bootstrapRs(:,candJ); + cand2refSimsJ = cand2refSims(:,candJ); + [~, ~, pairWisePs(candI, candJ)] = bootConfint( ... + cand2refSimsI - cand2refSimsJ, bsI - bsJ, ... + 'two-tailed', userOptions); + pairWisePs(candJ, candI) = pairWisePs(candI, candJ); + end + end stats_p_r.candDifferences_p = pairwisePs(sortedIs,sortedIs); stats_p_r.SEs = bootstrapEs(sortedIs); hold on errorbar(1:numel(candRDMs),y_sorted,bootstrapEs(sortedIs),'Color',[0 0 0],'LineWidth',... 2,'LineStyle','none'); - + otherwise fprintf('Not performing any test for comparing pairs of candidate RDMs.\n'); end @@ -799,18 +846,18 @@ %% add the p-values from RDM relatedness tests if ~isequal(userOptions.RDMrelatednessTest,'none') ps = stats_p_r.candRelatedness_p(1,:); - + pStringCell = cell(1, numel(candRDMs)); - + if isequal(userOptions.RDMrelatednessMultipleTesting,'none') || isequal(userOptions.RDMrelatednessMultipleTesting,'FWE') thresh = userOptions.RDMrelatednessThreshold; elseif isequal(userOptions.RDMrelatednessMultipleTesting,'FDR') thresh = FDRthreshold(ps,userOptions.candRDMdifferencesThreshold); end - + if isequal(userOptions.RDMrelatednessMultipleTesting,'FWE') ps_fwe = stats_p_r.candRelatedness_p(2,sortedIs); - + for test = 1:numel(candRDMs) if isequal(userOptions.plotpValues,'*') if ps(test) <= ps_fwe(test) @@ -852,7 +899,7 @@ end%switch:pIndication end%for:test end % isequal(userOptions.RDMrelatednessTest.... - + for test = 1:numel(candRDMs) switch userOptions.plotpValues case '*' @@ -881,7 +928,7 @@ pMat(logical(eye(numel(candRDMs)))) = 0; allPairwisePs = squareform(pMat); threshold = FDRthreshold(allPairwisePs,thresh); - + case 'FWE' nTests = numel(candRDMs)*(numel(candRDMs)-1)/2; threshold = thresh/nTests; @@ -1056,4 +1103,4 @@ exportCurrentFigAsPostscript([userOptions.analysisName,'.ps'],userOptions); end -end%function \ No newline at end of file +end%function diff --git a/Engines/exhaustivePermutations.m b/Engines/exhaustivePermutations.m new file mode 100644 index 0000000..89154c5 --- /dev/null +++ b/Engines/exhaustivePermutations.m @@ -0,0 +1,115 @@ +% permutations = exhaustivePermutations(n[, cap]) +% +% n number of items to be permuted +% cap (optional) stop after finding this many +% permutations each row is a permutation; they'll be in lexicographic order +% +% Cai Wingfield 3-2010 + +function permutations = exhaustivePermutations(varargin) + + % The following algorithm generates the next permutation lexicographically after + % a given permutation. It changes the given permutation in-place. + % + % 1. Find the largest index j such that a[j] < a[j + 1]. If no such index + % exists, the permutation is the last permutation. + % 2. Find the largest index l such that a[j] < a[l]. Such a l exists and + % satisfies j < l, since j + 1 is such an index. + % 3. Swap a[j] with a[l]. + % 4. Reverse the sequence from a[j + 1] up to an including the final element + % a[n]. + % + % After step 1, one knows that all of the elements strictly after position j + % form a weakly decreasing sequence, so no permutation of these elements will + % make it advance in lexicographic order; to advance one must increase a[j]. + % Step 2 finds the smallest value a[l] to replace a[j] by, and swapping them in + % step 3 leaves the sequence after position j in weakly decreasing order. + % Reversing this sequence in step 4 then produces its lexicographically minimal + % permutation, and the lexicographic successor of the initial state for the + % whole sequence. + + switch nargin + case 1 + n = varargin{1}; + cap = inf; + case 2 + n = varargin{1}; + cap = varargin{2}; + otherwise + error('exhaustivePermutations:nargin', 'Only one or two arguments, please.'); + end%switch:nargin + + absoluteMaximum = 21; % Larger than this and n! is not accurately represented by a double. Does this matter? + + initialLexOrder = 1:n; + + %if n > absoluteMaximum, permutations = []; return; end%if + if n > absoluteMaximum, warning('exhaustivePermutations:tooMany', ['21! is the largest number accurately represented by a double.\nAttempting ' num2str(n) '!, therefore, may give\nunpredictable results.']); end%if + + % The first in the list is just the initial lexical order + permutations = initialLexOrder; + i = 1; + allPermutationsFound = false; + + while ~allPermutationsFound && i < cap + + currentPermutation = permutations(i,:); + [nextPermutation, allPermutationsFound] = getNextPermutation(currentPermutation); + permutations = [permutations; nextPermutation]; + i = i + 1; + + end%while:~allPermutationsFound + +end%function + +%%%%%%%%%%%%%%%%%% +%% Subfunctions %% +%%%%%%%%%%%%%%%%%% + +function [nextPermutation, allPermutationsFound] = getNextPermutation(currentPermutation) + + n = numel(currentPermutation); + + allPermutationsFound = false; + + % Search through the current permutation, beginning to end, to find the largest index j such that currentPermutation(j) < currentPermutation(j+1) + j = 0; + for jSearch = 1:n-1 + if currentPermutation(jSearch) < currentPermutation(jSearch+1) + j = jSearch; + end%if + end%for:jSearch + + % If no such j exists, the current permutation is the last permutation in lexicographical order. + % Otherwise, find the largest index k such that currentPermutation(k) > currentPermutation(j). + if j == 0 + + nextPermutation = []; + allPermutationsFound = true; + + else + + k = 0; + for kSearch = j:n + if currentPermutation(kSearch) > currentPermutation(j) + k = kSearch; + end%if + end%for:kSearch + + nextPermutation = currentPermutation; + + % Swap currentPermutation(j) and currentPermutation(k) + + nextPermutation_j = nextPermutation(j); + nextPermutation_k = nextPermutation(k); + + nextPermutation(j) = nextPermutation_k; + nextPermutation(k) = nextPermutation_j; + + % Reverse order of currentPermutation(j+1:end) + + nextPermutation(j+1:end) = nextPermutation(end:-1:j+1); + + end%if:j == 0 + +end%function