%% COMPLEMENTO ALLE LEZIONI DI TEORIA - METODI STATISTICI PER LA BIOINGEGNERIA
%% Canale B - A.A. 2026/2027
%% Docente: Enrico Longato

% Semplice esercizio motivazionale sulla necessità della correzione per
% confronti multipli.

clc
clear all
close all

% Genero un campione "artificiale": so tutti i parametri, ma faccio finta
% di non saperli, ovvero che i dati che seguono mi siano stati dati e mi
% sia stata posta una domanda di interesse statistico.
rng(42) % Setto un seed per riproducibilità
data = randn(50, 1) + 10; % Media 10

figure % Velocemente disegno l'istogramma delle frequenze assolute
histogram(data) % Uso questa funzione per farvi vedere che va bene

% Adesso faccio una cosa "anomala": invece di fare un test statistico
% ragionevole, con un set di ipotesi ben chiaro in mente e prespecificato
% per protocollo, "sparo a caso" e provo una batteria di test statistici
% con medie target tutte diverse

medie_target_tutte_diverse = 9.5:0.05:10.5; % Provo tanti valori intorno a 10

p_value_corrispondenti = [];
for media_target_corrente=medie_target_tutte_diverse
    [h, p, ci, stats] = ttest(data, media_target_corrente);
    p_value_corrispondenti = [p_value_corrispondenti p]; % Concatenazione
end

% Adesso ordiniamo i p-value e controlliamo quanti sono <0.05
n_confronti = length(medie_target_tutte_diverse);
n_p_value_piccoli = sum(p_value_corrispondenti < 0.05);
disp(' ')
disp(['I p value minori di 0.05 sono ' num2str(n_p_value_piccoli) ...
    ' su ' num2str(n_confronti)])

% Proviamo ora a correggere per i confronti multipli.
% Correggerò mostrando in un grafico l'andamento dei p-value contro delle
% righe che corrispondono alle soglie delle varie correzioni.

% Ordino i pvalue in senso crescente
[p_ordinati, i_p] = sort(p_value_corrispondenti);

% Disegno la "riga" dei test non corretti (si dice il p value nominale)
x_nominale = [1 n_confronti];
y_nominale = [0.05 0.05]; % Sto solo disegnando la retta per punti

% Disegno la "riga" di Bonferroni
x_B = [1 n_confronti];
y_B = [0.05 0.05] / n_confronti; % Sto solo disegnando la retta per punti

% Disegno la spezzata (o scala) di Holm-Bonferroni
m = 1:n_confronti; % questo è "1, 2, ..., m-1, m" nella teoria
y_HB = 1 ./ fliplr(m) * 0.05; % fliplr(m) sarebbe "m, m-1, ecc..."

% Disegno la retta di Benjamini-Hochberg
y_BH = m * 0.05 / n_confronti;


figure
plot(m, p_ordinati, '.', 'MarkerSize', 15)
hold on
plot(x_nominale, y_nominale, 'LineWidth', 2)
plot(x_B, y_B, 'LineWidth', 2)
plot(m, y_HB, 'LineWidth', 2)
plot(m, y_BH, 'LineWidth', 2)
legend('p value', 'Non corretto', 'Bonferroni', 'HB', 'BH')
set(gca, 'Xtick', m, 'Xticklabel', medie_target_tutte_diverse(i_p))
xlabel('Valore della media contro cui testavamo')
ylabel('p value')
ylim([.1*0.05/n_confronti .06])
title(['p value contro metodi di correzione (n confronti = ' num2str(n_confronti) ')'])