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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% La prima parte è stata svolta nel laboratorio 7 %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Parte 1 
clc        
clear all  
close all  

%% Caricamento dati
load bodyfat

%% Definizione del problema di regressione lineare
% Cerco gli indici per costruire Y e X
i_y_mask = strcmp('BodyFat', labels);
i_X_mask = ~i_y_mask;

% Oppure "a occhio"
% i_y = 1
% i_X = 2:end % Sarebbero tutte le colonne, saltando la 1

Y = data(:, i_y_mask); 
X_no_intercept = data(:, i_X_mask);
X = [X_no_intercept ones(length(Y), 1)]; % Aggiungo l'intercetta

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

% Predizione del modello
Y_hat = X*beta_hat;

% Calcolo dei residui
residuals = Y - Y_hat;

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

% Intervallo di confidenza dei parametri
% Limite inferiore: due SE sotto il valore stimato
ci_lower_beta_hat = beta_hat - 1.96*se_beta_hat;
% Limite superiore: due SE sopra il valore stimato
ci_upper_beta_hat = beta_hat + 1.96*se_beta_hat;

% Coefficiente di variazione dei parametri in %
cv_beta_hat = 100 * se_beta_hat ./ abs(beta_hat); % Ricordarsi il valore assoluto!

% Il disp diventa complicato: fate prima a guardare sul workspace.
% Il mio consiglio è costruire una matrice di appoggio e ricordarvi che
% la prima colonna sono i parametri, la seconda il loro SE, la terza il CV;
% e che l'ultima riga è l'intercetta (in questo caso).
% Poi, fate doppio click sul nome della variabile nel workspace e vedete
% e/o copiate quello che vi interessa
regression_results = [beta_hat se_beta_hat ...
    ci_lower_beta_hat ci_upper_beta_hat ...
    cv_beta_hat];

% RMSE
% Ricordatevi: in questo caso (e in molti altri) MATLAB si scrive "quasi come si legge"
% cioè: ROOT = sqrt; MEAN = mean; SQUARED ERROR = residuals.^2
rmse = sqrt(mean(residuals.^2)); 

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

R2 = 1 - SSE/SST;
disp(' ')
disp(['RMSE = ' num2str(rmse) units{i_y_mask}]) % Unità di misura!
disp('NB: in questo caso il "%" è l''unità di misura della percentuale di grasso corporeo!')
disp(['R^2 = ' num2str(R2)])

%% Analisi dei residui
% Plot (v. istruzioni che iniziano per "subtitle")
figure
subplot(221) 
plot(Y, Y_hat, '.')
subtitle('Valore predetto contro valore vero')
xlabel(['Valore vero: ' labels{i_y_mask} ' (' units{i_y_mask} ')'])
ylabel(['Valore predetto: ' labels{i_y_mask} ' (' units{i_y_mask} ')'])

subplot(222) 
hold on
scatter(Y_hat, residuals, 'o') 
plot([min(Y_hat) max(Y_hat)], [0 0], '--k') % Retta che passa per i punti (minimo di Y_hat, 0) e (massimo di Y_hat, 0)
subtitle('Residui')
xlabel(['Valore predetto: ' labels{i_y_mask} ' (' units{i_y_mask} ')']) 
ylabel(['Residui: ' labels{i_y_mask} ' (' units{i_y_mask} ')'])

subplot(223)
N_bins = 10; % Per l'istogramma
[counts, centres] = hist(residuals, N_bins);
bar(centres, counts)
xlabel(units{i_y_mask})
ylabel('Frequenze assolute (count)')
subtitle('Istogramma dei residui')

subplot(224)
qqplot(residuals)
subtitle('QQ-plot dei residui')

% Verifica della media nulla
residuals_mean = mean(residuals);

% Skeweness e curtosi
residuals_skewness = skewness(residuals);
residuals_kurtosis = kurtosis(residuals);

% Test di gaussianità
[h_gauss, p_gauss, kstat_gauss, critval_gauss] = ... % ... = a capo riga
                           lillietest(residuals);

disp(' ')
disp('I residui hanno')
disp(['  - Media: ' num2str(residuals_mean)])
disp(['  - Skeweness: ' num2str(residuals_skewness)])
disp(['  - Curtosi: ' num2str(residuals_kurtosis)])
disp(['  - p-value del Lilliefors test: ' num2str(p_gauss)])

% Bianchezza dei residui con l'autocorrelazione
% Ordino i residui secondo il valore predetto
[Y_hat_sort, i_sort] = sort(Y_hat);

% Disegno l'autocorrelazione dei residui ordinati secondo il valore
% predetto
figure
autocorr(residuals(i_sort))
xlabel('Delta secondo il valore predetto')
title('Autocorrelazione dei residui')

% Riassunto delle conclusioni sui residui
disp('QQ-plot, istogramma, sknewness e curtosi sono compatibili con la gaussianità dei residui')
disp('Il test di lilliefors è inconclusivo (p value alto)')
disp('Residui evidentemente a media nulla')
disp('La funzione di autocorrelazione non evidenzia correlazioni inopportune => Bianchezza dei residui')
disp('Il plot dei residui contro il valore **predetto** è una banda/nuvola omogenea intorno alla retta orizzontale passante per 0, quindi possiamo considerarli a varianza costante')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Il laboratorio 8 parte da qui %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Confronto con fitlm
% Questa istruzione salva in "mdl" tantissime informazioni utili
mdl = fitlm(X, Y, 'Intercept', false);
% Oppure mdl = fitlm(X_no_intercept, Y, 'Intercept', true);

% Il disp fa "il riassunto"
% disp(mdl) % Lo commento per brevità

% Le informazioni sono dentro a mdl: l'help di MATLAB lista esattamente
% cosa si può trovare: https://it.mathworks.com/help/stats/linearmodel.html 
% Recuperiamone alcune che sappiamo calcolare a mano
beta_hat_fitlm = mdl.Coefficients.Estimate;
se_beta_hat_fitlm = sqrt(diag(mdl.CoefficientCovariance));
sigma2_hat_fitlm = mdl.MSE;
R2_fitlm = mdl.Rsquared.Ordinary;

% ATTENZIONE! L'RMSE non è proprio quello a cui siamo abituati: al
% denominatore ha i gradi di libertà n-m-1 (o n-numero_parametri). Noi
% continuiamo a usare "il nostro" con la media "semplice"
rmse_fitlm_diverso_da_quello_che_usiamo_noi = mdl.RMSE;

% Recuperiamo ora quelle "utili", cioè quelle per cui è meglio affidarsi a
% fitlm: i p value dei test statistici

% F test
p_F_test = coefTest(mdl);

disp(' ')
disp(['Il p-value dell''F-test era ' num2str(p_F_test)])
disp(['Siccome è minore del livello di significatività 0.05, possiamo dire' ...
    'che almeno uno dei coefficienti (oltre all''intercetta) è diverso da 0'])

% Test sui valori dei parametri
p_values = mdl.Coefficients.pValue;
disp(' ')
disp(['Ignorando, per ora, alcune sottigliezze interpretative, possiamo ' ...
    'dire che le variabili i cui p value corrispondenti sono <0.05 sono' ...
    ' significativamente legate alla variabile dipendente'])
disp('Segue lista') % Ignorate pure il codice qui sotto; serve solo per stampare
labels_and_intercept = [labels(i_X_mask) 'intercept'];
significant_betas = labels_and_intercept(p_values < 0.05);
disp(significant_betas)
disp('Se nella lista sopra appare il termine "intercept", un''interpretazione migliore è che l''intercetta del modello è significativamente diversa da zero')

%% Calcolo del numero di condizionamento
% Calcolo gli autovalori di X'*X
eigenvalues = eig(X'*X);
condition_number_kappa = max(eigenvalues) / min(eigenvalues); 

disp(' ')
disp(['Il numero di condizionamento è ' num2str(condition_number_kappa)])
disp('Poiché si tratta di un valore molto alto, è presente collinearità.')

%% Quantificazione e riduzione della collinearità tramite VIF
% NB: la soluzione "non bonus" è quella del ciclo for interno
% NB2: non preoccupatevi dell'altro ciclo; era solo per mostrarvi il
% procedimento "completo"
% NB3: il messaggio è che, in teoria, bisognerebbe calcolare il VIF per
% tutte le variabili, toglierne eventualmente una e poi ricalcolare tutto.
% NB4: all'esame potrebbe esserci il calcolo del VIF (ciclo for), 
% ma NON ci sarà *sicuramente* l'intero procedimento (ciclo while)

% Inizializzo un contenitore per i valori del VIF
VIF_trace = [];

% Inizializzo un vettore con gli indici delle variabili
indices_vif = 1:size(X_no_intercept, 2);

while true % Continuo a ciclare finché non do l'istruzione "break"

    VIF_trace_column = NaN(size(X_no_intercept, 2),1); % Uso il valore NaN per simboleggiare un VIF non calcolato

    % Scandisco tutte le variabili ancora in gioco
    for i_variable_curr=indices_vif % Scandisco gli indici delle variabili
    
        % La variabile corrente sarà la Y per il calcolo del "suo" VIF
        i_Y_vif = i_variable_curr; 
    
        % Le altre variabili saranno la X
        i_X_vif = indices_vif(indices_vif ~= i_Y_vif); % Controllate nel workspace

        % NB: Convincetevi, usando il workspace, di cosa sta facendo
        % l'istruzione; suggerimento: sta togliendo un indice da una
        % sequenza che parte da 1 e va al numero di variabili
        
        % Suddivido opportunamente e aggiungo l'intercetta
        Y_vif = X_no_intercept(:, i_Y_vif);
        X_vif = [X_no_intercept(:, i_X_vif) ones(size(X_no_intercept, 1), 1)];
    
        % Stimo i parametri
        beta_hat_vif = (X_vif'*X_vif)\X_vif'*Y_vif;
    
        % Calcolo R^2 e, di conseguenza, il VIF
        residuals_vif = Y_vif - X_vif*beta_hat_vif;
        SSE_vif = sum(residuals_vif.^2);
        SST_vif = sum((Y_vif - mean(Y_vif)).^2);

        R2_vif = 1 - SSE_vif/SST_vif;
        VIF = 1/(1 - R2_vif);

        % Inserisco il VIF al punto giusto
        VIF_trace_column(i_variable_curr) = VIF;
    
    end % for

    % Inserisco nella traccia completa
    VIF_trace = [VIF_trace VIF_trace_column];

    % Cerco il VIF massimo
    [max_vif, i_max_vif] = max(VIF_trace_column);

    if max_vif > 5 % Se c'è almeno una variabile correlata...
        indices_vif = indices_vif(indices_vif ~= i_max_vif); % ... la tolgo per la prossima iterazione
    else % Altrimenti...
        break % ... esco dal ciclo perché "ho finito"
    end
    

end % while

disp(' ')
disp('Ciascuna colonna della matrice qui sotto rappresenta una iterazione del calcolo dei VIF')
disp('Dove è scritto "NaN" si intende che la variabile era stata eliminata alla iterazione precedente')
disp(VIF_trace)


%% Confrontare le stime dei parametri e la loro variabilità dopo la rimozione della variabile
% Faccio tutto con fitlm; per casa, fatelo a mano. 
% NB: all'esame, potrebbe essere richiesto di farlo a mano.

% Ri-alleno il modello iniziale con tutte le variabili
mdl_before_vif = fitlm(X_no_intercept, Y, 'Intercept', true);

% Cerco gli indici delle variabili rimaste
i_variables_left_after_vif = find(VIF_trace(:, end) > 0);

% Alleno un nuovo modello con le sole variabili post-VIF
mdl_after_vif = fitlm(X_no_intercept(:, i_variables_left_after_vif), Y, 'Intercept', true);

% Mi preparo le etichette
labels_before_vif = labels(i_X_mask);
labels_after_vif = labels_before_vif(i_variables_left_after_vif);

disp(' ')
disp(['Con gli strumenti che abbiamo a disposizione, purtroppo, non possiamo' ...
    ' direttamente includere i nomi delle variabili nel disp di fitlm.'])
disp(['Si potrebbe fare abbastanza facilmente, ma non complichiamoci la' ...
    ' vita: faccio un disp dei nomi delle variabili, in ordine, sotto al' ...
    ' disp del summary del modello.'])
disp(mdl_before_vif)
disp(labels_before_vif)
disp(mdl_after_vif)
disp(labels_after_vif)

% Calcolo i CV con quanto mi viene restituito da fitlm
% (v. sopra per la corrispondenza tra calcoli a mano e campi dell'oggetto
% in uscita da fitlm)
cv_before_vif = 100 * sqrt(diag(mdl_before_vif.CoefficientCovariance)) ...
    ./ abs(mdl_before_vif.Coefficients.Estimate);
cv_after_vif = 100 * sqrt(diag(mdl_after_vif.CoefficientCovariance)) ...
    ./ abs(mdl_after_vif.Coefficients.Estimate);

% Faccio il disp "completo" di quanto trovato.
% NB: soluzione un po' verbosa, non necessaria all'esame; mi serve solo per
% poter commentare meglio il risultato ottenuto. A voi basta fare come ho
% fatto sopra, mostrando le label corrispondenti e guardando "a occhio".

% Recupero le stime dei parametri e le inserisco in vettori
% "corrispondenti"
beta_hat_before_vif = mdl_before_vif.Coefficients.Estimate;
beta_hat_after_vif = NaN(length(beta_hat_before_vif), 1);
beta_hat_after_vif([1; i_variables_left_after_vif+1]) = mdl_after_vif.Coefficients.Estimate;

% Faccio lo stesso per i CV
cv_after_vif_table = NaN(length(cv_before_vif), 1);
cv_after_vif_table([1; i_variables_left_after_vif+1]) = cv_after_vif;

% Creo il vettore di stringhe con i nomi delle variabili
var_names_table = ['Intercept', labels_before_vif]';

beta_table = strcat(var_names_table, {' --> beta_hat: '}, ...
   cellstr(num2str(beta_hat_before_vif)), ...
    {' vs. '}, cellstr(num2str(beta_hat_after_vif)));

cv_table = strcat(var_names_table, {' --> CV: '}, ...
     cellstr(num2str(cv_before_vif)), ...
    {'% vs. '}, cellstr(num2str(cv_after_vif_table)), '%');

disp(' ')
disp('Confronto pre vs. post esclusione delle variabili col VIF')
for i=1:length(var_names_table)
    disp(beta_table{i})
    disp(cv_table{i})
end

disp(' ')
disp('Come ci aspettavamo, il CV percentuale dei parametri stimati diminuisce')