clc
clear all
close all

% Caricamento dati
load dati_protesi

% Outcome
y = data(:,1);

% Istogramma dei dati per la variabile di outcome
figure
histogram(y,'Normalization','probability')
title('Outcome Y')
xlabel('Diametro testa acetabolare [mm]')
ylabel('Frequenze relative')

n = length(y); 

% Stima dei parametri del modello di regressione lineare
% Variabili indipendenti X1-X5 (colonne 2-6 della matrice data)
X = [ones(n,1) data(:,2:6)];
beta_hat = inv(X'*X)*X'*y;

% Stima della varianza dell'errore
y_hat = X*beta_hat;
residui = y-y_hat;
SSE = sum(residui.^2);
sigma2 = SSE/(n-size(X,2));

% Plot outcome predetto vs outcome
figure
subplot(1,2,1)
plot(y,y_hat,'bo')
hold on
plot([min(y) max(y)],[min(y) max(y)],'r--','linewidth',2)
title('Predizioni')
xlabel('Diametro componente acetabolare reale [mm]')
ylabel('Diametro componente acetabolare predetto [mm]')
set(gca,'fontsize',14)

% Ordino i residui in base al valore predetto dell'outcome
[y_hat_ord,ind_ord]=sort(y_hat);
residui_ord = residui(ind_ord); 

% Plot dei residui ordinati vs. outcome predetta
subplot(1,2,2)
stem(y_hat_ord,residui_ord,'b-o')
hold on
plot([min(y_hat_ord) max(y_hat_ord)],[0 0],'r--','linewidth',2)
xlabel('Diametro componente acetabolare predetto [mm]')
ylabel('Residui [mm]')
title('Residui')
set(gca,'fontsize',14)

% Valutazione dell'errore di predizione del modello
RMSE = sqrt(SSE/n);

% Calcolo di R2
SST = sum((y-mean(y)).^2);
R2 = 1- SSE/SST;

% F test per il confronto del modello completo vs. modello nullo
m = size(X,2)-1;
F = ((SST-SSE)/m)/(SSE/(n-m-1)); % Statistica del test
alpha = 0.05; % Livello di significatività
Fc = finv(1-alpha,m,n-m-1); % Valore critico

disp("F test:")
if F>Fc
    disp("  Rifiutiamo H0")
else
    disp("  Non possiamo rifiutare H0")
end

% Analisi dei residui

% Check distribuzione normale

% Istogramma e Q-Q plot
figure
subplot(1,2,1)
histogram(residui,'normalization','probability')
title('Distribuzione dei residui')
xlabel('Residui [mm]')
ylabel('Frequenze relative')
set(gca,'fontsize',14)
subplot(1,2,2)
qqplot(residui)
title('Q-Q plot dei residui')
xlabel('Quantili della normale standard')
ylabel('Quantili dei residui')
set(gca,'fontsize',14)
box on

% Test di Lilliefors
[h, p] = lillietest(residui);
disp('Test di normalità dei residui: ')
disp(['  p-value = ' num2str(p)])

% Indici di forma campionari
gamma1 = skewness(residui);
gamma2 = kurtosis(residui);

disp(['Asimmetria: ' num2str(gamma1)])
disp(['Curtosi: ' num2str(gamma2)])

% Verifica della media nulla
media = mean(residui);
[h p] = ttest(residui);
disp(['Media dei residui: ' num2str(media) ' (p-value = ' num2str(p) ')'])

% Verifica dell'autocorrelazione dei residui
figure
subplot(1,2,1)
stem(y_hat_ord,residui_ord,'b-o')
hold on
plot([min(y_hat_ord) max(y_hat_ord)],[0 0],'r--','linewidth',2)
xlabel('Diametro componente acetabolare predetto [mm]')
ylabel('Residui [mm]')
title('Residui')
set(gca,'fontsize',14)
subplot(1,2,2)
autocorr(residui_ord,'NumSTD',2)
title('Funzione di autocorrelazione dei residui')
set(gca,'fontsize',14)

% Verifica varianza omogenea, outlier, trend rispetto al valore predetto
figure
plot(y_hat_ord,residui_ord,'bo')
xlabel('Diametro componente acetabolare predetto [mm]')
ylabel('Residui [mm]')
set(gca,'fontsize',14)

% Significatività statistica dei coefficienti della regressione

% Standard error
covB = sigma2*inv(X'*X);
SE = sqrt(diag(covB));
% Coefficiente di variazione
CV = SE./abs(beta_hat)*100;
% Intervallo di confidenza al 95%
CI = [beta_hat-1.96*SE beta_hat+1.96*SE]; 
% Possiamo anche mettere +/-2*SE per una stima più approssimata 

% Vettore degli Z-score
zscore = beta_hat./SE;
% Vettore dei p-value del test statistico sui coefficienti
p = 2*(1-tcdf(abs(zscore),n-m-1));
% Valore critico
zc = tinv(1-alpha/2,n-m-1);

disp('Test statistico sulla significatività dei coefficienti:')
disp(['  Valore critico: ' num2str(zc)])
zscore