%% plot roc curves

clear

close all

load('\multi-lstm-summary-clinical_avg.mat','mean_curve_all','mean_curve_seperate_all','auc_micro','auc_micro_std','prob_all','label_all');

time_step = 6;

time_range = 12:time_step:72;

intervals= linspace(0, 1, 100);

colormaps = [{'[128,128,128]/255'};{'[0, 213, 255]/255'};{'[149, 0, 255]/255'};{'[123, 255, 0]/255'};{'[1,0,0]'};{'[0,0,1]'}];

count = 1;

for time_points = 66%12:12:72

index = find(time_range==time_points);

auc_micro(index);

mean_curve = mean_curve_all{index};

mean_curve_seperate = mean_curve_seperate_all{index};

ySEM = std(mean_curve_seperate)/sqrt(size(mean_curve_seperate,1));                              % Compute ‘Standard Error Of The Mean’ Of All Experiments At Each Value Of ‘x’

CI95 = tinv([0.025 0.975], size(mean_curve_seperate,1)-1);                    % Calculate 95% Probability Intervals Of t-Distribution

yCI95 = bsxfun(@times, ySEM, CI95(:));              % Calculate 95% Confidence Intervals Of All Experiments At Each Value Of ‘x’

% x is a vector, matrix, or any numeric array of data. NaNs are ignored.

% p is the confidence level (ie, 95 for 95% CI)

% The output is 1x2 vector showing the [lower,upper] interval values.

CIFcn = @(x,p)prctile(x,abs([0,100]-(100-p)/2));

auc_boostrap = [];

num_boostrap = 1000;

for iboost = 1:num_boostrap

    [labels, idx] = datasample(label_all{index},length(label_all{index}));

    preds = prob_all{index}(idx);

    [X_tmp,Y_tmp,T_tmp,AUC_tmp] = perfcurve(labels-1,preds,false);

    auc_boostrap = [auc_boostrap; AUC_tmp];

end

p = 95; 

CI = CIFcn(auc_boostrap,p);

% figure,

eval(['[l',num2str(count),', p] = boundedline(intervals, mean_curve, ySEM,''cmap'',' colormaps{count}, ',''alpha'')']);% -b r 'cmap', [123, 255, 0]/255 [149, 0, 255]/255 [0, 213, 255]/255 [128,128,128]/255

% h1 = outlinebounds(l,p);

hold on

%     

%     plot(intervals, mean_curve,'-r','linewidth',2);

    axis square

    box on

     grid on

    title('Average ROC Curves at Different Time Intervals')

    xlabel('1-Specificity')

    ylabel('Sensitivity')

% hold off

temp = [num2str(time_points),' h ',num2str(round(auc_micro(index),2)),' (',num2str(round(CI(1),2)),',' num2str(round(CI(2),2)),')'];

legend_disp{count} = temp;

% eval(['legend(l',num2str(index),',','''',legend_disp,'''',')'])

% legend boxoff 

count = count+1;

end

legend([l1 l2 l3 l4 l5 l6],legend_disp)

legend boxoff 

plot(intervals, intervals, 'k--')

%%

num_pts =[];

for i = 1:11

num_pts = [num_pts; size(label_all{i},1)];

end

figure, bar(num_pts);

set(gca,'xticklabel',num2cell(12:6:72));

xlabel('Hours After Cardiac Arrest');

ylabel('Number of Patients')

%% calibration curves

close all

clear

load('D:\Research\Cardiac_arrest_EEG\Codes\ComaPrognosticanUsingEEG-master\multiscale-lstm\multi-lstm-summary-clinical_avg.mat','mean_curve_all','mean_curve_seperate_all','auc_micro','auc_micro_std','prob_all','label_all');

time_step = 6;

time_range = 12:time_step:72;

intervals= linspace(0, 1, 100);

colormaps = [{'[128,128,128]/255'};{'[0, 213, 255]/255'};{'[149, 0, 255]/255'};{'[123, 255, 0]/255'};{'[1,0,0]'};{'[0,0,1]'}];

count = 1;

for time_points = 12:12:72

index = find(time_range==time_points);

yp = prob_all{index};

y = label_all{index};

y(y==2)=0;

% x is a vector, matrix, or any numeric array of data. NaNs are ignored.

% p is the confidence level (ie, 95 for 95% CI)

% The output is 1x2 vector showing the [lower,upper] interval values.

CIFcn = @(x,p)prctile(x,abs([0,100]-(100-p)/2));

M = 10;

xx = linspace(0, 1, M+1); 

yy = [];     

for k = 1:M

    y_k = y(yp>=xx(k) & yp<xx(k+1));

    yy(k) = sum(y_k)/length(y_k);

end

xx = linspace(0, 1, M+1);

xx = (xx(1:M)+xx(2:end))/2;

Brier_mean(index) = nansum(abs(xx-yy))/sum(~isnan(yy)); %% Brier's Score

hold on

eval(['l', num2str(count),' = plot(xx, yy,''Color'',',colormaps{count}, ',''linewidth'',2);']);

xx = linspace(0, 1, M+1); 

Brier_boostrap = [];

num_boostrap = 1000;

for iboost = 1:num_boostrap

    [labels, idx] = datasample(label_all{index},length(label_all{index}));

    yp = prob_all{index}(idx);

    yy = [];

    xx = linspace(0, 1, M+1); 

for k = 1:M

    y_k = y(yp>=xx(k) & yp<xx(k+1));

    yy(k) = sum(y_k)/length(y_k);

end

xx = linspace(0, 1, M+1);

xx = (xx(1:M)+xx(2:end))/2;

Brier_boostrap(iboost) = nansum(abs(xx-yy))/sum(~isnan(yy));

end

p = 95; 

CI = CIFcn(Brier_boostrap,p);

    axis square

    box on

    grid on

    title('Calibration')

    xlabel('Prediction')

    ylabel('Proportion of the positive')

temp = [num2str(time_points),' h ',num2str(round(Brier_mean(index),2)),' (',num2str(round(CI(1),2)),',' num2str(round(CI(2),2)),')'];

legend_disp{count} = temp;

% eval(['legend(l',num2str(index),',','''',legend_disp,'''',')'])

% legend boxoff 

count = count+1;

end

xx = linspace(0, 1, M+1); 

plot(xx, xx, 'k--')

legend([l1 l2 l3 l4 l5 l6],legend_disp)

legend boxoff