%% Script that evaluates intelligibility scores
clear,clc,close all

%% Read raw data
T = readtable('../Behavioural.xlsx');
T.Properties.VariableNames = {'ID','Order','Condition','RT','Intelligibility'};
NSubj = max(T.ID);

%% Mean Intelligibility per subject and condition (WC)
% Condition: 1.SRT80-QOmni;  2.SRT80-DM;  3.SRT95-QOmni;  4.SRT95-DM
WC = nan(NSubj,4);     % Words correct per subject (rows) and condition (columns)
for Subj=1:NSubj
    for Test=1:4
        cond = T.ID==Subj & T.Condition==Test;
        WC(Subj,Test) = nanmean(T.Intelligibility(cond))*20;
    end
end

%% Descriptive stats
fprintf(sprintf('%s','Mean intelligibility (%%-correct)\n\n'))
fprintf('\t\tSRT-80 [Q-Omni]\tSRT-80 [DM]\t\tSRT-95 [Q-Omni]\tSRT-95 [DM]\n')
for Subj=1:NSubj
    fprintf(sprintf(['P' num2str(Subj) '\t\t%1.1f\t\t\t%1.1f\t\t\t%1.1f\t\t\t%1.1f\n'],...
        WC(Subj,1),WC(Subj,2),WC(Subj,3),WC(Subj,4)))
end
WC_Av = mean(WC);
fprintf(sprintf('MEAN\t%1.1f\t\t\t%1.1f\t\t\t%1.1f\t\t\t%1.1f\n',...
    WC_Av(1),WC_Av(2),WC_Av(3),WC_Av(4)))

%% Figure
figure('PaperSize',[20,8])
% Kelly's 22 colours of maximum constrast
cmap = ['#F2F3F4';'#222222';'#F3C300';'#875692';'#F38400';'#A1CAF1';...
    '#BE0032';'#C2B280';'#848482';'#008856';'#E68FAC';'#0067A5';...
    '#F99379';'#604E97';'#F6A600';'#B3446C';'#DCD300';'#882D17';...
    '#8DB600';'#654522';'#E25822';'#2B3D26'];
cmap = cmap([2:16 18:22],:);    % Selection of 20 adequate colours
FS = 12;    LW = 1;   MS = 6;
subplot('Position',[0.05 0.12 0.93 0.85]),hold on, grid on, box on
for Subj=1:NSubj
    plot((1:2)+0.10*randn(1,2),WC(Subj,1:2),'-o','Color',cmap(Subj,:),...
        'MarkerFaceColor',cmap(Subj,:),'MarkerEdgeColor','w','MarkerSize',MS,'LineWidth',0.2);
    plot((3:4)+0.10*randn(1,2),WC(Subj,3:4),'-o','Color',cmap(Subj,:),...
        'MarkerFaceColor',cmap(Subj,:),'MarkerEdgeColor','w','MarkerSize',MS,'LineWidth',0.2);
end
bp=boxplot(WC,'notch','off','whisker',1000);
set(bp,'LineWidth',LW,'Color','k')
axis([0.5 4.5 45 120])
shift = 5;
plot([1.1 1.9],[113 113]-shift,'k','LineWidth',LW)
plot([1.1 1.1],[110 113]-shift,'k','LineWidth',LW)
plot([1.9 1.9],[110 113]-shift,'k','LineWidth',LW)
text(1.5,113,'p = 9·10^{-52}','FontSize',FS,'HorizontalAlignment','center')
plot([2.1 2.9]+1,[113 113]-shift,'k','LineWidth',LW)
plot([2.1 2.1]+1,[110 113]-shift,'k','LineWidth',LW)
plot([2.9 2.9]+1,[110 113]-shift,'k','LineWidth',LW)
text(3.5,113,'p = 2·10^{-13}','FontSize',FS,'HorizontalAlignment','center')
text(1,59.5,'e^{\beta_{Int}} = 80.4%','FontSize',FS,'HorizontalAlignment','center')
text(2,59.5,'e^{[\beta_{Int}+\beta_{DM}]} = 85.2%','FontSize',FS,'HorizontalAlignment','center')
text(3,59.5,'e^{\beta_{Int}} = 89.6%','FontSize',FS,'HorizontalAlignment','center')
text(4,59.5,'e^{[\beta_{Int}+\beta_{DM}]} = 92.0%','FontSize',FS,'HorizontalAlignment','center')
text(1.5,50,'SRT-80','FontSize',FS,'HorizontalAlignment','center','fontweight','bold')
text(3.5,50,'SRT-95','FontSize',FS,'HorizontalAlignment','center','fontweight','bold')
set(gca,'xtick',1:4,'xticklabel',{'Q-Omni','DM','Q-Omni','DM'},'FontSize',FS)
set(gca,'ytick',0:20:100,'FontSize',FS)
ylabel('Intelligibility (%-correct)','FontSize',FS)
orient tall
print('-dpng','-r600','Figure3')
print('-depsc','Figure3')
close all

%% Statistical analysis
switch 'SRT80'
    case 'SRT80', cond = T.Condition==1 | T.Condition==2; % SRT80-QOmni or SRT80-DM
    case 'SRT95', cond = T.Condition==3 | T.Condition==4; % SRT95-QOmni or SRT95-DM
end
SiN = T.Intelligibility(cond)*20;           % Words correct (Intelligibility)
Run = nominal(T.Order(cond));               % Run
Sbj = nominal(T.ID(cond));                  % Subject
HA = nominal(T.Condition(cond));            % Hearing Aid program
TBL = table(Sbj,HA,SiN,Run,'VariableNames',{'Sbj','HA','SiN','Run'});
GLME = fitglme(TBL,'SiN ~ HA + Run + (1|Sbj)','Distribution','Poisson','Link','log');
disp(GLME)