-
Notifications
You must be signed in to change notification settings - Fork 1
/
spmrt_corr.m
174 lines (146 loc) · 6.35 KB
/
spmrt_corr.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
function [rP,CIP,rC,CIC]=spmrt_corr(image1,image2,mask,metric,figout,threshold,alpha_level)
% Computes the Pearson (vector-wise) and concordance correlations
% between image1 and image2 for data included in the mask
%
% FORMAT [rP,rC,CIP,CIC]=spmup_corr(image1,image1,'both',mask,threshold,alpha_level)
% [rP,CIP]=spmup_corr(image1,image1,'Pearson',mask,threshold,alpha_level)
% [~,~,rC,CIC]=spmup_corr(image1,image1,'Concordance',mask,threshold,alpha_level)
%
% INPUT if no input the user is prompted
% image1 is the filename of an image (see spm_select)
% image2 is the filename of an image (see spm_select)
% mask is the filename of a mask in same space as image1 and image2
% metric is 'Pearson', 'Concordance', or 'both'
% figout 1/0 (default) to get correlation figure out
% threshold (optional) if mask is not binary, threshold to apply
% alpha_level is the level used to compute the confidence interval (default is 5%)
%
% OUTPUT rP is the Pearson correlation coefficient
% rC is the concordance correlation coefficient
% CIP is the 95% confidence interval of rP (if 0 not included then significant)
% CIC is the 95% confidence interval of rC (if 0 not included then significant)
%
% Concordance correlation is more useful for reliability because it estimates
% how much variation from the 45 degree line we have (by using the cov)
% see Lin, L.I. 1989. A Corcordance Correlation Coefficient to Evaluate
% Reproducibility. Biometrics 45, 255-268.
%
% Cyril Pernet
% --------------------------------------------------------------------------
% Copyright (C) spmrt
rP = []; CIP = [];
rC = []; CIC = [];
nboot = 1000; % a thousand bootstraps
%% check inputs
spm('defaults', 'FMRI');
if nargin < 7; threshold = []; alpha_level = 5/100; end
if nargin < 5; figout = 0; end
if nargin < 4; metric = 'both'; end
if nargin == 0
[image1,image2,mask]=spmrt_pickupfiles;
if size(mask,1) > 1
warning('only one mask image allowed, using the 1st frame only')
mask = mask(1,:);
end
figout = 1; % if user if prompt to enter data then return a figure
end
%% Get the data
if exist('threshold','var') && ~isempty(threshold)
X = spmrt_getdata(image1,image2,mask,threshold);
else
X = spmrt_getdata(image1,image2,mask);
end
n = size(X,1);
if any(sum(X,1) == 0)
error('at least one image is empty')
end
%% Pearson correlation
if strcmpi(metric,'Pearson') || strcmpi(metric,'Both')
rP = nansum(detrend(X(:,1),'constant').*detrend(X(:,2),'constant')) ./ ...
(nansum(detrend(X(:,1),'constant').^2).*nansum(detrend(X(:,2),'constant').^2)).^(1/2);
if nargout > 1
disp('computing Pearson''s CI'); go = 1;
while go == 1
table = randi(n,n,nboot);
rPB = nansum(detrend(reshape(X(table,1),[n nboot]),'constant').*detrend(reshape(X(table,2),[n nboot]),'constant')) ./ ...
(nansum(detrend(reshape(X(table,1),[n nboot]),'constant').^2).*nansum(detrend(reshape(X(table,2),[n nboot]),'constant').^2)).^(1/2);
rPB = sort(rPB);
rPB(isnan(rPB)) = [];
adj_nboot = length(rPB);
low = round((alpha_level*adj_nboot)/2); % lower bound
high = adj_nboot - low; % upper bound
CIP = [rPB(low) rPB(high)];
if rP>CIP(1) && rP<CIP(2)
go = 0; % this can happen that resampling is weird and rp not in the bounds
end
end
end
end
%% Concordance
if strcmpi(metric,'Concordance') || strcmpi(metric,'Both')
S = cov(X,1);
ybar = mean(X);
shiftC = (ybar(1)-ybar(2))/(sqrt(sqrt(S(1,1))*sqrt(S(2,2))));
scaleC = sqrt(S(1,1))/sqrt(S(2,2));
biasFactorC = ((scaleC+1/scaleC+shiftC^2)/(2))^-1;
Var1 = S(1,1); Var2 = S(2,2); S = S(1,2);
rC = (2.*S) ./ ( Var1 +Var2 + (ybar(1)-ybar(2))^2);
% Assign these values to the caller workspace
assignin('caller','shiftC',shiftC);
assignin('caller','scaleC',scaleC);
assignin('caller','biasFactorC',biasFactorC);
if nargout > 1
disp('computing Concordance correlation CI'); go = 1;
while go == 1
if strcmpi(metric,'Concordance')
table = randi(n,n,nboot); % otherwise reuse the one from above = same sampling scheme
end
rCB = zeros(1,nboot);
for b=1:nboot
S = cov(X(table(:,b),:),1); Var1 = S(1,1); Var2 = S(2,2); S = S(1,2);
rCB(b) = (2.*S) ./ ( Var1 +Var2 + (mean(X(table(:,b),1))-mean(X(table(:,b),2)))^2);
end
rCB = sort(rCB);
adj_nboot = nboot - sum(isnan(rCB));
low = round((alpha_level*adj_nboot)/2); % lower bound
high = adj_nboot - low; % upper bound
rCB(isnan(rCB)) = [];
CIC = [rCB(low) rCB(high)];
if rC>CIC(1) && rC<CIC(2)
go = 0;
end
end
end
end
%% figure
if figout == 1
figure('Name','images correlation')
set(gcf,'Color','w','InvertHardCopy','off', 'units','normalized','outerposition',[0 0 1 1])
scatter(X(:,1),X(:,2),50); grid on % plot pairs of observations
xlabel('img1','FontSize',14); ylabel('img2','FontSize',14); % label
box on; set(gca,'Fontsize',14); axis square; hold on
v = axis;
if strcmp(metric,'Pearson') || strcmp(metric,'both')
h=lsline; set(h,'Color','b','LineWidth',4); % add the least square line
end
if strcmp(metric,'Concordance') || strcmp(metric,'both')
idxMax = find(v==max(abs(v)));
idxMin = find(v==min(abs(v)));
identity = v(idxMin):v(idxMax);
idplot = plot(identity,identity,'k--','LineWidth',2); % Identity line
concplot = plot(identity,identity*scaleC + shiftC,'r','LineWidth',4);
end
if strcmpi(metric,'Pearson')
title(['Pearson corr =' num2str(rP)],'FontSize',16);
legend(h, ...
{'Pearson'});
elseif strcmpi(metric,'Concordance')
title(['Concordance corr =' num2str(rC) 'Bias factor: ' num2str(biasFactorC)],'FontSize',16);
legend([idplot,concplot], ...
{'Identity Line', 'Concordance'});
else
title(sprintf('Pearson corr =%g \n Concordance corr =%g Bias factor=%g ', rP,rC,biasFactorC),'FontSize',16);
legend([idplot,concplot,h], ...
{'Identity Line', 'Concordance','Pearson'});
end
end