% 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?

% 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('  Usando la funzione anova1: ')
disp(['    F = ' num2str(F)])
disp(['    p-value = ' num2str(p)])

% Conto a mano

n_tot = length(bmi); % numero totale di campioni

% Calcolo varianze campionarie e medie campionarie degli 4 gruppi
activity_levels = [0 1 2 3]; 
N = length(activity_levels); % Numero di gruppi
varianze_gruppi = zeros(1,N); % Varianze degli N gruppi
n_gruppi = zeros(1,N);  % Numero di campioni negli N gruppi
medie_gruppi = zeros(1,N);  % Medie degli N gruppo

for i = 1:N
    gruppo = activity_levels(i);
    ind = find(activity==gruppo);
    bmi_gruppo = bmi(ind);
    varianze_gruppi(i) = var(bmi_gruppo);
    n_gruppi(i) = length(ind);
    medie_gruppi(i) = mean(bmi(ind));
end

% Calcolo di MSE
MSE = sum((n_gruppi-1).*varianze_gruppi)/(n_tot-N);

% Calcolo di MSB
media_tot = mean(bmi);
MSB = sum(n_gruppi.*(medie_gruppi-media_tot).^2)/(N-1);

% Calcolo della statistica F
F= MSB/MSE;

alpha = 0.05;  % Livello di significatività
Fc = finv(1-alpha,N-1,n_tot-N);  % Soglia critica

disp(' Conto a mano:')
disp(['  F = ' num2str(F)])
disp(['  Fc = ' num2str(Fc)])
if F>Fc
    disp('  Rifiuto H0')
else
    disp('  Non posso rifiutare H0')
end

% Conto con la regressione lineare
model = fitlm(activity,bmi,'CategoricalVars',[1]);

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

disp(' Usando la regressione lineare:')
disp(['  F = ' num2str(F)])
disp(['  p-value = ' num2str(p)])

%% Anova a due vie usando come fattori activity e sesso

% Usando la function di Matlab
factors = {activity, sesso};
[p,anovatab,~,~] = anovan(bmi, factors);

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

disp(' ')
disp('Anova a due vie su bmi usando activity e sesso come fattori:')
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))])

% Usando la regressione lineare
model = fitlm([activity sesso],bmi,'CategoricalVars',[1 2]);
R2 = model.Rsquared.ordinary; 
m = 4; 

% Modello senza activity
model = fitlm(sesso,bmi,'CategoricalVars',[1]);
R2_red = model.Rsquared.ordinary; 
m_red = 1;

% Statistica F dell' F test per confrontare il modello con/senza activity
F_activity = ((R2 - R2_red)/(m-m_red)) / ((1-R2)/(n_tot-m-1));

% Modello senza sesso
model = fitlm(activity,bmi,'CategoricalVars',[1]);
R2_red = model.Rsquared.ordinary; 
m_red = length(activity_levels)-1;

% Statistica F dell' F test per confrontare il modello con/senza sesso
F_sesso = ((R2 - R2_red)/(m-m_red)) / ((1-R2)/(n_tot-m-1));

disp(' Usando la regressione lineare:')
disp('  Fattore activity:')
disp(['    F = ' num2str(F_activity)])
disp('  Fattore sesso:')
disp(['    F = ' num2str(F_sesso)])

%% Anova a una via usando la variabile sesso

% Anova usando variabile sesso
[p_anova anovatab stats] = anova1(bmi,sesso);
ind = find(strcmp(anovatab,'F'));
F_anova = anovatab{ind+1};

% T test usando variabile sesso
bmi_maschi = bmi(find(sesso==1));
bmi_femmine = bmi(find(sesso==0));
[h, p_ttest] = ttest2(bmi_femmine,bmi_maschi,'Vartype','equal');

% Regressione lineare usando sesso
model = fitlm(sesso,bmi);
F_reg = model.ModelFitVsNullModel.Fstat;  % Statistica F
p_reg = model.ModelFitVsNullModel.Pvalue; % P-value

disp(' ')
disp('Anova a una via su bmi usando come fattore la variabile sesso:')
disp('  Usando la function anovan:')
disp(['    F = ' num2str(F_anova)])
disp(['    p-value = ' num2str(p_anova)])
disp('  Usando il t test:')
disp(['    p-value = ' num2str(p_ttest)])
disp('  Usando la regressione lineare:')
disp(['    F = ' num2str(F_reg)])
disp(['    p-value = ' num2str(p_reg)])