%% LABORATORIO 10 - METODI STATISTICI PER LA BIOINGEGNERIA
%% Canale B - A.A. 2024/2025
%% Docente: Enrico Longato

clc
clear all
close all

%% Definizione del seed del generatore di numeri casuali
% Serve perché le parti "casuali" degli algoritmi implementati da MATLAB
% diano sempre lo stesso risultato
rng default
rng(42) 

%% Caricamento dati
load data_lab10

%% Parte 1 - Kaplan-Meier
[km_values, km_t] = ecdf(SurvTime, ...
                         'Censoring', ~C, ... % Attenzione che C va passata "al contrario", ovvero la function vuole 1 se il dato è censored
                         'Function', 'survivor'); % Questo serve per farsi dare i valori della S(t)

% Preparo un vettore per inserire le "crocette" che convenzionalmente si
% vedono sul plot della KM
censoring_times = SurvTime(C == 0);
S_of_t_at_censoring_times = interp1(km_t, km_values, censoring_times, 'previous');

% Plot (si poteva semplicemente copiare l'istruzione partendo da "ecdf"
figure
hold on
stairs(km_t, km_values, 'k') % KM
scatter(censoring_times, S_of_t_at_censoring_times, 'k+', 'LineWidth', 1.1) % Crocette
title('Survival function sull''intera popolazione (metodo KM)')
xlabel('Tempo (anni)')
ylabel('S(t)')
ylim([0 1])
xlim([0 20])

%% Parte 2 - Cox proportional hazards model "sul valore continuo del BMI"
% Fitto il modello di Cox
[b_cont, logl_cont, H_cont, stats_cont] = coxphfit(BMI, SurvTime, ...
    'censoring', ~C); % Attenzione sempre alla convenzione "strana"

% Calcolo l'hazard ratio
hr_b_cont = exp(b_cont);

disp(' ')
disp(['Considerando il BMI come variabile continua, ottengo un HR di ' num2str(hr_b_cont),  ...
    ' e relativo pvalue di ' num2str(stats_cont.p)])
disp('Concludo che aumentare di un''unità il valore del BMI aumenta l''hazard in maniera significativa')

%% Parte 3 - Kaplan-Meier per gruppi
% Divisione in due gruppi
i_group_1 = BMI < threshold;
i_group_2 = ~i_group_1;

% Si poteva fare un ciclo for, ma non vale la pena.
[km_values_1, km_t_1] = ecdf(SurvTime(i_group_1), ...
                         'Censoring', ~C(i_group_1), ... 
                         'Function', 'survivor'); 
[km_values_2, km_t_2] = ecdf(SurvTime(i_group_2), ...
                         'Censoring', ~C(i_group_2), ... 
                         'Function', 'survivor'); 

% Crocette
censoring_times_1 = SurvTime(i_group_1 & (C == 0)); % Notare l'AND logico
S_of_t_at_censoring_times_1 = interp1(km_t_1, km_values_1, censoring_times_1, 'previous');

censoring_times_2 = SurvTime(i_group_2 & (C == 0)); % Notare l'AND logico
S_of_t_at_censoring_times_2 = interp1(km_t_2, km_values_2, censoring_times_2, 'previous');

figure
hold on
stairs(km_t_1, km_values_1, 'r')
stairs(km_t_2, km_values_2, 'b')
scatter(censoring_times_1, S_of_t_at_censoring_times_1, 'r+', 'LineWidth', 1.1) 
scatter(censoring_times_2, S_of_t_at_censoring_times_2, 'b+', 'LineWidth', 1.1) 
title('Survival function sui due gruppi')
xlabel('Tempo (anni)')
ylabel('S(t)')
ylim([0 1])
xlim([0 20])
legend('Gruppo 1: BMI < soglia', 'Gruppo 2: BMI >= soglia')

%% Parte 4 - Cox proportional hazards model "sul gruppo"
% Fitto il modello di Cox
[b, logl, H, stats] = coxphfit(double(i_group_1), SurvTime, ...
    'censoring', ~C); % Attenzione sempre alla convenzione "strana"

% Calcolo l'hazard ratio
hr_b = exp(b);

disp(' ')
disp(['Coerentemente a quanto osservavo nel confronto tra le KM,' ...
     ' l''HR è significativamente diverso da zero (' num2str(stats.p) ').'])
disp(['In particolare, siccome l''HR è ' num2str(hr_b), ', quindi minore di 1, ' ...
    'posso ulteriormente dire che avere un BMI inferiore alla soglia è significativamente ' ...
    'associato ad un hazard inferiore, rispetto ad averlo sopra alla soglia.'])