%% Script that processes Self-reported effort measures
clear,clc,close all

%% Read raw data
% Self-reported effort (Effort): 1. No effort, 2. Very little, 3. Little
% 4. Moderate, 5. Considerable, 6. Much effort, 7. Extreme effort
T = readtable('../Self_Reported.xlsx');
T.Properties.VariableNames = {'ID','Order','Condition','Effort'};
NSubj = max(T.ID);

%% Mean Self-reported effort per participant and condition
% Condition: 1.SRT80-QOmni;  2.SRT80-DM;  3.SRT95-QOmni;  4.SRT95-DM
SR = nan(NSubj,4);          % Self-report per subject (rows) and condition (columns)
for Subj=1:NSubj
    for Test=1:4
        cond = T.ID==Subj & T.Condition==Test;
        SR(Subj,Test) = nanmean(T.Effort(cond));
    end
end

%% Descriptive stats
fprintf(sprintf('%s','Self-reported effort (1.No effort -- 7.Extreme effort)\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.Effort(T.ID==Subj & T.Condition==1)),...
        nanmean(T.Effort(T.ID==Subj & T.Condition==2)),...
        nanmean(T.Effort(T.ID==Subj & T.Condition==3)),...
        nanmean(T.Effort(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.Effort(T.Condition==1)),nanmean(T.Effort(T.Condition==2)),...
    nanmean(T.Effort(T.Condition==3)),nanmean(T.Effort(T.Condition==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.1 0.12 0.88 0.85]),hold on, grid on, box on
for Subj=1:NSubj
    plot((1:2)+0.1*randn(1,2),SR(Subj,1:2),'-o','Color',cmap(Subj,:),...
        'MarkerFaceColor',cmap(Subj,:),'MarkerEdgeColor','w','MarkerSize',MS,'LineWidth',0.2);
    plot((3:4)+0.1*randn(1,2),SR(Subj,3:4),'-o','Color',cmap(Subj,:),...
        'MarkerFaceColor',cmap(Subj,:),'MarkerEdgeColor','w','MarkerSize',MS,'LineWidth',0.2);
end
bp=boxplot(SR,'notch','off','whisker',1000);
set(bp,'LineWidth',LW,'Color','k')
axis([0.5 4.5 0 9])
sep = 0.4;
plot([1.1 1.9],[8 8]-sep,'k','LineWidth',LW)
plot([1.1 1.1],[7.8 8]-sep,'k','LineWidth',LW)
plot([1.9 1.9],[7.8 8]-sep,'k','LineWidth',LW)
text(1.5,8.3,'p = 0.0024','FontSize',FS,'HorizontalAlignment','center')
plot([2.1 2.9]+1,[8 8]-sep,'k','LineWidth',LW)
plot([2.1 2.1]+1,[7.8 8]-sep,'k','LineWidth',LW)
plot([2.9 2.9]+1,[7.8 8]-sep,'k','LineWidth',LW)
text(3.5,8.3,'p = 0.2375','FontSize',FS,'HorizontalAlignment','center')
text(1,1.5,'e^{\beta_{Int}} = 4.81','FontSize',FS,'HorizontalAlignment','center')
text(2,1.5,'e^{[\beta_{Int}+\beta_{DM}]} = 4.28','FontSize',FS,'HorizontalAlignment','center')
text(3,1.5,'e^{\beta_{Int}} = 4.10','FontSize',FS,'HorizontalAlignment','center')
text(4,1.5,'e^{[\beta_{Int}+\beta_{DM}]} = 3.91','FontSize',FS,'HorizontalAlignment','center')
text(1.5,0.5,'SRT-80','FontSize',FS,'HorizontalAlignment','center','fontweight','bold')
text(3.5,0.5,'SRT-95','FontSize',FS,'HorizontalAlignment','center','fontweight','bold')
set(gca,'ytick',1:7,'yticklabel',{'No effort - 1','Very little - 2','Little - 3',...
    'Moderate - 4','Considerable - 5','Much - 6','Extreme effort - 7'},'FontSize',FS-1)
set(gca,'xtick',1:4,'xticklabel',{'Q-Omni','DM','Q-Omni','DM'},'FontSize',FS-1)
orient tall
print('-dpng','-r600','Figure5')
print('-depsc','Figure5')
close all

%% Statistical analysis
switch 'SRT80'
    case 'SRT80', cond = T.Condition==1 | T.Condition==2;   % SRT80-QOmni vs SRT80-DM
    case 'SRT95', cond = T.Condition==3 | T.Condition==4;   % SRT95-QOmni vs SRT95-DM
end
SR = T.Effort(cond);                    % Self-reported effort
Run = nominal(T.Order(cond));           % Run
Sbj = nominal(T.ID(cond));              % Subject
HA = nominal(T.Condition(cond));        % Hearing aid program (RES or SLN)
TBL = table(Sbj,HA,SR,Run,'VariableNames',{'Sbj','HA','SR','Run'});
GLME = fitglme(TBL,'SR ~ HA + Run + (1|Sbj)','Distribution','Poisson','Link','log');
disp(GLME)