-
Notifications
You must be signed in to change notification settings - Fork 1
/
demo_projpca.m
154 lines (139 loc) · 4.12 KB
/
demo_projpca.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
run('../env.m');
dataset1 = load([GEN_DATAPATH filesep '3749595figshare/' filesep 'empirical_HCP']);
dataset2 = load([GEN_DATAPATH filesep '3749595figshare/' filesep 'empirical_Paris']);
%dataset2 = load([GEN_DATAPATH filesep '3749595figshare/' filesep 'empirical_HCP']);
% Create data matrix Y and covariate matrix X;
p = length(dataset1.TS) + length(dataset2.TS);
n_rois = size(dataset1.TS{1},1);
n_volumes = min([size(dataset1.TS{1},2) size(dataset2.TS{1},2) 300]);
Y = zeros(p,n_rois,n_volumes);
X = zeros(p,1,n_volumes);
n_sites = 2;
standardize_rois = true;
for vol_no = 1:n_volumes;
tmp_Y = zeros(p,n_rois);
p1 = length(dataset1.TS);
for subj_no = 1:p1;
tmp_Y(subj_no,:) = dataset1.TS{subj_no}(:,vol_no);
end
if(standardize_rois)
std_y = std(tmp_Y(1:p1,:));
tmp_Y(1:p1,:) = bsxfun(@rdivide,tmp_Y(1:p1,:),std_y);
end
p2 = length(dataset2.TS);
for subj_no = 1:p2
tmp_Y(p1+subj_no,:) = dataset2.TS{subj_no}(:,vol_no);
end
if(standardize_rois)
std_y = std(tmp_Y(p1+1:p,:));
tmp_Y(p1+1:p,:) = bsxfun(@rdivide,tmp_Y(p1+1:p,:),std_y);
end
tmp_X(1:p1,1) = 1;
tmp_X(p1+1:p,1) = 2;
Y(:,:,vol_no) = tmp_Y;
X(:,1:n_sites,vol_no) = dummyvar(tmp_X);
clear tmp_Y tmp_X;
end
% Correction per volume
results = {};
mu = zeros(n_volumes,n_rois);
B = zeros(n_volumes,2,n_rois);
Yres = zeros(size(Y));
U = zeros(size(Y));
for vol_no=1:n_volumes;
results{vol_no} = projpca(Y(:,:,vol_no),X(:,1:2,vol_no));
B(vol_no,:,:) = results{vol_no}.beta;
mu(vol_no,:) = results{vol_no}.mu;
Yres(:,:,vol_no) = results{vol_no}.Y;
U(:,:,vol_no) = results{vol_no}.U;
end
figure;
subplot(3,2,1);
imagesc(corr(squeeze(Y(1,:,:))')); axis image equal;
title('Original Site 1')
subplot(3,2,2);
imagesc(corr(squeeze(Y(end,:,:))')); axis image equal;
title('Original Site 2')
subplot(3,2,3);
imagesc(corr(squeeze(Yres(1,:,:))')); axis image equal;
title('Res + mean Site 1')
subplot(3,2,4);
imagesc(corr(squeeze(Yres(end,:,:))')); axis image equal;
title('Res + mean Site 2')
subplot(3,2,5);
imagesc(corr(squeeze(U(1,:,:))')); axis image equal;
title('Res Site 1')
subplot(3,2,6);
imagesc(corr(squeeze(U(end,:,:))')); axis image equal;
title('Res Site 2')
% Plot mean site effects
figure;
roi_idx = 2
subplot(3,2,1);
plot(squeeze(Y(1,roi_idx,:)));
title('Original Site 1')
subplot(3,2,2);
plot(squeeze(Y(end,roi_idx,:)));
title('Original Site 2')
subplot(3,2,3);
plot(squeeze(Yres(1,roi_idx,:)));
title('Res + mean Site 1')
subplot(3,2,4);
plot(squeeze(Yres(end,roi_idx,:)));
title('Res + mean Site 2')
subplot(3,2,5);
plot(squeeze(U(1,roi_idx,:)));
title('Res Site 1')
subplot(3,2,6);
plot(squeeze(U(end,roi_idx,:)));
title('Res Site 2')
% % Shared correction over all volumes
% % This approach is not recommended.
% % When rois are scaled to have unit variance, this also does nothing.
% results = projpca(reshape(Y,[p*n_volumes n_rois]),reshape(X,[p*n_volumes n_sites]));
% mu = results.mu;
% B = results.beta;
% Yres = reshape(results.Y,[p n_rois n_volumes]);
% U = reshape(results.U,[p n_rois n_volumes]);
%
% figure;
% subplot(3,2,1);
% imagesc(corr(squeeze(Y(1,:,:))')); axis image equal;
% title('Original Site 1')
% subplot(3,2,2);
% imagesc(corr(squeeze(Y(end,:,:))')); axis image equal;
% title('Original Site 2')
% subplot(3,2,3);
% imagesc(corr(squeeze(Yres(1,:,:))')); axis image equal;
% title('Res + mean Site 1')
% subplot(3,2,4);
% imagesc(corr(squeeze(Yres(end,:,:))')); axis image equal;
% title('Res + mean Site 2')
% subplot(3,2,5);
% imagesc(corr(squeeze(U(1,:,:))')); axis image equal;
% title('Res Site 1')
% subplot(3,2,6);
% imagesc(corr(squeeze(U(end,:,:))')); axis image equal;
% title('Res Site 2')
%
% % Plot mean site effects
% figure;
% roi_idx = 2
% subplot(3,2,1);
% plot(squeeze(Y(1,roi_idx,:)));
% title('Original Site 1')
% subplot(3,2,2);
% plot(squeeze(Y(end,roi_idx,:)));
% title('Original Site 2')
% subplot(3,2,3);
% plot(squeeze(Yres(1,roi_idx,:)));
% title('Res + mean Site 1')
% subplot(3,2,4);
% plot(squeeze(Yres(end,roi_idx,:)));
% title('Res + mean Site 2')
% subplot(3,2,5);
% plot(squeeze(U(1,roi_idx,:)));
% title('Res Site 1')
% subplot(3,2,6);
% plot(squeeze(U(end,roi_idx,:)));
% title('Res Site 2')