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

% NB2: Siccome inizia a diventare oneroso, non faccio il disp di tutto;
% guardate nel workspace e nella command window

%% 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
% La Y rimane la pp
i_pp = find(strcmp('pulse pressure', elsa_labels));

% Tutte tranne la Y (pp) e le pressioni
i_other_variables = find(~strcmp('pulse pressure', elsa_labels) ...
    & ~strcmp('systolic blood pressure', elsa_labels) ...
    & ~strcmp('diastolic blood pressure', elsa_labels));

pp = elsa_reduced(:, i_pp);
X_no_intercept = elsa_reduced(:, i_other_variables);

%% Correlazione tra le variabili

% \/ Lo commento perché potrebbe essere lento. Decommentatelo per vedere. \/

% figure
% corrplot(elsa_reduced, 'VarNames', elsa_labels) 

% /\ Lo commento perché potrebbe essere lento. Decommentatelo per vedere. /\

%% Regressione lineare 
% Scrivo la design matrix, ricordandomi di aggiungere l'intercetta alla
% fine
X = [X_no_intercept ones(length(pp), 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(121)
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(122)
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} ')'])

%% 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 (tranne le unità di misura)')
disp(beta_hat')
disp('SE dei beta stimati (tranne le unità di misura)')
disp(se_beta_hat')

%% 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 in %')
disp(cv_beta_hat')