%%
%% ESERCIZIO I
%  Caricare in memoria il dataset  kidney tramite l'istruzione
%  load kidney.mat
%  (attenzione si richiede che sia stata installato  il toolbox FSDA versione >=8.6.1)
%  Questo dataset contiene 203 righe e 12 colonne.
%
%  La prima colonna si riferisce alla classificazione data dai medici sulla  presenza di nefropatia (malattia renale).

% Le modalità della prima colonna sono
% 1= presenza di nefropatia
% 2= assenza di nefropatia

% Le variabili 2:12 contengono
%  informazioni sulle caratteristiche dei pazienti (età, pressione,
%  glicemia...).
%  Osservazione: i dati sono già stati standardizzati.
% 0) Calcolare il boxplot dell'età utilizzando come variabile di raggruppamento la classificazione fornita dai medici. Aggiungere all'asse y l'etichetta appropriata

load kidney.mat

% idxmedici = variabile classificatoria
idxmedici=kidney.cclass;

boxplot(kidney{:,2},idxmedici)
ylabel('Età (valori standardizzati)')

% Commentare il boxplot.

% Le persone che soffrono di nefropatia tendono ad avere un'età più elevata
% di quelli che non soffrono di nefropatia.

% Costruire la matrice dei diagrammi di disperisone con le variabili 2:4 utilizzando come variabile di raggruppamento  la classificazione fornita dai medici.

Xtable=kidney(:,2:4);
X=kidney{:,2:4};
spmplot(Xtable,idxmedici);



%  Commentare la matrice dei diagrammi di dispersione.

% Le persone che non soffrono di nefropatia tendono ad avere valori bassi
% di pressione e di glucosio
% Le variabili pressione e livello di glucosio sembrano scarsamente
% correlate. In generale, se non si considera la classificazione dei medici
% le tre variabili sembrano scarsamente correlate.

%  Applicare alle colonne 2:4 del dataset
%
%  1) la metodologia di cluster gerarchico del legame medio. Mostrare il
%  dendrogramma visualizzando tutte le aggregazione con l'opzione
%  'Orientation','right'. Colorare nel dendrogramma i rami utilizzando la
%  'ColorThreshold' di default. Quanti gruppi questa metodologia suggerisce
%  se si utilizza la 'ColorThreshold' di default? Utilizzando il numero di
%  gruppi suggerito dalla 'ColorThreshold' di default, assegnare ogni unità
%  ad un determinato gruppo. Utilizzando la funzione spmplot, commentare ed
%  interpretare i gruppi che vengono ottenuti. Illustrare se i gruppi
%  trovati corrispondono alla suddivisione fornita dai medici.
%  Calcolare l'indice di Rand aggiustato tra la classificazione fornita dai medici
%  e la classificazione ottenuta. Commentare l'indice di Rand.


% Metodo del legame medio
tree = linkage(X,'average');
% Motrare il dendrogramma visualizzando tutte le aggregazione
% con l'opzione 'Orientation','right'.
% Colorare nel dendrogramma i rami utilizzando la 'ColorThreshold' di default
dendrogram(tree,0,'Orientation','right','ColorThreshold','default');

% Quanti gruppi questa metodologia suggerisce
% se si utilizza la 'ColorThreshold' di default? 
% Se si utilizza la 'ColorThreshold' di default questa metodologia
% suggerisce la presenza di tre gruppi

% Utilizzando il numero di
% gruppi suggerito dalla 'ColorThreshold' di default, assegnare ogni unità
% ad un determinato gruppo. 

% La funzione cluster viene chiamata con il name/value MaxClust',3 per
% utilizzare 3 gruppi
idxLegameMedio = cluster(tree,'MaxClust',3);


% Utilizzando la funzione spmplot, commentare ed
% interpretare i gruppi che vengono ottenuti. 


% Utilizzando la funzione spmplot, commentare ed
%  interpretare i gruppi che vengono ottenuti.

% Commentare i gruppi che vengono ottenuti.
spmplot(Xtable,idxLegameMedio);
% I tre gruppi trovati a mio avviso non presentano alcun significato dal
% punto di vista medico

% Distribuzione di frequenza delle unità che appartengono ai diversi gruppi
tabulate(idxLegameMedio)
% Il gruppo 3 contiene 6 unità che presentano valori eccezionalmente
% elevati di glucosio
% Il gruppo 1 contiene 12 unità che presentano una valori di
% pressione/glucosio generalmente superiori alla media.
% Il gruppo 2 contiene tutte le rimanenti 185 unità. Chiaramente questa
% classificazione non è soffisfacente.

% Illustrare se i gruppi  trovati corrispondono alla suddivisione fornita dai medici.
%  Calcolare l'indice di Rand aggiustato tra la classificazione fornita dai medici
%  e la classificazione ottenuta. Commentare l'indice di Rand.

rindexLegameMedio=RandIndexFS(idxLegameMedio,idxmedici);
disp('Indice di Rand aggiustato tra la specie e la partizione derivante dall''applicazione del legame medio')
disp(rindexLegameMedio)
% Il valore molto basso dell'Indice di Rand conferma che la classificazione
% trovata è molto distante da quella dei medici

%  2) la metodologia delle k medie adattando come numero di gruppi k=2.
%  Commentare i risultati ottenuti. Illustrare se i due gruppi trovati
%  corrispondono alla suddivisione fornita dai medici.

% Metodologia delle k medie con k=2
[groupKM]=kmeans(X,2,"Replicates",100);
spmplot(Xtable,groupKM);
tabulate(groupKM)
% Il gruppo composto da 53 unità presenta valori generalmente elevanti di
% età e livello di glucosio
rindexgropKM=RandIndexFS(groupKM,idxmedici)
disp('Indice di Rand aggiustato tra la specie e la partizione derivante dall''applicazione di kmeans con k=2')
disp(rindexgropKM)
% Il grado di accordo tra la classificazione dei medici e quella prodotta dal 
% metodo delle k medie è intorno al 46 per cento

%  Utilizzando la funzione evalcluster determinare il numero ottimo di
%  gruppi tramite il criterio di Calinski e Harabaz. Come valore di k da
%  investigare scegliere 1,2, ..., 8.
%  Mostrare la matrice dei diagrammi di dispersione utilizzando come
%  variabile di raggruppamento l'assegnazione ai diversi gruppi in base al valore ottomale di k.
%  Calcolare l'indice di Rand aggiustato tra la classificazione fornita dai
%  medici e quella ottenuta. Commentare l'indice di Rand.

figure
eva = evalclusters(X,'kmeans','CalinskiHarabasz','KList',1:8);
plot(eva);
% Questo grafico mostra chiaramente che il numero ottimo dei gruppi è pari
% a 4


groupKMopt=eva.OptimalY;
% Distribuzione di frequenza dei gruppi che sono stati ottenuti
tabulate(groupKMopt)
spmplot(X,groupKMopt);
% I gruppi trovati presentano una forma sferica 
rindexgropKMopt=RandIndexFS(groupKMopt,idxmedici)
disp('Indice di Rand aggiustato tra la classificazione dei medici e il numero ottimo di gruppi utilizzando kmeans')
disp(rindexgropKMopt)

% Non c'è alcuna relazione tra la classificazione dei medici e la
% classificazione in base al numero ottimo di gruppi in base al metodo
% delle k-medie

%  3) la metodologia di tclust utilizzando come numero di gruppi k=2, un
%  livello di trimming pari a 0 (terzo argomento di input) ed imponendo
%  il vincolo che il rapporto tra il massimo e il minimo autovalore delle
%  matrici di covarianze non sia superiore a 2 (quarto argomento di
%  input). Calcolare l'indice di Rand aggiustato tra la
%  classificazione ottenuta e quella dei medici. Commentare l'indice di Rand.
%

livellotrimming=0;
rapportomaxmineig=2;
outTC=tclust(X,2,livellotrimming,rapportomaxmineig,'plots',1,'equalweights',0);
idxTC=outTC.idx;

rindexTCLUST=RandIndexFS(idxTC,idxmedici)
disp('Indice di Rand aggiustato tra la classificazione dei medici e l''ouptut di tclust')
disp(rindexTCLUST)
% Anche in questo caso non c'è relazione tra le due classificazioni


%% ESERCIZIO II
%
% La zona A1:B194 del file sigarette.xlsx
% contiene i  valori delle seguenti due variabili
% * posizione occupata all'interno dell'azienda
% * grado di dipendenza rispetto alle sigarette
% 0) Creare il grafico a barre per la variabile posizione occupata e calcolare la relativa distribuzione di frequenza.
%

%% Caricamento dati
clear 
close all
Xca=readtable('sigarette.xlsx','Range','A1:B194');

%% Distribuzione di frequenze variabile Posizione e grafico a barre
% Modo 1
% Trasformo la prima colonna di Xca in variabile categorica e applico la
% funzione summary e la funzione histogram
X1cat=categorical(Xca{:,1});
summary(X1cat);
histogram(X1cat)

% Modo 2
% Applico la funzione tabulate (a1 diventa un cellarray)
a1=tabulate(Xca{:,1});
% La seconda colonna di a1 viene trasformata in double
% freq è un vettore di doubles.
freq=cell2mat(a1(:,2));
bar(freq)
% Aggiungo le etichette sull'asse x (e le ruoto)
set(gca,'xtick',1:length(freq),'xticklabel',a1(:,1),'XTickLabelRotation',60)

% Modo 3
% Applico la funzione tabulate
a1=tabulate(Xca{:,1});
% La seconda colonna di a1 viene trasformata in double
freq=cell2mat(a1(:,2));
% L'istruzione categorical(a1(:,1),a1(:,1)) significa
% trasforma in categorico il vettore a(:,1) e prendi come ordine delle
% modalità da rappresentare nel grafico l'ordine specificato da a1(:,1))
xlabels=categorical(a1(:,1),a1(:,1));
bar(xlabels,freq)

%% PARTE NON PRESENTE NEL TESTO DELL'ESERCITAZIONE
% esploro la differenza tra l'istruzione 
% categorical(a1(:,1))
% e 
% categorical(a1(:,1),a1(:,1))
subplot(2,1,1)
xlabels=categorical(a1(:,1));
bar(xlabels,freq)
title('L''ordine delle barre è associato all''ordine alfabetico di xlabels')
subplot(2,1,2)

xlabels=categorical(a1(:,1),a1(:,1));
bar(xlabels,freq)
title('L''ordine delle barre è associato all''ordine specificato dentro a1(:,1)')

% 1) Calcolare e commentare l'indice tauyx  di riduzione proporzionale dell'eterogeneità  (con il rispettivo intervallo di confidenza) e l'indice di Cramer. Commentare i risultati.

%% Calcolo indici Cramer e tauyx
Indici=corrNominal(Xca(:,1:2),'datamatrix',true);
disp('Indice di Cramer')
disp(Indici.CramerV)
% Associazione pari al 17 per cento circa del valore massimo possibile
disp(['Indice tauyx è:' num2str(Indici.tauyx(1))])
disp(Indici.ConfLimtable('tauyx',:))
% Indice tauyx molto basso ma significativo.


% 2) Calcolare la tabella di contingenza tra le due variabili (variabile Posizione sulle righe) e l'inerzia totale.
% Denominare la variabile che contiene la tabella di contingenza tabCont;

%% Analisi delle corrispondenze
Acor=CorAna(Xca(:,1:2),'datamatrix',true');
disp('Tabella di contingenza tra posizione lavorative e grado di dipendenza dalle sigarette')
tabCont=Acor.Ntable;
disp(tabCont)

% 4) Creare e commentare il grafico di analisi delle corrispondenze.
% 5) Interpretare le due dimensioni latenti e fornire un giudizio sulla bontà dell'analisi

% Interpretazione delle due dimensioni latenti
% interpretazione dim 1 (da non fumatore a fumatore incallito)
% ==> propensione verso il fumo
% interpretazione dim 2 (manager in basso impiegati in alto)
% ==> Relazione inversa con il grado di importanza
% Inerzia spiegata vicino al 100 per cento. Rappresentazione ottima.

% 3) Calcolare la distanza (al quadrato) tra il profilo dei Senior Managers ed il profilo medio ponderata per la rispettiva massa (Inerzia). Calcolare manualemente questa distanza.


% Overview del profilo riga dei Senior Managers
Acor.OverviewRows('Senior_Managers',:)
disp('Distanza al quadrato ponderata per la massa per i SM')
Acor.OverviewRows('Senior_Managers','Inertia')

% Implementazione manuale di questa quantità
% Osservazione il profilo medio di riga è dato dalle masse di colonna
ProfiloRigaSM=Acor.ProfilesRows(1,:);
% ProfiloRigaMedio (come vettore riga) = vettore delle masse di colonna
% trasposto
ProfiloRigaMedio=Acor.c';

Distanza2=sum((ProfiloRigaSM-ProfiloRigaMedio).^2./ProfiloRigaMedio);
MassaRigaSM=Acor.OverviewRows{1,1};
% InerziaSM = inerzia dei senior manager = distanza euclidea al quadrato
% dal profilo medio, ponderata con le masse di colonna per la rispettiva massa di riga
% Si può interepretare come (x_i-M)^2 *f_i
InerziaSM=Distanza2*MassaRigaSM;
% La somma delle inerzie dei diverse punti riga è l'inerzia totale ossia
% la variabilità totale della tabella di contingenza

diffImplementazione=Acor.OverviewRows{'Senior_Managers','Inertia'}-InerziaSM;
assert(abs(diffImplementazione)<1e-12,"Le due implementazioni della distanza per la massa differiscono")


% 6) Utilizzando l'ordine per le righe Junior Employee,  Medium Employee, Senior Employee, Junior Manager, Senior Manager e per le colonne None, Light Medium, Heavy, calcolare l'indice gamma ed il relativo intervallo di confidenza.


tabCont1=tabCont([4 5 3 2 1],:);
corrOrdinal(tabCont1)
% L'intervallo di confidenza di gamma contiene il valore zero==>
% associazione tra le due variabili ordinali non significativa


%% ESERCIZIO III
% Utilizzando i dati presenti nella zona A1:F21 del file regioni.xlsx
% Standardizzare i dati.
% Sui dati standardizzati
% Calcolare i boxplot delle variabili in esame
% Costruire la matrice dei diagrammi di dispersione e la matrice di correlazione (con intestazionI di di righe e colonne)
% Costruire il biplot e commentarlo
% Utilizzando il metodo delle k-medie (con l'opzione numero di repliche pari a 100) e k=3 suddividere le regioni in 3 gruppi. Calcolare i centroidi.
% Proiettare i centroidi nello spazio delle prime due componenti principali. Commentare i gruppi che sono stati trovati

%% Caricamento dati
Xpca=readtable('regioni.xlsx','Range','A1:F21',...
    'ReadRowNames',true);
[Z,medie,sigma]=zscore(Xpca{:,:});
figure
boxplot(Z,'Labels',Xpca.Properties.VariableNames)
biplotFS(Xpca)
figure
outPCA=pcaFS(Xpca);

%% Matrice dei diagrammmi di dispersione
spmplot(Xpca);
% Calcolo della matrice di correlazione (in formato table)
cor=corr(Xpca{:,:});
disp(array2table(cor,'RowNames',Xpca.Properties.VariableNames,...
    'VariableNames',Xpca.Properties.VariableNames))

%% Analisi con 3 gruppi (metodo delle k medie)
[classif,Cent]=kmeans(Z,3,'Replicates',100);
% Osservazione. Cent contiene già i valori dei centroidi nello spazio delle
% variabili standardizzate

% Y = score non standardizzati (cov(Y) = matrice che contiene gli
% autovalori sulla diagonale principale)
Y=outPCA.score;
etirighe=Xpca.Properties.RowNames;

% Proietto i centroidi nello spazio delle prime due componenti principali
% In altri termini, applico ai centroidi le stesse trasformazioni che ho
% applicato alla matrice X per ottenere le componenti principali
% ossia standardizzo e moltiplico per i primi due autovettori
centroidiPCA=Cent*outPCA.coeff;

%% Rappresentazione degli score con etichette delle province e centroidi dei gruppi
figure
% Y sono le componenti principali non standardizzate
plot(Y(:,1),Y(:,2),'o')
text(Y(:,1),Y(:,2),etirighe)
% Aggiungo assi cartesiani
xline(0)
yline(0)

hold('on')
plot(centroidiPCA(:,1),centroidiPCA(:,2),'s','MarkerFaceColor','r','LineWidth',5)
xlabel('Prima componente principale (Dimensione della regione)');
ylabel('Seconda componente principale (Superficie (+) e densità (-))');


%% Rappresentazione degli score + centroidi + gruppo di appartenenza
close all
plot(Y(:,1),Y(:,2),'o')
text(Y(:,1),Y(:,2),num2str(classif))
% Aggiungo assi cartesiani
xline(0)
yline(0)

hold('on')
plot(centroidiPCA(:,1),centroidiPCA(:,2),'s','MarkerFaceColor','r','LineWidth',5)
xlabel('Prima componente principale (Dimensione della regione)');
ylabel('Seconda componente principale (Superficie (+) e densità (-))');


%% Rappresentazione degli score + centroidi + (gruppi di appartenenz etichetta provincia)
close all
idxAndRowLabels=strcat(num2str(classif),'--',etirighe);
plot(Y(:,1),Y(:,2),'o')
axislim=axis;
xlabel('Prima componente principale (Dimensione della regione)');
ylabel('Seconda componente principale (Superficie (+) e densità (-))');
text(Y(:,1),Y(:,2),idxAndRowLabels)
% Aggiungo assi cartesiani
xline(0)
yline(0)

hold('on')
plot(centroidiPCA(:,1),centroidiPCA(:,2),'s','MarkerFaceColor','r','LineWidth',5)


%
%