clear
clc
close all
load('Dr_Namazi.mat')
%% Divided by Gender - Application of Switch command +Case
j=1;
m=1;

for i=1:size(DrNamazi.Gender,1)
    A = DrNamazi.Gender(i);
    switch A
        case 1
            Women(j,1)=DrNamazi.CA(i);
            Women(j,2)=DrNamazi.CB(i);
            Women(j,3)=DrNamazi.CC(i);
            Women(j,4)=DrNamazi.CD(i);
            Women(j,5)=DrNamazi.Ctot(i);
            Women(j,6)=DrNamazi.DTotal(i);
            j=j+1;
        case 2
            men(m,1)=DrNamazi.CA(i);
            men(m,2)=DrNamazi.CB(i);
            men(m,3)=DrNamazi.CC(i);
            men(m,4)=DrNamazi.CD(i);
            men(m,5)=DrNamazi.Ctot(i);
            men(m,6)=DrNamazi.DTotal(i);
            m=m+1;
    end
end
%% Investigation associasion - Correlation
% Pearson Correlation between Men in different groups (CA,CB,CC,CD,CTot,DTot)
for i=1:size(men,2)
    a=men(:,i);
    for j=1:size(men,2)
        [rho,pval]=corr(a,men(:,j));% a = men(:,i)
        corrmen(i,j)=rho;
        sigmen(i,j)=pval;
%         corrmen(i,j)=corr(a,men(:,j));% a = men(:,i)
    end
end
figure
imagesc(corrmen)
set(gca,'YDir','normal') % Normal kardan Mehvar Y, Agar nabashe zahere badi dare
caxis([-1 1]) % Range of Color bar (Range of Correlation)
colormap jet % Colors of the plot / Best colors that we could select
title('Men ')
set(gca,'fontsize',14)
asami={'CA','CB','CC','CD','CTot','D/B'};
set(gca, 'XTick', 1:6, 'XTickLabel', asami);
set(gca, 'YTick', 1:6, 'YTickLabel', asami);
axis square
colorbar
%%%%%% For Women
for i=1:size(Women,2)
    b=Women(:,i);
    for j=1:size(Women,2)
        [rho,pval]=corr(b,Women(:,j));% b = Women(:,i)
        corrwomen(i,j)=rho;
        sigwomen(i,j)=pval;
    end
end
figure
imagesc(corrwomen)
set(gca,'YDir','normal')
caxis([-1 1])
colormap jet
title('Women ')
set(gca,'fontsize',14)
asami={'CA','CB','CC','CD','CTot','D/B'};
set(gca, 'XTick', 1:6, 'XTickLabel', asami);
set(gca, 'YTick', 1:6, 'YTickLabel', asami);
colorbar
axis square
% Clear extra variables
clear a A b i j m pval rho 
% Subplot fancy - Men/Women
figure
subplot(1,2,1)
imagesc(corrmen)
set(gca,'YDir','normal') % Normal kardan Mehvar Y, Agar nabashe zahere badi dare
caxis([-1 1]) % Range of Color bar (Range of Correlation)
colormap jet % Colors of the plot / Best colors that we could select
title('Men ')
set(gca,'fontsize',14)
asami={'CA','CB','CC','CD','CTot','D/B'};
set(gca, 'XTick', 1:6, 'XTickLabel', asami);
set(gca, 'YTick', 1:6, 'YTickLabel', asami);
axis square
colorbar
subplot(1,2,2)
imagesc(corrwomen)
set(gca,'YDir','normal')
caxis([-1 1])
colormap jet
title('Women ')
set(gca,'fontsize',14)
asami={'CA','CB','CC','CD','CTot','D/B'};
set(gca, 'XTick', 1:6, 'XTickLabel', asami);
set(gca, 'YTick', 1:6, 'YTickLabel', asami);
colorbar
axis square
%%%% Total population
Pop=table2array(DrNamazi(:,16:21));
for i=1:size(Pop,2)
    b=Pop(:,i);
    for j=1:size(Pop,2)
        [rho,pval]=corr(b,Pop(:,j));
        corrpop(i,j)=rho;
        sigpop(i,j)=pval;
    end
end
figure
imagesc(corrpop)
set(gca,'YDir','normal')
caxis([-1 1])
colormap jet
title('Population ')
set(gca,'fontsize',14)
asami={'CA','CB','CC','CD','CTot','D/B'};
set(gca, 'XTick', 1:6, 'XTickLabel', asami);
set(gca, 'YTick', 1:6, 'YTickLabel', asami);
colorbar
axis square
%% One way Anova
asami1={'CA','CB','CC','CD','CTot','DTot'};
[pp,tbl,stats]=anova1(Pop,asami1);
figure
boxplot(Pop,asami1);
hold on
plot(mean(Pop), 'dg')
hold off
figure
% [c,~,~,gnames] = multcompare(stats);
[c, m, h, gnames] = multcompare(stats);
title('Multi comparison-Post HOC test')
%%%
% x1=pop(:,1);
% x2=pop(:,2);
% x3=pop(:,3);
% x4=pop(:,4);
% x5=pop(:,5);
% x6=pop(:,6);
% data = [x1' x2' x3' x4' x5' x6']; %// Create row vector with your data
% name={'CA','CB','CC','CD','CTot','DTot'};
% C1 = cell(1, length(x1));
% C1(:) = name(1);
% C2 = cell(1, length(x2));
% C2(:) = name(2);
% C3 = cell(1, length(x3));
% C3(:) = name(3);
% C4 = cell(1, length(x4));
% C4(:) = name(4);
% C5 = cell(1, length(x5));
% C5(:) = name(5);
% C6 = cell(1, length(x6));
% C6(:) = name(6);
% C=[C1 C2 C3 C4 C5 C6];
% group = C; %// set the groups according to the data above
% [p1] = anovan(data, group,'model','interaction') % ,'off' // Use the 'off' option to prevent the table/box plot from showing up.
% hold on
% plot([mean(y)], 'dg')
% title('Compare Random Data from Different Distributions')
% hold off

%%Only Ctot & Dtot
y2=y(:,5:6);
asami2={'CTot','DTot'};
[pp,tbl,stats]=anova1(y2,asami2);
figure
boxplot(y2,asami2);
hold on
plot(mean(y2), 'dg')
hold off
figure
% [c,~,~,gnames] = multcompare(stats);
[c, m, h, nms] = multcompare(stats); %(stats,'alpha',.05,'ctype','bonferroni');
title('Multi comparison-Post HOC test')
% xlabel('')
% https://uk.mathworks.com/help/stats/multcompare.html

%%%% For Men
[ppmen,tblmen,statsmen]=anova1(men,asami1);
figure
boxplot(men,asami1);
hold on
plot(mean(men), 'dg')
hold off
figure
[cmen, mmen, hmen, nms] = multcompare(statsmen); %(stats,'alpha',.05,'ctype','bonferroni');,'ctype','lsd'
title('Multi comparison-Post HOC test')
%%%% For Women
[ppwomen,tblwomen,statswomen]=anova1(Women,asami1);
figure
boxplot(Women,asami1);
hold on
plot(mean(Women), 'dg')
hold off
figure
[cwomen, mwomen, hwomen, nms] = multcompare(statswomen,'ctype','lsd'); %(stats,'alpha',.05,'ctype','bonferroni'); ,'ctype','lsd'
title('Multi comparison-Post HOC test')
%%%% For Ctot
y2=[Women(:,5);men(:,5)]';
[ppwomen,tblwomen,statswomen]=anova1(Women,asami1);
figure
boxplot(Women,asami1);
hold on
plot(mean(Women), 'dg')
hold off
figure
[cwomen, mwomen, hwomen, nms] = multcompare(statswomen,'ctype','lsd'); %(stats,'alpha',.05,'ctype','bonferroni'); ,'ctype','lsd'
title('Multi comparison-Post HOC test')
%%%% For Dtot

%% T student test
[h,pca,ci,stats] = ttest(DrNamazi.Ctot(:),DrNamazi.DTotal(:));
%pca=9.7375e-12,    tstat: -7.6952 , df: 101, sd: 25.1036

[h,p,ci,stats] = ttest2(DrNamazi.Ctot(:),DrNamazi.DTotal(:)); %(x,y,'Vartype','unequal')
% p = 1.6404e-18, tstat: -9.7017, df: 202 , sd: 14.0797
[h,p,ci,stats] = ttest2(DrNamazi.Ctot(:),DrNamazi.DTotal(:),'Vartype','unequal'); %(x,y,'Vartype','unequal')
% p = 1.7800e-18, tstat: -9.7017, df: 199.5561 , sd: [2×1 double]
%% Wilcoxon rank sum test - equal medians - different lengths -continuous distributions
[p,h]=ranksum(DrNamazi.Ctot(:),DrNamazi.DTotal(:))
% h=1 shows that they are significantly different
%% Kruskal-Wallis test
%Test the null hypothesis that the sample data from each column in x comes from the same distribution.
p = kruskalwallis(y,asami1);
%The returned value of p indicates that kruskalwallis rejects the null hypothesis that all columns data 
%samples come from the same distribution at a 1% significance level. The ANOVA table provides additional
%test results, and the box plot visually presents the summary statistics for each column in x.
% For Salamati
y1=y(:,1:4);
[p,tbl,stats] = kruskalwallis(y1,asami1(1:4)); % ,'Alpha',0.95
c = multcompare(stats)
%% Histogram with a distribution fit
histfit(y(:,5),10,'kernel') % 'beta','kernel'
histfit(DrNamazi.Ctot(:),10,'normal') % 'beta','kernel'
%% Multiple-sample tests for equal variances
[p,stats] = vartestn(Pop);
%% boxcox transforms nonnormally distributed to normal distribution
% data to a set of data that has approximately normal distribution
% The Box-Cox transformation is a family of power transformations.
[transdat,lambda] = boxcox(Women(:,5));
figure
histfit(transdat,10,'normal') 
figure
plot(Women(:,5))
hold on 
plot(transdat)
corr(Women(:,5),transdat)
figure
histogram(transdat)
figure
histogram(Women(:,5))
% For Pop
for i=1:size(Pop,2)
    a=Pop(:,i);
    a(find(a == 0))=eps;
    [transdat_pop,lambda_pop] = boxcox(a);
    transpop(:,i)=transdat_pop;
    Lambdamat(:,i)=lambda_pop;  
    clear transdat_pop lambda_pop a
end
% Anova
asami1={'CA','CB','CC','CD','CTot','DTot'};
[pp_tn,tb_tn,stats_Pop_Norm]=anova1(transpop,asami1);
figure
boxplot(transpop,asami1);
hold on
plot(mean(transpop), 'dg')
hold off
figure
% [c,~,~,gnames] = multcompare(stats);
[c, m, h, gnames] = multcompare(stats_Pop_Norm);
title('Multi comparison-Post HOC test')

%% New Structure - For Pop

% Sign test - Instead Pair T test
% the hypothesis that data in x has a continuous distribution with zero median 
% Note that a hypothesis of zero median for x – y is not equivalent to a hypothesis of equal median for x and y.
[p,h,stats] = signtest(Pop(:,5),Pop(:,6))
% h = 1 indicates that signtest reject to the null hypothesis of zero median

[p,h,stats]  = signrank(Pop(:,5),Pop(:,6))
% the null hypothesis that x – y comes from a distribution with zero median.
% The test assumes that the data in x come from a continuous distribution symmetric about its median.

% kruskalwallis test - Instead T-test & Anova
%Test the null hypothesis that the sample data from each column in x comes from the same distribution.
[p_k_pop,tbl_k_pop,stats_k_pop] = kruskalwallis(Pop,asami1); % ,'Alpha',0.95
%The returned value of p indicates that kruskalwallis rejects the null hypothesis that all columns data 
%samples come from the same distribution at a 1% significance level. The ANOVA table provides additional
%test results, and the box plot visually presents the summary statistics for each column in x.
% For Salamati
[c, m, h, gnames] = multcompare(stats_k_pop);
title('Multi comparison-Post HOC test')
% Only for GHQ questionnaire
[p_k_pop,tbl_k_pop,stats_k_pop1_4] = kruskalwallis(Pop(1:4),asami1(1:4)); % ,'Alpha',0.95
[c, m, h, gnames] = multcompare(stats_k_pop1_4);
title('Multi comparison-Post HOC test')

% Wilcoxon , Mann-Whithney = Ranksum - Instead T test
[p_r_pop,h_r_pop,stats_r_pop]=ranksum(Pop(:,5),Pop(:,6))

% Levene test
x1=Pop(:,1);
x1(:,2)=ones(length(x1),1);
x2=Pop(:,2);
x2(:,2)=2*ones(length(x2),1);
x3=Pop(:,3);
x3(:,2)=3*ones(length(x3),1);
x4=Pop(:,4);
x4(:,2)=4*ones(length(x4),1);
x5=Pop(:,5);
x5(:,2)=5*ones(length(x5),1);
x6=Pop(:,6);
x6(:,2)=6*ones(length(x6),1);
X=[x1;x2;x3;x4;x5;x6];
Levenetest(X)
clear x1 x2 x3 x4 x5 x6 X
% Levene's Test for Equality of Variances F=63.6752, df1= 5, df2=606
% Probability associated to the F statistic = 0.0000
% The associated probability for the F test is smaller than 0.05
% So, the assumption of homoscedasticity was not met.

%%%% Only for Ctot & Dtot
x5=Pop(:,5);
x5(:,2)=1*ones(length(x5),1);
x6=Pop(:,6);
x6(:,2)=2*ones(length(x6),1);
X=[x5;x6];
Levenetest(X)
clear x5 x6 X
% Levene's Test for Equality of Variances F=2.1118, df1= 1, df2=202
% Probability associated to the F statistic = 0.1477
% The associated probability for the F test is equal or larger than 0.05
% So, the assumption of homoscedasticity was met.

%% New Structure - For Men

% Sign test - Instead Pair T test
[p_men_s,h_men_s,stats_men_s] = signtest(men(:,5),men(:,6))
% h = 1 indicates that signtest reject to the null hypothesis of zero median

[p_men_sr,h_men_sr,stats_men_sr]  = signrank(men(:,5),men(:,6))
% the null hypothesis that x – y comes from a distribution with zero median.
% The test assumes that the data in x come from a continuous distribution symmetric about its median.

% kruskalwallis test - Instead T-test & Anova
%Test the null hypothesis that the sample data from each column in x comes from the same distribution.
[p_k_men,tbl_k_men,stats_k_men] = kruskalwallis(men,asami1); % ,'Alpha',0.95
%The returned value of p indicates that kruskalwallis rejects the null hypothesis that all columns data 
%samples come from the same distribution at a 1% significance level. The ANOVA table provides additional
%test results, and the box plot visually presents the summary statistics for each column in x.
% For Salamati
figure
[c, m, h, gnames] = multcompare(stats_k_men);
title('Multi comparison-Post HOC test')
% Only for GHQ questionnaire
[p_k_men,tbl_k_men,stats_k_men1_4] = kruskalwallis(men(1:4),asami1(1:4)); % ,'Alpha',0.95
[c, m, h, gnames] = multcompare(stats_k_men1_4);
title('Multi comparison-Post HOC test')

% Wilcoxon , Mann-Whithney = Ranksum - Instead T test
[p_r_men,h_r_men,stats_r_men]=ranksum(men(:,5),men(:,6))

% Levene test
x1=men(:,1);
x1(:,2)=ones(length(x1),1);
x2=men(:,2);
x2(:,2)=2*ones(length(x2),1);
x3=men(:,3);
x3(:,2)=3*ones(length(x3),1);
x4=men(:,4);
x4(:,2)=4*ones(length(x4),1);
x5=men(:,5);
x5(:,2)=5*ones(length(x5),1);
x6=men(:,6);
x6(:,2)=6*ones(length(x6),1);
X=[x1;x2;x3;x4;x5;x6];
Levenetest(X)
clear x1 x2 x3 x4 x5 x6 X
% Levene's Test for Equality of Variances F=33.9834, df1= 5, df2=258
% Probability associated to the F statistic = 0.0000
% The associated probability for the F test is smaller than 0.05
% So, the assumption of homoscedasticity was not met.

%%%% Only for Ctot & Dtot
x5=men(:,5);
x5(:,2)=1*ones(length(x5),1);
x6=men(:,6);
x6(:,2)=2*ones(length(x6),1);
X=[x5;x6];
Levenetest(X)
clear x5 x6 X
% Levene's Test for Equality of Variances F=3.9405, df1= 1, df2=86
% Probability associated to the F statistic = 0.0503
% The associated probability for the F test is equal or larger than 0.05
% So, the assumption of homoscedasticity was met.

%% New Structure - For Women

% Sign test - Instead Pair T test
[p_Women_s,h_Women_s,stats_Women_s] = signtest(Women(:,5),Women(:,6))
% h = 1 indicates that signtest reject to the null hypothesis of zero median

[p_Women_sr,h_Women_sr,stats_Women_sr]  = signrank(Women(:,5),Women(:,6))
% the null hypothesis that x – y comes from a distribution with zero median.
% The test assumes that the data in x come from a continuous distribution symmetric about its median.

% kruskalwallis test - Instead T-test & Anova
%Test the null hypothesis that the sample data from each column in x comes from the same distribution.
[p_k_Women,tbl_k_Women,stats_k_Women] = kruskalwallis(Women,asami1); % ,'Alpha',0.95
%The returned value of p indicates that kruskalwallis rejects the null hypothesis that all columns data 
%samples come from the same distribution at a 1% significance level. The ANOVA table provides additional
%test results, and the box plot visually presents the summary statistics for each column in x.
% For Salamati
figure
[c, m, h, gnames] = multcompare(stats_k_Women);
title('Multi comparison-Post HOC test')
% Only for GHQ questionnaire
[p_k_Women,tbl_k_Women,stats_k_Women1_4] = kruskalwallis(Women(1:4),asami1(1:4)); % ,'Alpha',0.95
figure
[c, m, h, gnames] = multcompare(stats_k_Women1_4);
title('Multi comparison-Post HOC test')

% Wilcoxon , Mann-Whithney = Ranksum - Instead T test
[p_r_Women,h_r_Women,stats_r_Women]=ranksum(Women(:,5),Women(:,6))

% Levene test
x11=Women(:,1);
x11(:,2)=ones(length(x11),1);
x21=Women(:,2);
x21(:,2)=2*ones(length(x21),1);
x31=Women(:,3);
x31(:,2)=3*ones(length(x31),1);
x41=Women(:,4);
x41(:,2)=4*ones(length(x41),1);
x51=Women(:,5);
x51(:,2)=5*ones(length(x51),1);
x61=Women(:,6);
x61(:,2)=6*ones(length(x61),1);
X=[x11;x21;x31;x41;x51;x61];
[levenetest] = Levenetest(X)
clear x11 x21 x31 x41 x51 x61 X
% Levene's Test for Equality of Variances F=33.5898, df1= 5, df2=342
% Probability associated to the F statistic = 0.0000
% The associated probability for the F test is smaller than 0.05
% So, the assumption of homoscedasticity was not met.

%%%% Only for Ctot & Dtot
x51=Women(:,5);
x51(:,2)=1*ones(length(x51),1);
x61=Women(:,6);
x61(:,2)=2*ones(length(x61),1);
X=[x51;x61];
Levenetest(X)
clear x51 x61 X
% Levene's Test for Equality of Variances F=0.0019, df1= 1, df2=114
% Probability associated to the F statistic = 0.9656
% The associated probability for the F test is equal or larger than 0.05
% So, the assumption of homoscedasticity was met.

%% Interactions between Men & Women
x1=men(:,1);
x2=men(:,2);
x3=men(:,3);
x4=men(:,4);
x5=men(:,5);
x6=men(:,6);
x11=Women(:,1);
x21=Women(:,2);
x31=Women(:,3);
x41=Women(:,4);
x51=Women(:,5);
x61=Women(:,6);
data = [x1' x2' x3' x4' x5' x6' x11' x21' x31' x41' x51' x61']; %// Create row vector with your data
name1={'CA_M','CB_M','CC_M','CD_M','CTot_M','DTot_M','CA_W','CB_W','CC_W','CD_W','CTot_W','DTot_W'};
C1 = cell(1, length(x1));
C1(:) = name1(1);
C2 = cell(1, length(x2));
C2(:) = name1(2);
C3 = cell(1, length(x3));
C3(:) = name1(3);
C4 = cell(1, length(x4));
C4(:) = name1(4);
C5 = cell(1, length(x5));
C5(:) = name1(5);
C6 = cell(1, length(x6));
C6(:) = name1(6);
C11 = cell(1, length(x11));
C11(:) = name1(7);
C21 = cell(1, length(x21));
C21(:) = name1(8);
C31 = cell(1, length(x31));
C31(:) = name1(9);
C41 = cell(1, length(x41));
C41(:) = name1(10);
C51 = cell(1, length(x51));
C51(:) = name1(11);
C61 = cell(1, length(x61));
C61(:) = name1(12);
C=[C1 C2 C3 C4 C5 C6 C11 C21 C31 C41 C51 C61];
group = C; %// set the groups according to the data above

[p,tbl,stats_MW] = kruskalwallis(data, group)
figure
[c, m, h, gnames] = multcompare(stats_MW);
title('Multi comparison-Post HOC test')

%%%%% Only Ctot & Dtot
x5=men(:,5);
x6=men(:,6);
x51=Women(:,5);
x61=Women(:,6);
data_1 = [x5' x6' x51' x61']; %// Create row vector with your data
name1={'CTot_M','DTot_M','CTot_W','DTot_W'};
C1 = cell(1, length(x5));
C1(:) = name1(1);
C2 = cell(1, length(x6));
C2(:) = name1(2);
C3 = cell(1, length(x51));
C3(:) = name1(3);
C4 = cell(1, length(x61));
C4(:) = name1(4);
C=[C1 C2 C3 C4];
group = C; %// set the groups according to the data above

[p,tbl,stats_MW_CD] = kruskalwallis(data_1, group)
figure
[c, m, h, gnames] = multcompare(stats_MW_CD);
title('Multi comparison-Post HOC test')

%%%%% Only CA, CB, CC, CD for men & Women
x1=men(:,1);
x2=men(:,2);
x3=men(:,3);
x4=men(:,4);
x11=Women(:,1);
x21=Women(:,2);
x31=Women(:,3);
x41=Women(:,4);
data_2 = [x1' x2' x3' x4' x11' x21' x31' x41']; %// Create row vector with your data
name1={'CA_M','CB_M','CC_M','CD_M','CA_W','CB_W','CC_W','CD_W'};
C1 = cell(1, length(x1));
C1(:) = name1(1);
C2 = cell(1, length(x2));
C2(:) = name1(2);
C3 = cell(1, length(x3));
C3(:) = name1(3);
C4 = cell(1, length(x4));
C4(:) = name1(4);
C5 = cell(1, length(x11));
C5(:) = name1(5);
C6 = cell(1, length(x21));
C6(:) = name1(6);
C7 = cell(1, length(x31));
C7(:) = name1(7);
C8 = cell(1, length(x41));
C8(:) = name1(8);
C=[C1 C2 C3 C4 C5 C6 C7 C8];
group = C; %// set the groups according to the data above

[p,tbl,stats_MW_C] = kruskalwallis(data_2, group)
[c, m, h, gnames] = multcompare(stats_MW_C);

%%%% 
[p_r_WM_CA,h_r_WM_CA,stats_r_WM_CA]=ranksum(men(:,1),Women(:,1))
[p_r_WM_CB,h_r_WM_CB,stats_r_WM_CB]=ranksum(men(:,2),Women(:,2))
[p_r_WM_CC,h_r_WM_CC,stats_r_WM_CC]=ranksum(men(:,3),Women(:,3))
[p_r_WM_CD,h_r_WM_CD,stats_r_WM_CD]=ranksum(men(:,4),Women(:,4))
[p_r_WM_C,h_r_WM_C,stats_r_WM_C]=ranksum(men(:,5),Women(:,5))
[p_r_WM_D,h_r_WM_D,stats_r_WM_D]=ranksum(men(:,6),Women(:,6))