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

%% Parte 1
clc
clear all
close all

%% Caricamento dati
load IVGTT_subset

%% Creazione del sotto-dataset di interesse
% Ci sono vari modi per farlo; all'esame va bene "a occhio", qui lo faccio,
% per esercizio, in maniera programmatica
i_to_keep = endsWith(labels, '_basal') | endsWith(labels, '_ss');

% Tolgo dai dati
data = data(:, i_to_keep);

% Mi ricordo di togliere anche dai vettori delle etichette e delle unità
% di misura (IMPORTANTE! SPESSO FONTE DI ERRORI)
labels = labels(:, i_to_keep);
units = units(:, i_to_keep);

%% Individuazione e eventuale rimozione dati mancanti
n_nan = sum(isnan(data(:))); % Questo il numero
any_nans_with_sum = n_nan > 0;
any_nans = any(isnan(data(:))); % Altro modo

disp('Analisi dei dati mancanti per tutta la tabella')
disp(['  - Ci sono ' num2str(n_nan) ' dati mancanti (risultato op. logica: ' num2str(any_nans) ')'])
disp('       => Poiché non ci sono dati mancanti, non serve fare nulla.')

% Per casa: verificare se ce ne sono di infiniti o non fisiologici

%% Parte 2
% Verificare le ipotesi di gaussianità su ciascuna variabile e sulle
% differenze
% NB: in laboratorio, abbiamo visto solo sulle differenze; qui, per
% esercizio, lascio lo svolgimento completo nello stesso ciclo for.

% Prima di procedere, volendo usare un unico ciclo for "per tutto", creo le
% variabili differenza (Accrescimento dinamico della matrice con end+1)
% Anche qui: soluzione abbastanza elegante; all'esame non serve.
variables_to_differentiate = {'glucose', 'insulin', 'c_reactive_protein'};
units_to_differentiate = {'mg/dL', 'pmol/mL', 'mg/L'};
for variable_name_cell = variables_to_differentiate
    
    % Estraggo la stringa dal cell array unitario (guardarlo sulla command
    % window o sul workspace per capire cosa intendo)
    variable_name_string = variable_name_cell{1}; 
    
    % Aggiungo i suffissi
    ss = [variable_name_string '_ss']; % Concatenazione
    basal = [variable_name_string '_basal'];
    
    % Accrescimento dinamico della matrice
    data(:, end+1) = data(:, strcmp(ss, labels)) ...
        - data(:, strcmp(basal, labels));

    % Accrescimento dinamico dei cell array di etichette e unità di misura
    labels = [labels [variable_name_string '_diff']];
    units = [units units_to_differentiate];

end

% Creata la matrice, faccio lo stesso ciclo for del laboratorio 5 es 1
N_variables = length(labels);
figure
for i=1:N_variables

    %% Calcolo di statistiche descrittive per ogni variabile
    
    % Seleziono la colonna
    variable_curr = data(:, i);

    % Individuo e conto i dati mancanti
    i_nan_curr = isnan(variable_curr);
    n_nan_curr = sum(i_nan_curr);

    % Tolgo i dati mancanti (altrimenti la risposta è NaN, v. lab 3)
    variable_curr = variable_curr(~i_nan_curr);

    % Calcolo media, sd, mediana e IQR
    mean_curr = mean(variable_curr);
    std_curr = std(variable_curr);
    median_curr = median(variable_curr);
    iqr_curr = iqr(variable_curr);

    % Skeweness e curtosi
    skewness_curr = skewness(variable_curr);
    kurtosis_curr = kurtosis(variable_curr);

    % !NEW! Test di gaussianità (Lilliefors perché è quello che c'è in
    % MATLAB)
    [h_gauss, p_gauss, kstat_gauss, critval_gauss] = ... % ... = a capo riga
                               lillietest(variable_curr);

    % "Per abitudine" faccio il display sulla command window tutto insieme.
    % Si poteva fare man mano
    disp(' ')
    disp(['La variabile ' labels{i} ' (' units{i} ') ha:'])
    disp(['  - Numero di valori mancanti: ' num2str(n_nan_curr)])
    disp(['  - Media (senza NaN): ' num2str(mean_curr)])
    disp(['  - SD (senza NaN): ' num2str(std_curr)])
    disp(['  - Mediana (senza NaN): ' num2str(median_curr)])
    disp(['  - IQR (senza NaN): ' num2str(iqr_curr)])
    disp(['  - Skeweness (senza NaN): ' num2str(skewness_curr)])
    disp(['  - Curtosi (senza NaN): ' num2str(kurtosis_curr)])
    disp(['  - p-value del Lilliefors test: ' num2str(p_gauss)])

    %% Rappresentazione grafica delle distribuzioni
    % Istogrammi nella prima riga
    subplot(2, N_variables, i) 
    N_bins = 20; % Per l'istogramma
    [counts_curr, centres_curr] = hist(variable_curr, N_bins);
    bar(centres_curr, counts_curr)
    xlabel(units{i})
    ylabel('Frequenze assolute (count)')
    subtitle(labels{i})

     % !NEW! QQ plot nella seconda riga
     subplot(2, N_variables, N_variables+i)
     qqplot(variable_curr)
     subtitle(labels{i})
     % N.B.: Non complichiamo troppo le cose: le ordinate avrebbero le
     % unità di misura, le ascisse no (domanda per casa: perché?);
     % all'esame va bene se il QQ-plot ha le etichette di default (tranne
     % eventualmente il titolo)
    
end


%% Parte 3
% Differenze tra stato stazionario e basale per glucosio
[h_g, p_g, ci_g, stats_g] = ...
    ttest(data(:, strcmp('glucose_diff', labels)), 0);

% Mostro i risultati e faccio una deduzione
disp(' ') % cosmesi
disp('Glicemia allo stato stazionario vs. glicemia basale')
disp(' ') 
disp('Ipotesi di gaussianità')
disp('  - Istogramma della differenza unimodale e senza particolari asimmetrie')
disp('  - Test di gaussianità inconclusivo (H0 non rifiutata => non posso dire nulla)')
disp('  - QQ-plot: sommariamente giace sulla rette, possibile coda destra')
disp('  - Skewness e curtosi intorno, rispettivamente, a 0 e 3')
disp('     => Situazione un po'' borderline: decidiamo di assumere per buona la gaussianità')
disp(' ') 
disp('t-test appaiato a due campioni (i dati sono due istanti di tempo dello stesso individuo)')
disp('  - H0 era che la media delle differenze fosse nulla')
disp('  - Alfa era stato lasciato di default a 0.05.')
disp('  - Il test era a due code.')
disp(['  - Il p-value è risultato ' num2str(p_g)])
disp('     => Rifiuto H0: con livello di significatività 0.05, posso affermare che le due medie sono diverse')

%% Parte 4
% Differenze tra stato stazionario e basale per insulina
[p_i, h_i, stats_i] = ... % Ordine diverso! Controllare sempre l'help
    signtest(data(:, strcmp('insulin_diff', labels)));

% Mostro i risultati e faccio una deduzione
disp(' ') % cosmesi
disp('Insulina allo stato stazionario vs. insulina basale')
disp(' ') 
disp('Ipotesi di gaussianità')
disp('  - Istogramma della differenza unimodale, ma con una probabile asimmetria')
disp('  - Test di gaussianità rifiuta H0: con significatività 0.05 posso dire che la variabile non è normale')
disp('  - QQ-plot: vista la forte deviazione dalla retta, conferma la non gaussianità')
disp('  - Skewness e curtosi: confermano ancora la non gaussianità perché lontane, rispettivamente, da 0 e 3')
disp('     => Non possiamo assumere gaussianità della differenza')
disp(' ') 
disp('Per decidere per il sign test o il signed rank test, devo verificare la simmetria o meno della distribuzione della differenza')
disp('  - Avevamo già visto l''asimmetria dall''istogramma')
disp('     => Sign test della differenza contro 0.')
disp(' ') 
disp('Sign test appaiato a due campioni (i dati sono due istanti di tempo dello stesso individuo)')
disp('  - H0 era che la mediana delle differenze fosse nulla')
disp('  - Alfa era stato lasciato di default a 0.05.')
disp('  - Il test era a due code.')
disp(['  - Il p-value è risultato ' num2str(p_i)])
disp('     => Rifiuto H0: con livello di significatività 0.05, posso affermare che la mediana della differenze è diversa da 0 (da questo, posso dedurre che le due popolazioni abbiano distribuzioni diverse; altrimenti la mediana della differenza sarebbe uguale).')

% Analogamente per la proteina c-reattiva; lasciato in preparazione
% all'esame.