%% Demographics score plots

clear;clc;close all

% Change the directory to where the main model is,
% abcdAll_noA\framework\spls_fair20-5_fair_modeA_pub_noAnom_new_common_bmi.

base_dir = 'C:\Users\Konstantinos Tsirlis\Dropbox\My Desktop\Mourao - Miranda\ABCD results\abcdAll_noA';
framework_dir = [base_dir '\framework\spls_fair20-5_fair_modeA_pub_noAnom_new_common_bmi'];

load([framework_dir '\outmat_1.mat'])


% In our case we created a table with a few different demographics that are
% presented in the paper. If you rerun the example with a different subject
% list keep in mind that you need to respect this structure if you want the
% script to run appropriately. The data in there are not real and are there
% for illustration purposes only. If you want the actual demographics of
% the subjects you need to request permission to download them directly
% from the ABCD.
load('C:\Users\Konstantinos Tsirlis\Dropbox\My Desktop\Mourao - Miranda\ABCD results\ABCD Github Repo\demographics_checkup.mat')
sex = pdem_dummy.sex; male = sex==0; female = sex==1;
white = pdem_dummy.white; white = logical(white);
black = pdem_dummy.black; black = logical(black);
other = pdem_dummy.other; other = logical(other);
socioecon_top = pdem_dummy.socioecon==2;
socioecon_bot = pdem_dummy.socioecon==1;



% NEW CHANGES to the groupings above.
load([base_dir '\data\no_anom_new_common_bmi\Y_Proj.mat'])

% 1) we will not use other in general from the analysis.

% 2) we build the hispanic group.
hispanic = Y(:,281);
hispanic(hispanic==2) = 0;
hispanic(isnan(hispanic)) = 0;
hispanic(hispanic&black) = 0;
hispanic(hispanic&other) = 0;
hispanic = logical(hispanic);

% 3) so we exclude from the white group all the people who replied yes to the question "do you consider your child Hispanic?"
white = white & Y(:,281)==2;



x_lim = [-60 100; -135 100; -160 200; -300 100; -220 150; -160 140];
y_lim = [-14 20;  -9 7;     -7 11;    -6 3;     -15 10; -11 14];
location = ['nw'; 'nw'; 'nw'; 'nw'; 'nw'; 'nw'];




for i = 1:6 % number of signiticant modes, change accordingly in case of rerun.

    
    
% cd to the mode whose scores you want to plot - in our case mode1 of the main analysis
cd([framework_dir '\res\level' num2str(i)])
load('P.mat')




% SEX

% figure
figure('Units', 'Normalized', 'OuterPosition', [0 0 0.5 0.78]);
pos = [0.3 0.2 0.6 0.56]; % position same across different types of plots (e.g. 2d_group, 2d_cmap)
axes('Position', pos,'linewidth',2)
hold on;

% Plot data
scatter(P(male&oteid,1), P(male&oteid,2), 120, 'MarkerFaceColor', [0.6350 0.0780 0.1840], 'MarkerEdgeColor', 'k');

set(gca, 'FontSize', 40, 'FontName', 'Arial');
xlabel('Brain Score', 'FontSize', 43, 'fontname', 'Arial');
ylabel('Psychosocial Score', 'FontSize', 43, 'fontname', 'Arial');

h(1) = plot(nan, nan, 'o', 'MarkerSize', 16, 'MarkerFaceColor', [0.6350 0.0780 0.1840], 'MarkerEdgeColor', 'k', 'DisplayName', 'Male');
legend(h,'location', location(i,:)); clear h

xlim(x_lim(i,:));
ylim(y_lim(i,:));

fname = ['mode' num2str(i) '_proj2_sex_male'];

saveas(gcf, [fname '.png']);


% figure
figure('Units', 'Normalized', 'OuterPosition', [0 0 0.5 0.78]);
pos = [0.3 0.2 0.6 0.56]; % position same across different types of plots (e.g. 2d_group, 2d_cmap)
axes('Position', pos,'linewidth',2)
hold on;

% Plot data
scatter(P(female&oteid,1), P(female&oteid,2), 120, 'MarkerFaceColor', [0 0.4470 0.7410], 'MarkerEdgeColor', 'k');

set(gca, 'FontSize', 40, 'FontName', 'Arial');
xlabel('Brain Score', 'FontSize', 43, 'fontname', 'Arial');
ylabel('Psychosocial Score', 'FontSize', 43, 'fontname', 'Arial');

h(1) = plot(nan, nan, 'o', 'MarkerSize', 16, 'MarkerFaceColor', [0 0.4470 0.7410], 'MarkerEdgeColor', 'k', 'DisplayName', 'Female');
legend(h,'location', location(i,:)); clear h

xlim(x_lim(i,:));
ylim(y_lim(i,:));

fname = ['mode' num2str(i) '_proj2_sex_female'];

saveas(gcf, [fname '.png']);







% RACE

% figure
figure('Units', 'Normalized', 'OuterPosition', [0 0 0.5 0.78]);
pos = [0.3 0.2 0.6 0.56]; % position same across different types of plots (e.g. 2d_group, 2d_cmap)
axes('Position', pos,'linewidth',2)
hold on;


% Plot data
scatter(P(white&oteid,1), P(white&oteid,2), 120, 'MarkerFaceColor', [0.6350 0.0780 0.1840], 'MarkerEdgeColor', 'k'); hold on,

set(gca, 'FontSize', 40, 'FontName', 'Arial');
xlabel('Brain Score', 'FontSize', 43, 'fontname', 'Arial');
ylabel('Psychosocial Score', 'FontSize', 43, 'fontname', 'Arial');


h(1) = plot(nan, nan, 'o', 'MarkerSize', 16, 'MarkerFaceColor', [0.6350 0.0780 0.1840], 'MarkerEdgeColor', 'k', 'DisplayName', 'White');
legend(h,'location', location(i,:)); clear h

xlim(x_lim(i,:));
ylim(y_lim(i,:));

fname = ['mode' num2str(i) '_proj2_race_white'];

saveas(gcf, [fname '.png']);


% figure
figure('Units', 'Normalized', 'OuterPosition', [0 0 0.5 0.78]);
pos = [0.3 0.2 0.6 0.56]; % position same across different types of plots (e.g. 2d_group, 2d_cmap)
axes('Position', pos,'linewidth',2)
hold on;


% Plot data
scatter(P(black&oteid,1), P(black&oteid,2), 120, 'MarkerFaceColor', [0 0.4470 0.7410], 'MarkerEdgeColor', 'k');


set(gca, 'FontSize', 40, 'FontName', 'Arial');
xlabel('Brain Score', 'FontSize', 43, 'fontname', 'Arial');
ylabel('Psychosocial Score', 'FontSize', 43, 'fontname', 'Arial');


h(1) = plot(nan, nan, 'o', 'MarkerSize', 16, 'MarkerFaceColor', [0 0.4470 0.7410], 'MarkerEdgeColor', 'k', 'DisplayName', 'Black');
legend(h,'location', location(i,:)); clear h

xlim(x_lim(i,:));
ylim(y_lim(i,:));

fname = ['mode' num2str(i) '_proj2_race_black'];

saveas(gcf, [fname '.png']);


% figure
figure('Units', 'Normalized', 'OuterPosition', [0 0 0.5 0.78]);
pos = [0.3 0.2 0.6 0.56]; % position same across different types of plots (e.g. 2d_group, 2d_cmap)
axes('Position', pos,'linewidth',2)
hold on;


% Plot data
scatter(P(hispanic&oteid,1), P(hispanic&oteid,2), 120, 'MarkerFaceColor', [0.4660 0.6740 0.1880], 'MarkerEdgeColor', 'k');


set(gca, 'FontSize', 40, 'FontName', 'Arial');
xlabel('Brain Score', 'FontSize', 43, 'fontname', 'Arial');
ylabel('Psychosocial Score', 'FontSize', 43, 'fontname', 'Arial');


h(1) = plot(nan, nan, 'o', 'MarkerSize', 16, 'MarkerFaceColor', [0.4660 0.6740 0.1880], 'MarkerEdgeColor', 'k', 'DisplayName', 'Hispanic');
legend(h,'location', location(i,:)); clear h

xlim(x_lim(i,:));
ylim(y_lim(i,:));

fname = ['mode' num2str(i) '_proj2_race_hispanic'];

saveas(gcf, [fname '.png']);






% SOCIOECON

% figure
figure('Units', 'Normalized', 'OuterPosition', [0 0 0.5 0.78]);
pos = [0.3 0.2 0.6 0.56]; % position same across different types of plots (e.g. 2d_group, 2d_cmap)
axes('Position', pos,'linewidth',2)
hold on;


% Plot data
scatter(P(socioecon_top&oteid,1), P(socioecon_top&oteid,2), 120, 'MarkerFaceColor', [0.6350 0.0780 0.1840], 'MarkerEdgeColor', 'k'); hold on,


set(gca, 'FontSize', 40, 'FontName', 'Arial');
xlabel('Brain Score', 'FontSize', 43, 'fontname', 'Arial');
ylabel('Psychosocial Score', 'FontSize', 43, 'fontname', 'Arial');


h(1) = plot(nan, nan, 'o', 'MarkerSize', 16, 'MarkerFaceColor', [0.6350 0.0780 0.1840], 'MarkerEdgeColor', 'k', 'DisplayName', 'SocEc top');
legend(h,'location', location(i,:)); clear h

xlim(x_lim(i,:));
ylim(y_lim(i,:));

fname = ['mode' num2str(i) '_proj2_socioecon_top'];

saveas(gcf, [fname '.png']);



% figure
figure('Units', 'Normalized', 'OuterPosition', [0 0 0.5 0.78]);
pos = [0.3 0.2 0.6 0.56]; % position same across different types of plots (e.g. 2d_group, 2d_cmap)
axes('Position', pos,'linewidth',2)
hold on;


% Plot data
scatter(P(socioecon_bot&oteid,1), P(socioecon_bot&oteid,2), 120, 'MarkerFaceColor', [0 0.4470 0.7410], 'MarkerEdgeColor', 'k');


set(gca, 'FontSize', 40, 'FontName', 'Arial');
xlabel('Brain Score', 'FontSize', 43, 'fontname', 'Arial');
ylabel('Psychosocial Score', 'FontSize', 43, 'fontname', 'Arial');


h(1) = plot(nan, nan, 'o', 'MarkerSize', 16, 'MarkerFaceColor', [0 0.4470 0.7410], 'MarkerEdgeColor', 'k', 'DisplayName', 'SocEc bot');
legend(h,'location', location(i,:)); clear h

xlim(x_lim(i,:));
ylim(y_lim(i,:));

fname = ['mode' num2str(i) '_proj2_socioecon_bot'];

saveas(gcf, [fname '.png']);


end



%% Stats


clear;clc

% Change the paths accordingly.
base_dir = 'C:\Users\Konstantinos Tsirlis\Dropbox\My Desktop\Mourao - Miranda\ABCD results\abcdAll_noA';
framework_dir = [base_dir '\framework\spls_fair20-5_fair_modeA_pub_noAnom_new_common_bmi'];


load('C:\Users\Konstantinos Tsirlis\Dropbox\My Desktop\Mourao - Miranda\ABCD results\demographics_checkup.mat')
sex = pdem_dummy.sex; male = sex==0; female = sex==1;
white = pdem_dummy.white; white = logical(white);
black = pdem_dummy.black; black = logical(black);
other = pdem_dummy.other; other = logical(other);
socioecon_top = pdem_dummy.socioecon==2;
socioecon_bot = pdem_dummy.socioecon==1;



% NEW CHANGES to the groupings above.
load([base_dir '\data\no_anom_new_common_bmi\Y_Proj.mat'])

% 1) we will not use other in general from the analysis.

% 2) we build the hispanic group.
hispanic = Y(:,281);
hispanic(hispanic==2) = 0;
hispanic(isnan(hispanic)) = 0;
hispanic(hispanic&black) = 0;
hispanic(hispanic&other) = 0;
hispanic = logical(hispanic);

% 3) so we exclude from the white group all the people who replied yes to the question "do you consider your child Hispanic?"
white = white & Y(:,281)==2;







for j = 1:6 % number of significant modes, change accordingly




clear correl_pval


res.dir.frwork = framework_dir;
res.frwork.level = j;
res = res_defaults(res, 'load');


% cd to the framework to load the outmat file

cd(framework_dir)

load_outmat = dir('outmat*');
load(load_outmat.name);


% cd to the level of the framework you want to work with (j)
% Load param_* of each level

cd([framework_dir filesep 'res' filesep 'level' num2str(j)])

load_param = dir('param*');
load(load_param.name);


% each level has a different P so you have to reload it accordingly.
load('P.mat');
P = reshape(P,[length(P),2]);


% Change S.correl of model_*
model_name = dir('model*'); model_name = model_name.name;


% male
correl = corr(P(oteid & male,1),P(oteid & male,2));
save(model_name,'correl','-append')
res = stat_inference(res);
correl_pval(1,1:2) = [correl res.stat.split.pval];

% female
correl = corr(P(oteid & female,1),P(oteid & female,2));
save(model_name,'correl','-append')
res = stat_inference(res);
correl_pval(2,1:2) = [correl res.stat.split.pval];

% white
correl = corr(P(oteid & white,1),P(oteid & white,2));
save(model_name,'correl','-append')
res = stat_inference(res);
correl_pval(3,1:2) = [correl res.stat.split.pval];

% black
correl = corr(P(oteid & black,1),P(oteid & black,2));
save(model_name,'correl','-append')
res = stat_inference(res);
correl_pval(4,1:2) = [correl res.stat.split.pval];

% hispanic
correl = corr(P(oteid & hispanic,1),P(oteid & hispanic,2));
save(model_name,'correl','-append')
res = stat_inference(res);
correl_pval(5,1:2) = [correl res.stat.split.pval];

% socioecon_top
correl = corr(P(oteid & socioecon_top,1),P(oteid & socioecon_top,2));
save(model_name,'correl','-append')
res = stat_inference(res);
correl_pval(6,1:2) = [correl res.stat.split.pval];

% socioecon_bot
correl = corr(P(oteid & socioecon_bot,1),P(oteid & socioecon_bot,2));
save(model_name,'correl','-append')
res = stat_inference(res);
correl_pval(7,1:2) = [correl res.stat.split.pval];

% Change back the model_* to the original ones.
copyfile(['original_' model_name],model_name)


correl_pval_final{j,1} = correl_pval;


end




%% Some subgroup correlations

% Cognition Total composite score age-corrected standard score
% Total combined family income
[correl_iq_income_white, pval_iq_income_white] = corr(behVars_99to1_10mv(white,:).demo_comb_income_v2,behVars_99to1_10mv(white,:).nihtbx_totalcomp_agecorrected,'rows','complete')
[correl_iq_income_black, pval_iq_income_black] = corr(behVars_99to1_10mv(black,:).demo_comb_income_v2,behVars_99to1_10mv(black,:).nihtbx_totalcomp_agecorrected,'rows','complete')

% Crystallized Composite Age-corrected Standard score
% Fluid Composite Age-corrected Standard score
[correl_crystal_fluid_white, pval_crystal_fluid_white] = corr(behVars_99to1_10mv(white,:).nihtbx_cryst_agecorrected,behVars_99to1_10mv(white,:).nihtbx_fluidcomp_agecorrected,'rows','complete')
[correl_crystal_fluid_black, pval_crystal_fluid_black] = corr(behVars_99to1_10mv(black,:).nihtbx_cryst_agecorrected,behVars_99to1_10mv(black,:).nihtbx_fluidcomp_agecorrected,'rows','complete')