%% Script that processes reaction time (RT) results
clear,clc,close all

%% Read raw data
% Condition: 1.SRT80-QOmni;  2.SRT80-DM;  3.SRT95-QOmni;  4.SRT95-DM
T = readtable('../Behavioural.xlsx');
T.Properties.VariableNames = {'ID','Order','Condition','RT','Intelligibility'};
NSubj = max(T.ID);

%% Descriptive stats (mean values)
fprintf(sprintf('%s','Mean reaction time (ms)\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.2f\t\t\t%1.2f\t\t\t%1.2f\t\t\t%1.2f\n'],...
        nanmean(T.RT((T.ID==Subj) & T.Condition==1)),...
        nanmean(T.RT((T.ID==Subj) & T.Condition==2)),...
        nanmean(T.RT((T.ID==Subj) & T.Condition==3)),...
        nanmean(T.RT((T.ID==Subj) & T.Condition==4))))
end
fprintf(sprintf('MEAN\t%1.2f\t\t\t%1.2f\t\t\t%1.2f\t\t\t%1.2f\n',...
    nanmean(T.RT(T.Condition==1)),nanmean(T.RT(T.Condition==2)),...
    nanmean(T.RT(T.Condition==3)),nanmean(T.RT(T.Condition==4))))


fprintf(sprintf('\n\n%s','Mean reaction time differences in ms QOmni - DMNR\n\n'))
fprintf('\t\tSRT-80\t\tSRT-95\n')
for Subj=1:NSubj
    fprintf(sprintf(['P' num2str(Subj) '\t\t%1.0f\t\t\t%1.0f\n'],...
        nanmean(T.RT((T.ID==Subj) & T.Condition==2)) - ...
        nanmean(T.RT((T.ID==Subj) & T.Condition==1)),...
        nanmean(T.RT((T.ID==Subj) & T.Condition==4)) - ...
        nanmean(T.RT((T.ID==Subj) & T.Condition==3))))
end

%% 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 = 10;    LW = 1;   MS = 3;   sep = 0.04;

subplot('Position',[0.05 0.12 0.93 0.85]),hold on, grid on, box on
for Subj=1:NSubj
    cond = T.ID==Subj; 
    plot(T.Condition(cond)+0.005*randn(sum(cond),1)+(Subj-10)*sep,T.RT(cond),...
    'o','MarkerFaceColor',cmap(Subj,:),'MarkerEdgeColor','w','MarkerSize',MS,'LineWidth',0.2)
end
bp=boxplot([T.RT(T.Condition==1),T.RT(T.Condition==2),T.RT(T.Condition==3),...
    T.RT(T.Condition==4),],'notch','off','whisker',1000);
set(bp,'LineWidth',LW,'Color','k')
set(gca,'xtick',1:4,'xticklabel',{'Q-Omni','DM','Q-Omni','DM'},'FontSize',FS-1)
set(gca,'ytick',500:1000:6000,'FontSize',FS-1)
axis([0.5 4.5 -1000 6700])
shift = 2750;
plot([1.1 1.9],[3300 3300]+shift,'k','LineWidth',LW)
plot([1.1 1.1],[3160 3300]+shift,'k','LineWidth',LW)
plot([1.9 1.9],[3160 3300]+shift,'k','LineWidth',LW)
text(1.5,6350,'p = 0.0038','FontSize',FS,'HorizontalAlignment','center')
plot([2.1 2.9]+1,[3300 3300]+shift,'k','LineWidth',LW)
plot([2.1 2.1]+1,[3160 3300]+shift,'k','LineWidth',LW)
plot([2.9 2.9]+1,[3160 3300]+shift,'k','LineWidth',LW)
text(3.5,6350,'p = 0.0059','FontSize',FS,'HorizontalAlignment','center')
text(1,50,'\beta_{Int}^{-1} = 1677 ms','FontSize',FS,'HorizontalAlignment','center')
text(2,50,'(\beta_{Int}+\beta_{DM})^{-1} = 1624 ms','FontSize',FS,'HorizontalAlignment','center')
text(3,50,'\beta_{Int}^{-1} = 1527 ms','FontSize',FS,'HorizontalAlignment','center')
text(4,50,'(\beta_{Int}+\beta_{DM})^{-1} = 1490 ms','FontSize',FS,'HorizontalAlignment','center')
text(1.5,-700,'SRT-80','FontSize',FS,'HorizontalAlignment','center','fontweight','bold')
text(3.5,-700,'SRT-95','FontSize',FS,'HorizontalAlignment','center','fontweight','bold')
ylabel('Reaction time (ms)','FontSize',FS)

orient tall
print('-dpng','-r600','Figure4')
print('-depsc','Figure4')
close all

%% Statistical analysis
switch 'SRT80'
    case 'SRT80', cond = T.Condition==1 | T.Condition==2;   % SRT-80 [QOmni vs DM]
    case 'SRT95', cond = T.Condition==3 | T.Condition==4;   % SRT-95 [QOmni vs DM]
end
RT = T.RT(cond);                        % Reaction time
Run = nominal(T.Order(cond));           % Run
Sbj = nominal(T.ID(cond));              % Subject
HA = nominal(T.Condition(cond));        % Hearing aid program
TBL = table(Sbj,HA,Run,RT,'VariableNames',{'Sbj','HA','Run','RT'});
GLME = fitglme(TBL,'RT ~ HA + Run + (1|Sbj)','Distribution','Gamma','Link',-1);
disp(GLME)