%% LABORATORIO METODI STATISTICI PER LA BIOINGEGNERIA - LABORATORIO 7

close all;
clearvars;
clc;

%% 1) Caricamento dati

% NOTA
%Non viene chiesta l'analisi per la ricerca di valori non utilizzabili (tipo NaN) 
%perchè sono già stati testati e normalizzati tramite zscore. I dati sono quindi pronti all'uso

load('data_lab7.mat');

%% 2) Si richiede la stima di Y utilizzando tutte le variabili presenti in X.
% Calcolare: - residui (res)
%            - residui al quadrato (rss), 
%            - la stima a posteriori dell'errore di misura (signa_hat_2), 
%            - stadarderror (SE)
%            - coefficiente di variazione (CV), il
%            - coefficiente di determinazione (R2), R2 aggiustato (R_adjusted).
% Mostrare: 
%            - scatterplot Y vs Yhat
%            - residui
%            - qq-plot per valutare il fit

% Calcolo i parametri beta con il modello 
beta_hat = (X'*X)\X'*Y;

% Calcolo la predizione del modello
y_hat = X*beta_hat;

% Calcolo i residui
res = Y-y_hat;

% Calcolo la somma dei residui al quadrato
rss = (res'*res);

% calcolo la stima a posteriori dell'errore di misura 
sigma_hat_2 = (res'*res)/(length(Y)-length(beta_hat));

% calcolo lo stadard error
SE =sqrt((sigma_hat_2)*diag(inv(X'*X)));

% Calcolo i coefficienti di variazione
CV = SE./abs(beta_hat)*100;

% Calcolo il coeff. di determinazione R2
[rho,p_val_R2] = corr(Y,y_hat);

R2=rho.^2;

% calcolo R2 aggiustato
m=length(beta_hat);    
n=length(Y);
R_adjusted=(R2-(m/(n-1)))*((n-1)/(n-m-1));

% Mostrare scatterplot, residui e qq-plot per valutare il fit del modello
figure(1)
subplot(311)
scatter(Y,y_hat)
xlabel('Dati (Y)')
ylabel('Predizione (Yhat)')
title('Scatter Plot')

subplot(312)
plot(res)
hold on
plot([1 length(res)],[0 0],'k--')
xlim([1 length(res)])
ylim([-1 1])
title('Residui')

subplot(313)
qqplot(res)

%% 3) E' presente multicollinearità tra le variabili. Valutare il Fattore di
% Inflazione della Varianza (VIF) per ogni variabile di X.

for ii=1:size(X,2)

    G=X;

    G(:,ii)=[]; % elimino la i-esima variabile

    % stimo la variabile i-esima stessa mediante tutte le altre
    beta_i=(G'*G)\G'*X(:,ii);

    % valuto la correlazione tra la i-esima variabile e la sua predizione
    % mediante le altre variabili
    R_X(ii)=(corr(X(:,ii),G*beta_i))^2;
    
    % VIF per la i-esima variabile
    VIF(ii,1)=1/(1-R_X(ii));

end

%% 4) selezionare le variabili da eliminare in base alla condizione
% VIF>10 e mantenere solo le variabili X_new da mantenere.
selected_VIF=VIF>10;

X_new=X;
X_new(:,selected_VIF)=[];

labels_new=labels(not(selected_VIF));

%% 5) Si ripeta il punto 1) utilizzando solo X_new
% Calcolo i parametri beta con il modello 
beta_hat_new = (X_new'*X_new)\X_new'*Y;

% Calcolo la predizione del modello
y_hat_new = X_new*beta_hat_new;

% Calcolo i residui
res_new = Y-y_hat_new;

% Calcolo la somma dei residui al quadrato
rss_new = (res_new'*res_new);

% calcolo la stima a posteriori dell'errore di misura 
sigma_hat_2_new = (res_new'*res_new)/(length(Y)-length(beta_hat_new));

% calcolo lo stadard error
SE_new =sqrt((sigma_hat_2_new)*diag(inv(X_new'*X_new)));

% Calcolo i coefficienti di variazione
CV_new = SE_new./abs(beta_hat_new)*100;

% Calcolo il coeff. di determinazione R2
[rho_new,p_val_R2_new] = corr(Y,y_hat_new);

R2_new=rho_new.^2;

% calcolo R2 aggiustato
m=length(beta_hat_new);    
n=length(Y);
R_adjusted_new=(R2_new-(m/(n-1)))*((n-1)/(n-m-1));

%% 6) Confrontare il modello ottenuto da X_new con quello ottenuto da X_mod
% (usare le ultime due variabili di X_new) tramite R2_adjusted calcolando AIC e BIC.
X_mod=X_new(:, [end-1 end]);

% Calcolo i parametri beta con il modello 
beta_hat_mod = (X_mod'*X_mod)\X_mod'*Y;

% Calcolo la predizione del modello
y_hat_mod = X_mod*beta_hat_mod;

% Calcolo i residui
res_mod = Y-y_hat_mod;

% Calcolo la somma dei residui al quadrato
rss_mod = (res_mod'*res_mod);

% calcolo la stima a posteriori dell'errore di misura 
sigma_hat_2_mod = (res_mod'*res_mod)/(length(Y)-length(beta_hat_mod));

% calcolo lo stadard error
SE_mod =sqrt((sigma_hat_2_mod)*diag(inv(X_mod'*X_mod)));

% Calcolo i coefficienti di variazione
CV_mod = SE_mod./abs(beta_hat_mod)*100;

% Calcolo il coeff. di determinazione R2
[rho_mod,p_val_R2_mod] = corr(Y,y_hat_mod);

R2_mod=rho_mod.^2;

% calcolo R2 aggiustato
m=length(beta_hat_mod);    
n=length(Y);
R_adjusted_mod=(R2_mod-(m/(n-1)))*((n-1)/(n-m-1));

disp('  ')
disp(['R adjusted X_new = ',num2str(R_adjusted_new)])
disp(['R adjusted X_mod = ',num2str(R_adjusted_mod)])
disp('  ')

AIC1= length(Y)*log((res_new'*res_new)/length(Y))+2*length(beta_hat_new);
AIC2= length(Y)*log((res_mod'*res_mod)/length(Y))+2*length(beta_hat_mod);

disp('  ')
disp(['AIC X_new = ',num2str(AIC1)])
disp(['AIC X_mod = ',num2str(AIC2)])
disp('  ')


BIC1= length(Y)*log((res_new'*res_new)/length(Y))+length(beta_hat_new)*log(length(Y));
BIC2=length(Y)*log((res_mod'*res_mod)/length(Y))+length(beta_hat_mod)*log(length(Y));
disp('  ')
disp(['BIC X_new = ',num2str(BIC1)])
disp(['BIC X_mod = ',num2str(BIC2)])
disp('  ')

%modello migliore = ?


%% Mostrare scatterplot, residui e qq-plot per valutare il fit del modello
figure(3)
subplot(321)
scatter(Y,y_hat_new)
xlim([-3 3])
ylim([-3 3])
xlabel('Dati (Y)')
ylabel('Predizione (YhatNew)')
title('Scatter Plot')

subplot(323)
plot(res_new)
hold on
plot([1 length(res_new)],[0 0],'k--')
xlim([1 length(res_new)])
ylim([-1 1])
title('Residui')

subplot(325)
qqplot(res_new)

subplot(322)
scatter(Y,y_hat_mod)
xlim([-3 3])
ylim([-3 3])
xlabel('Dati (Y)')
ylabel('Predizione (YhatMod)')
title('Scatter Plot')

subplot(324)
plot(res_mod)
hold on
plot([1 length(res_mod)],[0 0],'k--')
xlim([1 length(res_mod)])
ylim([-1 1])
title('Residui')

subplot(326)
qqplot(res_mod)


%% 7) (Opzionale) Analizzare la presenza di multicollinearità tra variabili tramite CN
% (Numero di condizionamento) e stampare sulla command window il valore di
% CN

lambda = eig(X'*X); % calcolo gli autovalori della matrice X'X

CN = max(lambda)/min(lambda);

if CN>1000
    disp(['Sono presenti multicollinearità tra le variabili. Il valore di CN è ' num2str(CN)]);
else
    disp(['Non sono presenti multicollinearità tra le variabili. Il valore di CN è ' num2str(CN)]);
end

%% 8) Esaminare la presenza di multicollinearità tra le variabili tramite
% Indice di Condizionamento (CI) mediante la funzione collintest.
% Attenzione: per visualizzare i nomi delle variabili, collinrtest richiede
% in input una table (array2table)
table_x=array2table(X, 'VariableNames', labels);
[sValue,condIdx,VarDecomp] = collintest(X);

%% 9) Mostrare i risultati di CI in una figure (vedi help di collintest).
% Settare TolIdx=10.
figure
collintest(table_x, Plot="on", TolIdx=10, TolProp=0.5);

% collintest uses TolIdx to decide which indices are large enough to infer a near dependency in the data. TolIdx is used only when the Plot argument is "on".
% Plot='on': 	
% collintest plots critical rows of the output VarDecomp, specifically, rows with condition indices above the input tolerance TolIdx.
% If a group of at least two variables in a critical row have variance-decomposition proportions above the input tolerance TolProp, the group is identified with red markers.
% collintest uses TolIdx to decide which indices are large enough to infer a near dependency in the data. TolIdx is used only when the Plot argument is "on".