clc
clear all
close all

%% Modello non regolarizzato

load dati_protesi

% Outcome
y = data(:,1);
% Numero di osservazioni
n = length(y); 

% Creo le variabili dummy associate alla variabile patologia
frattura = zeros(n,1);
frattura(find(data(:,8)==1)) = 1;
necrosi = zeros(n,1);
necrosi(find(data(:,8)==2)) = 1;

% Costruzione della matrice X del modello di regressione
ind = [2 4 5 6 7];
X = [ones(n,1) data(:,ind) frattura necrosi];

% Normalizzazione con metodo min-max scaling
Xnorm = ones(size(X));
for k = 2:size(X,2)
    Xnorm(:,k) = (X(:,k)-min(X(:,k)))/(max(X(:,k))-min(X(:,k)));
end
X = Xnorm; 

% Stima dei coefficienti di regressione senza regolarizzazione
beta_hat = inv(X'*X)*X'*y;
y_hat = X*beta_hat;
residui = y-y_hat;
SSE = sum(residui.^2);
m = size(X,2)-1; % Numero di variabili nel modello
SST = sum((y-mean(y)).^2);

% R2 del modello non regolarizzato
R2_c = 1-SSE/SST; 
% R2 adjusted del modello non regolarizzato
R2_adj = 1- (n-1)/(n-m-1)*SSE/SST;
%RMSE del modello non regolarizzato
RMSE = sqrt(mean(residui.^2));

%% Regressione Ridge 

% Valore del parametro di regolarizzazione
lambda = 100;
% Matrice Q
Q = eye(size(X,2));
Q(1,1) = 0;
% Stima dei coefficienti con regolarizzazione Ridge
beta_hat_ridge = inv(X'*X+lambda*Q)*X'*y; 
% Predizioni
y_hat = X*beta_hat;
% Residui
residui = y-y_hat;

% Calcolo di R2, R2 adjusted e RMSE per il modello Ridge
SSE = sum(residui.^2);
SST = sum((y-mean(y)).^2);
R2_ridge = 1-SSE/SST;
R2_adj_ridge = 1- (n-1)/(n-m-1)*SSE/SST;
RMSE_ridge = sqrt(mean(residui.^2));

% Confronto tra i coefficienti del modello non regolarizzato e quelli del
% modello Ridge
[beta_hat beta_hat_ridge]

%% Regressione LASSO

% Valore del parametro di regolarizzazione
lambda = 1;
% Stima dei coefficienti di regressione
[beta_lasso stats] = lasso(X(:,2:end),y,'Alpha',1,'Lambda',lambda,'Standardize',false);
% Coefficienti stimati
beta_hat_lasso =[stats.Intercept; beta_lasso];

% Cross-validation per la scelta del valore ottimo di lambda
% Vettore delle potenze
p = -10:3:10;
% Vettore dei valori di lambda testati
lambda_vec = 10.^p;
% Regressione LASSO con 10-fold cross-validation sulla griglia lambda_vec
[lasso_model stats] = lasso(X(:,2:end),y,'Alpha',1,'Lambda',lambda_vec,'Standardize',false,'CV',10);

%% Regressione elastic net
% Vettore dei possibili valori di alpha
alpha_vec = [0.1 0.4 0.7 1];

% 5-fold cross-validation per la stima dei valori ottimi di alpha e lambda
mean_MSE_CV = zeros(length(alpha_vec),length(lambda_vec));
for k = 1:length(alpha_vec)
    rng(1)
    [lasso_model stats] = lasso(X(:,2:end),y,'Alpha',alpha_vec(k),'Lambda',lambda_vec,'Standardize',false,'CV',5);
    mean_MSE_CV(k,:) = stats.MSE;
    mse(k) = stats.MSE(stats.IndexMinMSE);
    lambda(k) = stats.LambdaMinMSE;
end

% Scelta dei valori ottimi di alpha e lambda
[m, ind] = min(mse); 
alpha_opt = alpha_vec(ind)
lambda_opt = lambda(ind)

% Modello elastic net con i parametri ottimi
[beta_enet stats] = lasso(X(:,2:end),y,'Alpha',alpha_opt,'Lambda',lambda_opt,'Standardize',false,'CV',5)

% Coefficienti del modello
beta_hat_enet = [stats.Intercept; beta_enet];
% Predizioni
y_hat = X*beta_hat_enet;
% Residui
residui = y-y_hat;

% Calcolo di R2, R2 adjusted e RMSE per il modello elastic net
SSE = sum(residui.^2);
SST = sum((y-mean(y)).^2);
R2_enet = 1-SSE/SST;
R2_adj_enet = 1- (n-1)/(n-m-1)*SSE/SST;
RMSE_enet = sqrt(mean(residui.^2));

% Confronto tra i coefficienti del modello non regolarizzato e quello
% regolarizzato con elastic net
[beta_hat beta_hat_enet]