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

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

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