-
Notifications
You must be signed in to change notification settings - Fork 0
/
NCP_08_SSD.m
132 lines (105 loc) · 5.4 KB
/
NCP_08_SSD.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
%% Script to calculate the Sum of Squared Diferences between diagnoses and groups (relatives, PE, FEP and chronic)
% Copyright (C) 2023 University of Seville
% Written by Natalia García San Martín (ngarcia1@us.es)
% This file is part of Neurobiology Centiles Psychosis toolkit.
%
% Neurobiological Centiles Mapping toolkit is free software:
% you can redistribute it and/or modify it under the terms of the
% GNU General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
%
% Neurobiological Centiles Mapping toolkit is distributed in the hope that
% it will be useful, but WITHOUT ANY WARRANTY; without even the implied
% warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Neurobiological Centiles Mapping toolkit. If not, see
% <https://www.gnu.org/licenses/>.
clear
close all
% Change as desired
location = 'D:\OneDrive - UNIVERSIDAD DE SEVILLA\Natalia\Code\Data\';
location = 'C:\Users\usuario\OneDrive - UNIVERSIDAD DE SEVILLA\Natalia\Code\Data\';
cd(location);
load("PE.mat",'-mat');
load("FEP.mat",'-mat');
load("SCZ.mat",'-mat');
load("Schizoaffective_Disorder.mat",'-mat')
dx_names = {'SCZ-relatives','SAD-relatives','PE-suspected','PE-definite','PE-severe', 'FEP','SCZ','SAD'};
names_group={'SCZ_Relative','Schizoaffective_Disorder_Relative','PE_1','PE_2','PE_3', ...
'FEP_1','SCZ','Schizoaffective_Disorder'};
nperm = 10000;
sse_real_total = nan(length(names_group));
sse_real_total = array2table(sse_real_total,'RowNames',names_group,'VariableNames',names_group);
PCT_pval = nan(length(names_group));
PCT_pval = array2table(PCT_pval,'RowNames',names_group,'VariableNames',names_group);
for i = 1:1:length(names_group)-1
for j = i+1:length(names_group)
[names_dx_1,names_dx_2] = names_group{[i,j]};
centiles_dx_1 = eval(['centiles_',names_dx_1]);
centiles_dx_2 = eval(['centiles_',names_dx_2]);
% SSE
sse_real = mse(mean(table2array(centiles_dx_1)),mean(table2array(centiles_dx_2)))*size(centiles_dx_1,2);
for iperm = 1:nperm
% Group membership is randomly reassigned across the compared diagnoses
[selection_dx_1,selection_dx_2,ratio_1(iperm),ratio_2(iperm)] = mix_dx(centiles_dx_1,centiles_dx_2);
selection_dx_1_m = mean(selection_dx_1);
selection_dx_2_m = mean(selection_dx_2);
% Permuted SSE
sse_perm(iperm) = mse(selection_dx_1_m,selection_dx_2_m)*size(selection_dx_1,2);
end
% Test if it explains more variance that expected by chance
PCT_pval{i,j} = sum(sse_perm > sse_real) / nperm;
sse_real_total{i,j} = sse_real;
end
end
figure;
sse_real_total_tri = triu(sse_real_total{:,:}) + triu(sse_real_total{:,:},1).';
heatmap(sse_real_total_tri,'XData',strrep(dx_names,'_',' '),'YData',strrep(dx_names,'_',' '));
title('SSE (centiles)')
figure;
PCT_pval_tri = triu(PCT_pval{:,:}) + triu(PCT_pval{:,:},1).';
heatmap(PCT_pval_tri ,'XData',strrep(dx_names,'_',' '),'YData',strrep(dx_names,'_',' '));
title('p-value (centiles)')
%% Groups %%
load("G0.mat",'-mat');
load("G1.mat",'-mat');
load("G1_5.mat",'-mat');
load("G2.mat",'-mat');
groups_names = {'SCZ_and_SAD-relatives','PE','FEP','SCZ_and_SAD-chronic'};
names_group = {'G0','G1','G1_5','G2'};
sse_real_total_groups = nan(length(names_group));
sse_real_total_groups = array2table(sse_real_total_groups,'RowNames',names_group,'VariableNames',names_group);
PCT_pval_groups = nan(length(names_group));
PCT_pval_groups = array2table(PCT_pval_groups,'RowNames',names_group,'VariableNames',names_group);
for i = 1:1:length(names_group)-1
for j = i+1:length(names_group)
[names_dx_1,names_dx_2] = names_group{[i,j]};
centiles_dx_1 = eval(['centiles_',names_dx_1]);
centiles_dx_2 = eval(['centiles_',names_dx_2]);
% SSE
sse_real_groups = mse(mean(table2array(centiles_dx_1)),mean(table2array(centiles_dx_2)))*size(centiles_dx_1,2);
for iperm = 1:nperm
% Group membership is randomly reassigned across the compared groups
[selection_dx_1,selection_dx_2,ratio_1(iperm),ratio_2(iperm)] = mix_dx(centiles_dx_1,centiles_dx_2);
selection_dx_1_m = mean(selection_dx_1);
selection_dx_2_m = mean(selection_dx_2);
% SSE
sse_perm_groups(iperm) = mse(selection_dx_1_m,selection_dx_2_m)*size(selection_dx_1,2);
end
% Test if it explains more variance that expected by chance
PCT_pval_groups{i,j} = sum(sse_perm_groups > sse_real_groups) / nperm;
sse_real_total_groups{i,j} = sse_real_groups;
end
end
figure;
sse_real_total_groups_tri = triu(sse_real_total_groups{:,:}) + triu(sse_real_total_groups{:,:},1).';
h = heatmap(sse_real_total_groups_tri,'XData',strrep(groups_names,'_',' '),'YData',strrep(groups_names,'_',' '),'FontSize',20);
h.CellLabelFormat = '%.3f';
title('SSE (centiles)')
figure;
PCT_pval_groups_tri = triu(PCT_pval_groups{:,:}) + triu(PCT_pval_groups{:,:},1).';
heatmap(PCT_pval_groups_tri,'XData',strrep(groups_names,'_',' '),'YData',strrep(groups_names,'_',' '));
title('p-value (centiles)')
'DONE'