%% 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

% NB: Analisi identiche/analoghe all'esercizio 1 che abbiamo visto in altri
% laboratori sono lasciate per casa senza soluzione per mantenere gli
% script focalizzati sul nuovo argomento.

%% 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;

%% 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

% 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

%% 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));
i_gluc = find(strcmp('glucose', elsa_labels));

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

%% Correlazione tra le tre variabili
figure
corrplot([pp, sbp, gluc], ... % Attenzione alla sintassi! help corrplot!
    'VarNames', elsa_labels([i_pp i_sbp i_gluc])) 
% 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'])
disp(['La variabile ' elsa_labels{i_gluc} ' non è correlata alle altre due'])

%% Regressione lineare di pp su SBP e glicemia
% Scrivo la design matrix, ricordandomi di aggiungere l'intercetta alla
% fine
X = [sbp gluc 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} ')'])



% NB: IL PIANO DI REGRESSIONE NON È PROGRAMMA D'ESAME. Serve solo per 
% farvi capire la differenza tra il grafico dei valori predetti e valori
% veri (che è sempre 2D) e il grafico della regressione (che, con più di 2
% variabili indipendenti diventa completamente in-visualizzabile)

% Plot della retta di regressione tra 70 e 220 mmHg sul terzo pannello
% Prima calcolo il piano
sbp_regr = 70:5:220; % Tra 70 e 220 con 5 di passo di campionamento
gluc_regr = 50:5:300; % Tra 50 e 300 con 5 di passo di campionamento

% Trasformo in matrici
[sbp_regr_mesh, gluc_regr_mesh] = meshgrid(sbp_regr, gluc_regr);

Y_regr = sbp_regr_mesh*beta_hat(1) ... % Calcolo "una variabile alla volta"
    + gluc_regr_mesh*beta_hat(2) ...   % perché sono matrici
    + beta_hat(3);

subplot(133)
hold on
plot3(sbp, gluc, Y, 'ok')
surf(sbp_regr, gluc_regr, Y_regr, ...
    'LineStyle', 'none')
view([-30 10]) % Cambia l'angolo da cui vediamo la superficie
subtitle('Piano di regressione')
zlabel([elsa_labels{i_pp} ' (' elsa_units{i_pp} ')'])
xlabel([elsa_labels{i_sbp} ' (' elsa_units{i_sbp} ')'])
ylabel([elsa_labels{i_gluc} ' (' elsa_units{i_gluc} ')'])


%% 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([' - ' elsa_labels{i_gluc} ': ' num2str(beta_hat(2)) ' (' num2str(se_beta_hat(2)) ') ' ...
    elsa_units{i_pp} '/' elsa_units{i_gluc}]) % Unità di misura di Y / unità di misura della variabile
disp([' - intercetta: ' num2str(beta_hat(3)) ' (' num2str(se_beta_hat(3)) ') ' ...
    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([' - ' elsa_labels{i_gluc} ': ' num2str(cv_beta_hat(2)) '%'])
disp([' - intercetta: ' num2str(cv_beta_hat(3)) '%']) % Unità di misura di Y