clc
clear all
close all

% Caricamento dati
load dati_protesi

% Outcome
y = data(:,1);

% Partiamo dal modello avente come variabili indipendenti X1-X5 (colonne 2-6 della matrice data)
n = length(y); 
X = [ones(n,1) data(:,2:6)];

%% Analisi di collinearità 

% Correlation plot
figure
corrplot(X(:,2:end),'VarNames',{'Altezza','Peso','Girovita','Lunghezza piede','Età'})
title('Correlation plot')
set(gca,'fontsize',12)

% Numero di condizionamento
lambda = eig(X'*X);
k = max(lambda)/min(lambda);

% Analisi VIF

% Variabili di partenza
Xred = X;
labels_red = {'altezza','peso','girovita','lunghezza piede','età'};

% Calcolo dell'indice VIF per ogni variabile nel modello
for j = 2:size(Xred,2)
    % Variabile j-esima considerata come outcome
    xj = Xred(:,j);
    % Matrice X delle variabili indipendenti usate per predire la variabile j-esima
    Xredj = Xred(:,[1:(j-1) (j+1):end]);
    % Stima dei coefficienti del modello
    beta_hat = inv(Xredj'*Xredj)*Xredj'*xj;
    % Calcolo di R2
    xj_hat = Xredj*beta_hat;
    SSE = sum((xj-xj_hat).^2);
    SST = sum((xj - mean(xj)).^2);
    R2j = 1 - SSE/SST;
    % Calcolo di VIF
    VIFj = 1/(1-R2j);
    disp(['VIF di ' labels_red{j-1} ': ' num2str(VIFj) ])
end
disp(' ')

% Rimuoviamo peso (X2, terza colonna di X)
Xred = X(:,[1:2 4:end]);
labels_red = {'altezza','girovita','lunghezza piede','età'};

% Ricalcolo gli indici VIF per tutte le variabili rimaste
for j = 2:size(Xred,2)
    xj = Xred(:,j);
    Xredj = Xred(:,[1:(j-1) (j+1):end]);
    beta_hat = inv(Xredj'*Xredj)*Xredj'*xj;
    xj_hat = Xredj*beta_hat;
    SSE = sum((xj-xj_hat).^2);
    SST = sum((xj - mean(xj)).^2);
    R2j = 1 - SSE/SST;
    VIFj = 1/(1-R2j);
    disp(['VIF di ' labels_red{j-1} ': ' num2str(VIFj) ])
end
% Sono tutti minori di 5 quindi non ci sono più multicollinearità importanti

%% Modello di regressione lineare senza la variabile peso
X = X(:,[1:2 4:end]);
% Stima dei coefficienti
beta_hat = inv(X'*X)*X'*y;
% Predizioni
y_hat = X*beta_hat;
% Residui
residui = y-y_hat;
% Somma dei residui al quadrato
SSE = sum(residui.^2);
% Stima a posteriori del sigma quadro
sigma2 = SSE/(n-size(X,2));
% Matrice di covarianza delle stime
covB = sigma2*inv(X'*X);
% Standard error
SE = sqrt(diag(covB));
% Coefficiente di variazione delle stime
CV = SE./abs(beta_hat)*100;
% Z-score dei coefficienti stimati
z_score = beta_hat./SE;

%% Aggiunta al modello delle variabili sesso e patologia

% Codifica della variabile patologia con 2 variabili dummy
frattura = zeros(n,1);
frattura(find(data(:,8)==1)) = 1;
necrosi = zeros(n,1);
necrosi(find(data(:,8)==2)) = 1;

% Costruzione della matrice X per il nuovo modello (senza peso, con sesso e patologia)
ind = [2 4 5 6 7];
X = [ones(n,1) data(:,ind) frattura necrosi];

% Stima dei coefficienti, standard error, CV, Z-score
beta_hat = inv(X'*X)*X'*y;
y_hat = X*beta_hat;
residui = y-y_hat;
SSE = sum(residui.^2);
sigma2 = SSE/(n-size(X,2));
covB = sigma2*inv(X'*X);
SE = sqrt(diag(covB));
CV = SE./abs(beta_hat)*100;
z_score = beta_hat./SE;
alpha = 0.05; % livello di significatività
zc = tinv(1-alpha/2,n-size(X,2)); % soglia critica

% Numero di parametri del modello completo
m_c = size(X,2);

% Calcolo di R2 e R2 adjusted
SST = sum((y-mean(y)).^2);
% R2 del modello completo
R2_c = 1-SSE/SST;
% R2 adjusted
R2_adj = 1- (n-1)/(n-(m_c-1)-1)*SSE/SST;

% Indici di parsimonia
AIC = n*log(SSE/n)+2*m_c;
BIC = n*log(SSE/n)+m_c*log(n);

disp(' ')
disp('Modello con variabile patologia: ')
disp([' R2 = ' num2str(R2_c)])
disp([' R2 adjusted = ' num2str(R2_adj)])
disp([' AIC = ' num2str(AIC)])
disp([' BIC = ' num2str(BIC)])

%% Rimozione della variabile patologia
X = [ones(n,1) data(:,ind)];

% Stima dei coefficienti, standard error, CV, Z-score
beta_hat = inv(X'*X)*X'*y;
y_hat = X*beta_hat;
residui = y-y_hat;
SSE = sum(residui.^2);
sigma2 = SSE/(n-size(X,2));
covB = sigma2*inv(X'*X);
SE = sqrt(diag(covB));
CV = SE./abs(beta_hat)*100;
z_score = beta_hat./SE;
alpha = 0.05; % livello di significatività
zc = tinv(1-alpha/2,n-size(X,2)); % soglia critica

% Numero di parametri del modello senza patologia
m_red1 = size(X,2);

% R2
R2_red1 = 1-SSE/SST;
% R2 adjusted
R2_adj = 1- (n-1)/(n-(m_red1-1)-1)*SSE/SST;
% Indici di parsimonia
AIC = n*log(SSE/n)+2*m_red1;
BIC = n*log(SSE/n)+m_red1*log(n);

disp(' ')
disp('Modello senza variabile patologia: ')
disp([' R2 = ' num2str(R2_red1)])
disp([' R2 adjusted = ' num2str(R2_adj)])
disp([' AIC = ' num2str(AIC)])
disp([' BIC = ' num2str(BIC)])

% Confronto tra i due modelli con l'F test
% Statistica del test
F = (R2_c - R2_red1)/(m_c-m_red1) / ((1-R2_c)/(n-m_c-1));
% Soglia critica
Fc = finv(1-alpha,m_c-m_red1,n-(m_c-1)-1);
% Decisione
disp('F test per il confronto del modello con/senza patologia:')
if F>Fc
    disp('  Rifiutiamo H0')
else
    disp('  Non possiamo rifiutare H0')
end

%% Rimozione anche della variabile lunghezza piede
ind =[2 4 6 7];
X = [ones(n,1) data(:,ind)];

% Stima dei coefficienti, standard error, CV, Z-score
beta_hat = inv(X'*X)*X'*y;
y_hat = X*beta_hat;
residui = y-y_hat;
SSE = sum(residui.^2);
sigma2 = SSE/(n-size(X,2));
covB = sigma2*inv(X'*X);
SE = sqrt(diag(covB));
CV = SE./abs(beta_hat)*100;
z_score = beta_hat./SE;
alpha = 0.05; % livello di significatività
zc = tinv(1-alpha/2,n-size(X,2)); % soglia critica

% Numero di parametri del modello senza patologia
m_red2 = size(X,2);

% R2
R2_red2 = 1-SSE/SST;
% R2 adjusted
R2_adj = 1- (n-1)/(n-(m_red2-1)-1)*SSE/SST;
% Indici di parsimonia
AIC = n*log(SSE/n)+2*m_red2;
BIC = n*log(SSE/n)+m_red2*log(n);

disp(' ')
disp('Modello senza variabile patologia: ')
disp([' R2 = ' num2str(R2_red2)])
disp([' R2 adjusted = ' num2str(R2_adj)])
disp([' AIC = ' num2str(AIC)])
disp([' BIC = ' num2str(BIC)])

% Confronto tra il modello senza patologia e quello senza patologia e lunghezza del piede mediante l'F test
% Statistica del test
F = (R2_red1 - R2_red2)/(m_red1-m_red2) / ((1-R2_red1)/(n-m_red1-1));
% Soglia critica
Fc = finv(1-alpha,m_red1-m_red2,n-(m_red1-1)-1);
% Decisione
disp('F test per il confronto del modello con/senza patologia:')
if F>Fc
    disp('  Rifiutiamo H0')
else
    disp('  Non possiamo rifiutare H0')
end