%% LABORATORIO 4 - 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_SS %% Pulizia dati: individuare NaN, valori infiniti, valori negativi % Qui come maschera; indifferente se preferite usare gli indici col find i_mask_nan = isnan(data); i_mask_inf = isinf(data); i_mask_neg = data < 0; %% Stampare il numero di dati di ciascuna categoria nell'intera tabella % Calcolo con sum (potevo fare length del vettore di indici col find) n_nan_tot = sum(i_mask_nan(:)); % oppure sum(sum(i_mask_nan) n_inf_tot = sum(i_mask_inf(:)); n_neg_tot = sum(i_mask_neg(:)); disp(' ') % Riga vuota per cosmesi disp('Analisi su tutta la tabella') disp([' - Numero di valori mancanti: ' num2str(n_nan_tot)]) disp([' - Numero di valori infiniti: ' num2str(n_inf_tot)]) disp([' - Numero di valori negativi: ' num2str(n_neg_tot)]) %% Pulizia dati: sostituire valori infiniti e valori negativi con NaN disp(' ') disp('Poiché non ci sono valori infiniti, né negativi, non faccio nulla.') % Risposta valida all'esame. %% Parte 2 % Uso un ciclo for "per esercizio"; all'esame vanno bene soluzioni anche % meno (o più) raffinate figure for i=1:length(labels) % Scandisco entrambe le variabili %% 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); % "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)]) %% Rappresentazione grafica delle distribuzioni % Boxplot nella prima riga della figura subplot(2, 2, i) boxplot(variable_curr) xlabel(labels{i}) ylabel(units{i}) % Istogrammi nella seconda riga (da cui 2 + ...) subplot(2, 2, 2 + 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}) end %% Parte 3 %% Mi concentro sulla variabile glucose_elderly e tolgo i NaN glucose_elderly = data(~isnan(data(:, 1)), 1); % Per casa: ricostruire la sequenza di operazioni operazioni, qui sopra % fatte in un'unica riga. %% t-test a un campione e a due code (1) % t-test per dimostrare che la media sia significativamente diversa da 100 target_mean = 100; % mg/dL % Se il secondo argomento è una costante, MATLAB capisce che il test è a un % campione solo. [h, p, ci, stats] = ... ttest(glucose_elderly, target_mean); % Mostro i risultati e faccio una deduzione disp(' ') % cosmesi disp(['H0 era che la media fosse ' num2str(target_mean) ' mg/dL.']) disp('Alfa era stato lasciato di default a 0.05.') disp('Il test era a due code.') disp(['Il p-value è risultato ' num2str(p)]) disp(['Rifiuto, dunque, H0 e dichiaro che la media della popolazione' ... ' da cui è estratto il campione è significativamente diversa da ' ... num2str(target_mean)]) %% t-test a un campione e a due code (2) % t-test per dimostrare che la media sia significativamente diversa da 82 target_mean_2 = 82; % mg/dL [h_2, p_2, ci_2, stats_2] = ... ttest(glucose_elderly, target_mean_2); % Mostro i risultati e faccio una deduzione disp(' ') % cosmesi disp(['H0 era che la media fosse ' num2str(target_mean_2) ' mg/dL.']) disp('Alfa era stato lasciato di default a 0.05.') disp('Il test era a due code.') disp(['Il p-value è risultato ' num2str(p_2)]) disp(['Non posso rifiutare H0, quindi non posso dedurre nulla sulla' ... ' media di popolazione.']) %% PARTE 4 %% t-test a un campione e a una coda % Uso la stessa media target del primo esempio (100 mg/dL), ma faccio il % test a una coda per capire se la media è superiore a 100 mg/dL. % Razionale: mi interessa individuare un forte diabete, cioè tale per cui % il glucosio a digiuno sia molto alto; non voglio "sprecare" potenza % statistica per l'altro caso. [h_one_tail, p_one_tail, ci_one_tail, stats_one_tail] = ... ttest(glucose_elderly, target_mean, 'Tail', 'right'); % Mostro i risultati e faccio una deduzione disp(' ') % cosmesi disp(['H0 era che la media fosse ' num2str(target_mean) ' mg/dL.']) disp('Alfa era stato lasciato di default a 0.05.') disp('Il test era a una coda (destra cioè l''ipotesi alternativa ha il >)') disp(['Il p-value è risultato ' num2str(p_one_tail)]) disp(['Non posso rifiutare H0, quindi non posso dedurre nulla sulla' ... ' media di popolazione.']) disp(['NOTA MOLTO BENE: a differenza del primo caso, in cui ' ... 'potevamo tranquillamente dire che la media fosse significativamente ' ... 'diversa da 100 mg/dL, qui non possiamo dire nulla perché, per protocollo, ' ... 'avendo usato un test a una coda destra, non ci interessava il caso in cui' ... ' la media fosse diversa, ma minore, di 100 mg/dL.']) disp(['NOTA ANCORA MEGLIO: questo non vuol dire che la media sia minore, ' ... 'perché non abbiamo fatto nessun calcolo che ci permettesse di dimostrarlo. ' ... 'Inoltre, *non vale* fare ora il test a una coda nell''altro senso.' ... ' Tutti i test vanno "pre-specificati" per protocollo']) %% PARTE 5 %% Esercizio sull'indimostrabilità dell'ipotesi nulla % Le due medie target sono target_mean_3 = 81.5; target_mean_4 = 82.5; % Test a due code come prima [h_3, p_3, ci_3, stats_3] = ... ttest(glucose_elderly, target_mean_3); [h_4, p_4, ci_4, stats_4] = ... ttest(glucose_elderly, target_mean_4); disp(' ') disp('Primo test') disp(['H0 era che la media fosse ' num2str(target_mean_3) ' mg/dL.']) disp('Alfa era stato lasciato di default a 0.05.') disp('Il test era a due code.') disp(['Il p-value è risultato ' num2str(p_3)]) disp(' ') disp('Secondo test') disp(['H0 era che la media fosse ' num2str(target_mean_4) ' mg/dL.']) disp('Alfa era stato lasciato di default a 0.05.') disp('Il test era a due code.') disp(['Il p-value è risultato ' num2str(p_4)]) disp(' ') disp(['Non posso rifiutare nessuna delle due ipotesi nulle e, quindi, ' ... 'non posso dire nulla in nessuno dei due casi.']) disp(' ') disp(['NOTA BENE: da questo esercizio è evidentissimo ' ... 'come l''ipotesi nulla non sia dimostrabile da un test' ... ' in cui non la possiamo rifiutare']) disp(['Infatti, se si potesse fare una cosa del genere, ci troveremmo ora' ... ' nel caso di aver dimostrato che la media è uguale a due numeri diversi, ' ... 'che è *chiaramente insensato*.']) %% PARTE 6 % Rifacciamo il primo test con livello di significatività 0.01 [h_001, p_001, ci_001, stats_001] = ... ttest(glucose_elderly, target_mean, 'Alpha', 0.01); % Stesse conclusioni di prima (omissis)