% This is how to calculate the standard error and 95% confidence limits on
% the Spearman rank correlation coefficient using the nonparametric bootstrap method in Matlab. 
% Note that by default Matlab reports bias corrected percentile 
% (BCA) confidence limits; these are more accurate than those
% obtained using the standard percentile method. 

var1=[6.1;3.2;5.5;2.3;3.8;7.2;4.1;7.6;4.2;10.1];  % first variable
var2=[3.9;2.8;8.2;3.9;6.2;4.8;2.1;9.0;5.5;8.6];   % second variable
var1rt=Tiedrank(var1); % to calculate the Spearman rank correlation we will calculate a Pearson correlation on rank transformed variables
var2rt=Tiedrank(var2); % Tiedrank() rank transforms the variables
B=100000; % desired number of bootstrap replicates
[bootstat,bootsam]=bootstrp(B,@corr,var1rt,var2rt); % calculate the correlation coefficients on 10 000 bootstrap resamplings of the original variables (returned as bootstat, bootsam gives the actual resampled data)
figure % create a new graphics window so that the previous figure is not overwritten
hist(bootstat,16) % make a histogram of the correlation coefficient calculated on the bootstrap resampled rank transformed variables
title('Distribution of correlation coefficients calculated on bootstrap resampled data') % histogram title
SpearmanR=corr(var1rt,var2rt) % the correlation between the original two rank transformed variables
sterror=std(bootstat) % the standard error on the correlation coefficient calculated as the standard deviation of the correlation coefficient calculated on the bootstrap resampled rank transformed variables
conflimsBCA=bootci(B,{@corr,var1rt,var2rt},'alpha',0.05,'type','bca')  % the BCA corrected 95% bootstrap confidence limits
conflimsPERC=bootci(B,{@corr,var1rt,var2rt},'alpha',0.05,'type','per') % the 95% bootstrap confidence limits, calculated using the basic percentile method (in this case the result is the same)

% If we would like to test whether the correlation coefficient is 
% significantly greater than zero we can simply count the proportion 
% of the cases in which the correlation coefficient calculated on the 
% bootstrapped variables was smaller than zero. We can do this as follows:

onesidedpboots=mean(bootstat<0); % one-sided p-value
twosidedpboots=2*onesidedpboots; % two-sided p-value
bootstrresult=[SpearmanR,onesidedpboots,twosidedpboots]

% This would have been the one or two-sided p-levels obtained for the 
% Spearman Rank correlation using the standard (asymptotic) normal approximation.

[SpearmanR,onesidedpasympt]=corr(var1,var2,'type','Spearman','tail','gt'); % one-sided p-value
[SpearmanR,twosidedpasympt]=corr(var1,var2,'type','Spearman','tail','ne'); % two-sided p-value
asymptresult=[SpearmanR,onesidedpasympt,twosidedpasympt]

% An alternative way for calculating the p-level could have been
% to randomly permute the original variables to create a null
% distribution and compare the obtained Spearman rank correlation
% against this null distribution. We can do this as follows:

P=100000; % desired number of random permutations
corrperm = zeros(P,1); % we first ceate an empty array to store the null distribution correlations
for i=1:P; corrperm(i) = corr(var1rt,var2rt(randperm(numel(var2rt)))); end % the correlation coefficients calculated after randomly permuting both variables
onesidedperm=mean(corrperm>SpearmanR);      % one-sided p-value
twosidedperm=mean(abs(corrperm)>SpearmanR); % two-sided p-value
permutationresult=[SpearmanR,onesidedperm,twosidedperm]
figure % create a new graphics window so that the previous figure is not overwritten
hist(corrperm,16) % make a histogram of the obtained null distribution
title('Correlation coefficient null distribution calculated using random permutation') % histogram title

% Permutation tests return exact p-levels, and should in general be preferred
% over bootstrap tests or asymptotic results.
% In this case you can see that the confidence limits returned by all
% three methods are very comparable but that the p-levels returned by
% the bootstrap method is overly optimistic.