% In questo esercizio verifichiamo l'equivalenza tra ANOVA e regressione
% lineare nel caso di analisi a una via o a due vie. 

clc
clear all
close all

load data_anova.mat

% Descrizione dei dati: 
% bmi --> indice di massa corporea [kg/m^2]
% activity --> frequenza di attività fisica 
%   0 = mai
%   1 = 1-3 volte al mese
%   2 = una volta a settimana
%   3 = più di una volta a settimana
% sesso --> sesso (0 = femmina, 1 = maschio) 

%% Anova sulla variabile bmi usando la variabile activity per definire i gruppi

% Domanda: esiste una differenza significativa tra il bmi medio 
% di diversi gruppi di pazienti che eseguono attività fisica 
% con diversa frequenza?

% 1. Soluzione con la function di Matlab 
[p, anovatab, stats] = anova1(bmi,activity);

ind = find(strcmp(anovatab,"F"));
F = anovatab{ind+1};

disp('Anova a una via su bmi usando come fattore la variabile activity:')
disp('  1. Soluzione usando la funzione anova1: ')
disp(['    F = ' num2str(F)])
disp(['    p-value = ' num2str(p)])

% Rifiutiamo H0 --> almeno uno dei gruppi ha bmi medio diverso dagli altri

% 2. Soluzione facendo il conto a mano

% Numero di gruppi
N = length(unique(activity));

% Numero di osservazioni
n_tot = length(bmi);

% Media campionaria globale
mean_tot = mean(bmi);

% Calcolare numero di osservazioni per gruppo, media e varianza campionaria
% per gruppo

% Diversi livelli di attività fisica (etichette dei gruppi)
activity_levels = unique(activity); 
% Numero di osservazioni per gruppo
n_gruppo = zeros(N,1);
% Media campionaria del BMI per gruppo
mean_gruppo = zeros(N,1);
% Varianza campionaria del BMI per gruppo
varianza_gruppo = zeros(N,1);

for i = 1:N

    % Vettore dei valori di bmi dei soggetti appartenenti ad un gruppo
    bmi_i = bmi(find(activity==activity_levels(i)));

    n_gruppo(i) = length(bmi_i);
    mean_gruppo(i) = mean(bmi_i);
    varianza_gruppo(i) = var(bmi_i);

end

% Calcolo MSE
MSE = (sum((n_gruppo-1).*varianza_gruppo))/(n_tot-N);

% Calcolo MSB
MSB = (sum(n_gruppo.*(mean_gruppo-mean_tot).^2))/(N-1);

% Statistica F del test ANOVA
F = MSB/MSE;

% Livello di significatività
alpha = 0.05;

% p-value
p = 1-fcdf(F,N-1,n_tot-N);

disp('  2. Soluzione facendo il conto a mano:')
disp(['  F = ' num2str(F)])
disp(['  p-value = ' num2str(p)])

% 3. Soluzione con la regressione lineare
model = fitlm(activity,bmi,'CategoricalVars',[1]);
% Note: fitlm aggiunge automaticamente l'intercetta. Se specifichiamo che la
% variabile activity è qualitativa (categorica), fitlm crea in automatico le 
% variabili dummy per la sua codifica. 

% Test F che confronta il modello completo vs. il modello nullo
F = model.ModelFitVsNullModel.Fstat;  % Statistica F
p = model.ModelFitVsNullModel.Pvalue; % P-value

disp('  3. Soluzione usando la regressione lineare:')
disp(['  F = ' num2str(F)])
disp(['  p-value = ' num2str(p)])

%% Anova a una via usando la variabile sesso
disp(' ')
disp('Anova a una via su bmi usando come fattore la variabile sesso:')

% 1. Soluzione con la function anova1 di Matlab
[p, anovatab, stats] = anova1(bmi, sesso);

ind = find(strcmp(anovatab,"F"));
F = anovatab{ind+1};

disp('  1. Soluzione usando la funzione anova1: ')
disp(['    F = ' num2str(F)])
disp(['    p-value = ' num2str(p)])

% Non rifiutiamo H0 --> non viene rilevata una differenza significativa nel 
% valore medio di bmi di maschi e femmine 

% 2. Soluzione facendo un t test per il confronto di 2 campioni
% omoschedastici
bmi_0 = bmi(find(sesso==0)); % bmi delle femmine
bmi_1 = bmi(find(sesso==1)); % bmi dei maschi
[h p ci stats] = ttest2(bmi_1,bmi_0); 

disp('  2. Soluzione usando il t test per due campioni omoschedastici: ')
disp(['    t = ' num2str(stats.tstat)])
disp(['    p-value = ' num2str(p)])

% 3. Soluzione facendo una regressione lineare con la variabile sesso in
% ingresso
model = fitlm(sesso,bmi,'CategoricalVars',[1]);

% Test F che confronta il modello completo vs. il modello nullo
F = model.ModelFitVsNullModel.Fstat;  % Statistica F
p_F = model.ModelFitVsNullModel.Pvalue; % P-value

% T test che valuta se il coefficiente beta associato a sesso è
% significativamente diverso da 0
t = model.Coefficients.tStat(2);
p_t = model.Coefficients.pValue(2);

disp('  3. Soluzione usando la regressione lineare:')
disp(['    F = ' num2str(F)])
disp(['    p-value = ' num2str(p_F)])
disp(['    t = ' num2str(t)])
disp(['    p-value = ' num2str(p_t)])

% Note: il t test per due campioni omoschedastici si può fare anche con una
% regressione lineare avente come unica variabile di ingresso la variabile
% binaria che definisce i due gruppi. La statistica t del t test è la
% stessa del t test che valuta se il coefficiente beta associato alla
% variabile indipendente è diverso da 0.

%% Anova a due vie usando come fattori activity e sesso
disp(' ')
disp('Anova a due vie su bmi usando come fattori le variabili activity e sesso:')

% 1. Soluzione usando la function anovan di Matlab
[p,anovatab,~,~] = anovan(bmi,{activity, sesso});

ind = find(strcmp(anovatab,"F"));
F_activity = anovatab{ind+1};
F_sesso = anovatab{ind+2};

disp(' Usando la funzione anovan:')
disp('  Fattore activity:')
disp(['    F = ' num2str(F_activity)])
disp(['    p-value = ' num2str(p(1))])
disp('  Fattore sesso:')
disp(['    F = ' num2str(F_sesso)])
disp(['    p-value = ' num2str(p(2))])

% Il livello di attività fisica ha un impatto significativo sul bmi medio.
% Il sesso non ha un impatto significativo sul bmi medio. 

% 2. Soluzione usando la regressione lineare
model = fitlm([activity sesso], bmi,'CategoricalVars',[1 2]); 
% Note: specifico che le variabili activity e sesso sono entrambe
% qualitativa (categoriche)

% R2 del modello che include entrambi i fattori
R2 = model.Rsquared.Ordinary;
% Numero di variabili indipendenti del modello che include entrambi i fattori
m = 4;
% oppure: m = size(model.Coefficients,1)-1 

% Modello senza activity
model = fitlm([sesso], bmi,'CategoricalVars',[1]);
% R2 del modello senza activity
R2_red1 = model.Rsquared.Ordinary;
% Numbero di variabili indipendenti del modello senza activity
m_red1 = 1; 

% F per il fattore activity (statistica F dell'F test che confronta il
% modello con activity vs il modello senza activity)
F_activity = (R2-R2_red1)/(m-m_red1) / ((1-R2)/(n_tot-m-1));
p_activity = 1-fcdf(F_activity,m-m_red1,n_tot-m-1);

% Modello senza sesso
model = fitlm([activity], bmi,'CategoricalVars',[1]);
% R2 del modello senza sesso
R2_red2 = model.Rsquared.Ordinary;
% Numbero di variabili indipendenti del modello senza sesso
m_red2 = 3; 
% oppure: m = size(model.Coefficients,1)-1 

% F per il fattore sesso (statistica F dell'F test che confronta il
% modello con sesso vs il modello senza sesso)
F_sesso = (R2-R2_red2)/(m-m_red2) / ((1-R2)/(n_tot-m-1));
p_sesso = 1-fcdf(F_sesso,m-m_red2,n_tot-m-1);

disp(' 2. Soluzione usando la regressione lineare:')
disp('  Fattore activity:')
disp(['    F = ' num2str(F_activity)])
disp(['    p-value = ' num2str(p_activity)])
disp('  Fattore sesso:')
disp(['    F = ' num2str(F_sesso)])
disp(['    p-value = ' num2str(p_sesso)])