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

%% Parte 1 
clc        
clear all  
close all  

%% Caricamento dati
load ELSAv2

%% Analisi dei dati mancanti/infiniti/non fisiologici in tutta la tabella
% Qui come maschera; indifferente se preferite usare gli indici col find
i_mask_nan = isnan(elsa);
i_mask_inf = isinf(elsa);
i_mask_neg = elsa < 0;

% Calcolo con sum (potevo fare length del vettore di indici col find)
n_nan = sum(i_mask_nan(:)); % oppure sum(sum(i_mask_nan)
n_inf = sum(i_mask_inf(:));
n_neg = sum(i_mask_neg(:));

disp(' ')
disp('Analisi su tutta la tabella')
disp(['  - Numero di valori mancanti: ' num2str(n_nan)])
disp(['  - Numero di valori infiniti: ' num2str(n_inf)])
disp(['  - Numero di valori negativi: ' num2str(n_neg)])

%% Rimozione dei soggetti con almeno un dato mancante/infinito/non fisiologico
% Procediamo sostituendo inf e valori non fisiologici con nan e poi
% togliendo i soggetti con nan
% (si poteva, altrimenti, togliere basandosi sull'OR delle condizioni
% logiche, ma, tipicamente, nella professione si procede come facciamo qui)
i_mask_nan_tot = i_mask_nan | i_mask_inf | i_mask_neg; % Tutti i valori da rendere NaN, inclusi quelli che già lo erano
elsa(i_mask_nan_tot) = nan; % Sostituisco

disp(' ')
disp('Dopo la sostituzione, correttamente:')
disp(['  - Numero di valori mancanti: ' num2str(sum(i_mask_nan_tot(:)))])

% Rimozione soggetti (v. laboratorio precedente)
elsa_reduced = elsa; % Teniamo sempre in memoria i dati originali quando la trasformazione potrebbe dover essere invertita

rows_with_at_least_a_nan = any(i_mask_nan_tot, 2); % Righe con almeno un nan

elsa_reduced = elsa_reduced(~rows_with_at_least_a_nan, :); % Tengo solo le righe senza nan

disp(' ')
disp('Dopo aver tolto i soggetti con almeno un dato mancante restano:')
disp(['  - Numero di righe: ' num2str(size(elsa_reduced, 1))])

%% Individuazione delle variabili di interesse
% Come al solito: qui uso il codice; all'esame va bene guardare a occhio
% l'indice di colonna
i_pp = find(strcmp('pulse pressure', elsa_labels));
i_sbp = find(strcmp('systolic blood pressure', elsa_labels));

pp = elsa_reduced(:, i_pp);
sbp = elsa_reduced(:, i_sbp);

%% Correlazione tra le due variabili
figure
corrplot([pp, sbp], ... % Attenzione alla sintassi! help corrplot!
    'VarNames', elsa_labels([i_pp i_sbp])) 
% NB: con questa funzione è più complicato aggiungere le unità di misura:
% per l'esame, non importa.

disp(' ')
disp(['Le variabili ' elsa_labels{i_pp} ' e ' elsa_labels{i_sbp} ' sono correlate'])

%% Verifica della gaussianità

% Skeweness e curtosi
skewness_pp = skewness(pp);
kurtosis_pp = kurtosis(pp);

skewness_sbp = skewness(sbp);
kurtosis_sbp = kurtosis(sbp);

% Test di gaussianità
[h_pp, p_pp, kstat_pp, critval_pp] = lillietest(pp);
[h_sbp, p_sbp, kstat_sbp, critval_sbp] = lillietest(sbp);

% QQ-plot
figure
subplot(121)
qqplot(pp)
subtitle(elsa_labels{i_pp})

subplot(122)
qqplot(sbp)
subtitle(elsa_labels{i_sbp})

disp(' ')
disp('Variabile heart rate')
disp([' - p-value test di gaussianità: ' num2str(p_pp)])
disp([' - Skewness: ' num2str(skewness_pp)])
disp([' - Curtosi: ' num2str(kurtosis_pp)])
disp('Sulla base dei plot e delle statistiche descrittive stampate sopra, possiamo dire che')
disp('  - Dall''istogramma, potrebbe esserci un''asimmetria. Proviamo a raccogliere ulteriori elementi.')
disp('  - Il test di gaussianità rifiuta l''ipotesi nulla. Può bastare così.')
disp('  - QQ-plot e skewness, comunque, confermano.')
disp('      => La variabile non è gaussiana.')

disp(' ')
disp('Variabile systolic blood pressure')
disp([' - p-value test di gaussianità: ' num2str(p_sbp)])
disp([' - Skewness: ' num2str(skewness_sbp)])
disp([' - Curtosi: ' num2str(kurtosis_sbp)])
disp('Sulla base dei plot e delle statistiche descrittive stampate sopra, possiamo dire che')
disp('  - Dall''istogramma non si rileva asimmetria grave.')
disp('  - Il test di gaussianità non ci permette di dire nulla.')
disp('  - QQ-plot sembra ragionevolmente su una retta.')
disp('  - Skewness e curtosi sono abbastanza in linea con i valori di 0 e 3.')
disp('      => Possiamo assumere la variabile come gaussiana.')

%% Regressione lineare di pp su SBP
% Scrivo la design matrix, ricordandomi di aggiungere l'intercetta alla
% fine
X = [sbp ones(length(sbp), 1)];

% Scrivo il vettore di outcome sincerandomi che sia una colonna 
% (size(Y) deve essere [N_righe, 1])
Y = pp;

% Utilizzo la formula risolutiva per trovare i parametri
beta_hat = (X'*X)\X'*Y;

%% Predizione del modello e calcolo dei residui
% Predizione
Y_hat = X*beta_hat;

% Calcolo dei residui
residuals = Y - Y_hat;

% RMSE: root mean squared error
% Notare come la sequenza di parole nel nome della metrica è anche l'ordine
% delle operazioni in MATLAB (tranne ovviamente il quadrato)
rmse = sqrt(mean(residuals.^2)); % "PUNTO elevamento a potenza" --> .^

disp(' ')
disp(['L''errore quadratico medio è: ' num2str(rmse) ' ' elsa_units{i_pp}])
% Mi raccomando le unità di misura

% Plot della predizione contro la realtà + plot dei residui
figure
subplot(131)
plot(Y, Y_hat, '.')
subtitle('Valore predetto contro valore vero')
xlabel(['Valore vero: ' elsa_labels{i_pp} ' (' elsa_units{i_pp} ')'])
ylabel(['Valore predetto: ' elsa_labels{i_pp} ' (' elsa_units{i_pp} ')'])

subplot(132)
stem(1:length(residuals), residuals, '.')
subtitle('Residui')
xlabel('Indice del soggetto') % O numero riga o qualunque cosa che mi faccia capire che avete capito
ylabel(['Residui: ' elsa_labels{i_pp} ' (' elsa_units{i_pp} ')'])

% Plot della retta di regressione tra 70 e 220 mmHg sul terzo pannello
% Prima calcolo la retta di regressione
sbp_regr = 70:.5:220; % Tra 70 e 220 con .5 di passo di campionamento
X_regr = [sbp_regr' ones(length(sbp_regr), 1)]; % Attenzione alle dimensioni
Y_regr = X_regr*beta_hat;

subplot(133)
hold on
plot(sbp, pp, 'o')
plot(sbp_regr, Y_regr, '--r')
subtitle('Retta di regressione')
xlabel([elsa_labels{i_pp} ' (' elsa_units{i_pp} ')'])
ylabel([elsa_labels{i_sbp} ' (' elsa_units{i_sbp} ')'])

%% Calcolare R^2
SSE = sum(residuals.^2);
SST = sum((Y - mean(Y)).^2);

R2 = 1 - SSE/SST;
disp(' ')
disp(['R^2 = ' num2str(R2)])

%% Standard error e CV dei parametri
% Standard error dei parametri
% 1. Calcolo inv(X'*X)
XtX = inv(X'*X);
% 2. Stimo la varianza a posteriori
sigma2_hat = SSE / (size(X, 1) - size(X, 2));
% 3. Moltiplico ed estraggo la diagonale
var_beta_hat = diag(sigma2_hat*XtX);
% 4. Faccio la radice
se_beta_hat = sqrt(var_beta_hat);

disp(' ')
disp('Valori stimati dei beta con relativi SE')
disp(' - variabile: beta_hat_i (SE) <-- Legenda')
disp([' - ' elsa_labels{i_sbp} ': ' num2str(beta_hat(1)) ' (' num2str(se_beta_hat(1)) ') ' ...
    elsa_units{i_pp} '/' elsa_units{i_sbp}]) % Unità di misura di Y / unità di misura della variabile
disp([' - intercetta: ' num2str(beta_hat(2)) ' (' num2str(se_beta_hat(2)) ') ' ...
    elsa_units{i_pp}]) % Unità di misura di Y

%% Coefficiente di variazione dei parametri in %
cv_beta_hat = 100 * se_beta_hat ./ abs(beta_hat); % Ricordarsi il valore assoluto!
disp(' ')
disp('Coefficienti di variazione')
disp(' - variabile: CV_hat_i <-- Legenda')
disp([' - ' elsa_labels{i_sbp} ': ' num2str(cv_beta_hat(1)) '%'])
disp([' - intercetta: ' num2str(cv_beta_hat(2)) '%']) % Unità di misura di Y