Analisi delle distribuzioni dei dati nei gruppi/cluster

Se trovate utile l'impiego delle tecniche di clustering (o analisi dei gruppi) [1] come metodo per l'analisi esplorativa dei dati – ma anche se le affrontate per la prima volta – potrebbe interessarvi questo esempio.

Copiate lo script, incollatelo nella Console di R e premete ↵ Invio:

# ANALIZZARE LE DISTRIBUZIONI DEI DATI NEI CLUSTER con ggpairs  

#  

library(cluster) # carica il pacchetto per il clustering

library(factoextra) # carica un pacchetto addizionale per il clustering

library(GGally) # carica il pacchetto che estende le funzioni di ggplot2

#

# preparazione dei dati

#

mydata <- read.table("c:/Rdati/ema.csv", header=TRUE, sep=";", dec=".") # importa i dati dei 643 casi

newdata <- subset(mydata, select=c("gr", "hb", "mcv", "hb_a2", "ferro")) # seleziona le variabili per il clustering

#

# clustering non gerarchico con il metodo di Rousseew (k-medoids)

#

myclust <- pam(newdata, 3, metric=c("euclidean"), stand=TRUE) # genera un oggetto 'pam' che contiene i risultati del clustering

fviz_cluster(myclust, myclust$cluster, labelsize=7, main = "Grafico dei cluster - metodo di Rousseew (k-medoids)") # grafico dei cluster per le prime due componenti principali

#

Vediamo di cosa si tratta. Abbiamo nel file ema.csv una raccolta di 643 casi dei quali abbiamo raccolto i risultati di tre analisi del sangue richieste a scopo diagnostico:

→ esame emocromocitometrico – del quale ci interessano in particolare i valori di: (i) gr la concentrazione degli eritrociti nel sangue (globuli rossi) espressa in migliaia di miliardi per litro di sangue 10¹²/L o in milioni per microlitro di sangue 10⁶/µL (le due espressioni sono numericamente equivalenti); (ii) hb la concentrazione dell'emoglobina espressa in grammi per decilitro di sangue g/dL; (iii) mcv il volume medio dei globuli rossi espresso in milionesimi di miliardesimo di litro (10-¹⁵) o femtolitri (fL);

→ emoglobina A2 – hb_A2 la variante fisiologica dell'emoglobina presente nell'adulto espressa come percentuale dell'emoglobina totale;

→ ferro – ferro la concentrazione del ferro espressa in microgrammi per decilitro di siero µg/dL. 

Sappiamo che i 643 casi comprendono:

→ soggetti sani, impiegati come controlli, con i risultati delle tre analisi all'interno dei valori fisiologici (valori "normali");

→ soggetti con alterazioni dell'esame emocromocitometrico e con carenza di ferro (confermata dalla diminuita concentrazione del ferro nel siero);

→  soggetti con alterazioni dell'esame emocromocitometrico e portatori del gene della beta-talassemia (confermato dall'aumento della percentuale di emoglobina A2);

→  soggetti con alterazioni dell'esame emocromocitometrico ma non classificabili ulteriormente con i dati disponibili.

Siamo interessati a una analisi esplorativa dei dati e, seguendo lo script, procediamo come segue.

1) Preparazione dei dati. 

Installiamo i tre pacchetti che ci servono, cluster, factoextra e GGally, quindi ci procuriamo i dati impiegati nell'esempio (per il file di dati trovate link e istruzioni alla pagina di Download):

→ effettuando il download del file di dati ema.csv

→ salvando il file nella cartella C:\Rdati\

Con la funzione read.table() importiamo i dati nell'oggetto mydata. Quindi con la funzione subset() riportiamo nell'oggetto newdata le sole variabili che ci interessano, cioè i risultati delle analisi indicate sopra (643 record), in cinque colonne:

> newdata

      gr   hb  mcv hb_a2 ferro

1   4.90 13.3 82.8   1.8   106

2   4.66 10.8 73.6   2.6   148

3   5.43 11.5 66.5   4.8   104

4   5.41 10.8 73.4   2.5    74

..........

642 5.28 11.0 65.4   1.9    11

643 4.03  7.2 55.0   1.9    21

2) Clustering non gerarchico (analisi dei gruppi). 

A questo punto eseguiamo il clustering non gerarchico, lo facciamo impiegando il metodo di Rousseew. La prima funzione utilizzata è pam() che deriva il suo nome da "Partitioning Around Medoids", impiega l'algoritmo di clustering dei dati  k-medoids (un versione robusta dell'algoritmo k-means) e prevede come argomenti:

→ i dati newdata da impiegare;

→ il numero di cluster, in questo caso 3, in cui classificare i dati;

→ la metrica da impiegare che è "euclidean" (in alternativa può essere "manhattan");

→ stand=TRUE che impone la standardizzazione dei dati, che viene eseguita automaticamente dalla funzione.

Con la successiva funzione fviz_cluster() è generato il grafico dei cluster (myclust$cluster) ottenuti con l'algoritmo k-medoids e contenuti nell'oggetto myclust.

Date le premesse riportate all'inizio abbiamo impostato 3 come numero dei gruppi (cluster) in cui classificare i dati allo scopo di separare (i) sani, (ii) carenza di ferro, (iii) beta-talassemia eterozigote. Consapevoli del fatto che data la casistica, la quantità e il tipo delle informazioni di cui disponiamo, le altre possibili e meno frequenti cause di alterazione dell'esame emocromocitometrico (e ve ne sono molte) rimarranno inestricabilmente sovrapposte tra loro e/o con qualcuno di questi tre gruppi (a causa dei vincoli imposti dalla informazione disponibile che è, come sempre, limitata).

3) Analisi delle distribuzioni dei dati nei cluster.

Eccoci arrivati al punto. Non ci basta avere, per ciascun cluster, l'indicazione dei soggetti (numerati in ordine progressivo in mydata) appartenenti a ciascun cluster. Vogliamo avere anche l'indicazione della corrispondenza tra i cluster e le distribuzioni dei risultati delle grandezze misurate. Lo facciamo in una nuova finestra grafica aperta con windows() associando a ciascun caso (record) della tabella newdata con la funzione cbind() il cluster di appartenenza myclust$clustering estratto dall'oggetto 'pam' myclust che contiene i risultati del clustering e che viene riportato aggiungendo una sesta colonna alla tabella.

#

# analisi delle distribuzioni dei dati nei cluster

#

windows() # apre e inizializza una nuova finestra grafica

newclust <- cbind(newdata, myclust$clustering) # riporta per ciascun caso il cluster di appartenenza

ggpairs(newclust, columns=1:5, aes(color=as.factor(myclust$clustering), alpha=0.5), title="Scatterplot e kernel density plot dei dati dei cluster") # scatterplot e kernel density plot dei dati nei cluster

# 

Ora abbiamo la tabella che possiamo impiegare per generare il grafico della distribuzione dei dati dei cluster con la funzione ggpairs() che impiega gli argomenti:

→ newclust per indicare la tabella che contiene i dati;

→ 1:5 per indicare che i dati da rappresentare sono contenuti nelle colonne da 1 a 5 (la sesta contiene il cluster di appartenenza);

→ as.factor() in quanto le variabili numeriche 1, 2 e 3 della colonna 6 che indicano il cluster di appartenenza devono essere trasformate in variabili fattore;

→ alpha=0.5 per rappresentare i colori con una trasparenza del 50%.

Da notare che al bisogno potete anche esportare la tabella newclust, che nelle prime cinque colonne contiene i dati e nella sesta il cluster di appartenenza come qui indicato 

> newclust

      gr   hb  mcv hb_a2 ferro myclust$clustering

1   4.90 13.3 82.8   1.8   106                  1

2   4.66 10.8 73.6   2.6   148                  1

3   5.43 11.5 66.5   4.8   104                  2

4   5.41 10.8 73.4   2.5    74                  3

..........

642 5.28 11.0 65.4   1.9    11                  3

643 4.03  7.2 55.0   1.9    21                  3

e salvarla sotto forma di file [2].

Ed ecco il grafico risultante.

La conclusione? 

1) l'analisi dei gruppi fornisce una sintesi nella quale le grandezze soggiacenti non sono più riconoscibili singolarmente in quanto le componenti principali sono una combinazione lineare di tutte le variabili originali (nel primo dei grafici la prima componente principale in ascisse è la variabile formata dalla combinazione lineare di tutte le variabili originali che spiega la maggior quantità di varianza, la seconda componente principale in ordinate è la variabile che spiega la maggior quantità di varianza in ciò che rimane una volta rimosso l'effetto della prima componente principale);

2) la distribuzione dei dati all'interno dei cluster fornisce una rappresentazione di dettaglio che consente di ricollegare questi ultimi con la distribuzione originaria delle variabili la cui combinazione ha dato loro origine, e consente di riconoscere come:

→ il cluster 1 (rosso) contiene i soggetti con emoglobina A2 e ferro normali;

→ il cluster 2 (verde) contiene i soggetti con emoglobina A2 aumentata; 

→ il cluster 3 (azzurro) contiene i soggetti con ferro diminuito.

Pur con distribuzioni che in parte si sovrappongono – un fatto determinato dalla incompletezza della informazione disponibile rispetto alla complessità dei casi – l'esempio illustra come:

→ il clustering fornisce una vista sintetica dell'informazione globale contenuta nei dati;

→ possiamo aggiungere all'analisi dei cluster una vista di dettaglio delle distribuzioni dei singoli dati all'interno dei cluster;

→ le due viste degli stessi dati si integrano sinergicamente l'una con l'altra. 

----------

[1] Fate click su clustering nelle Parole chiave.

[2] Vedere il post Esportazione dei dati in un file .csv.

Kernel density plot

I kernel density plot forniscono la rappresentazione grafica della densità di probabilità di una distribuzione o più precisamente della stima kernel di detta densità di probabilità e possono essere utilizzati in alternativa o accanto ai tradizionali istogrammi dei quali rappresentano una importante evoluzione. 

Il metodo non parametrico [1] impiegato per realizzare i kernel density plot – invece di suddividere i dati in classi e rappresentarli sotto forma di barre come negli istogrammi – prevede di collocare in corrispondenza di ciascun dato una piccola gobba (bump), determinata da un fattore di forma (kernel), e di sommare tutte le gobbe impiegando un fattore di smussamento detto ampiezza di banda (bandwidth). La curva risultante, il kernek densiy plot, riporta quindi sulle ascisse x i valori della variabile oggetto di studio espressi in una scala numerica continua [2] e sulle ordinate y la frequenza di tali valori espressa come [stima kernel della] densità di probabilità della distribuzione.

Da notare che gli esempi riportati sono stati realizzati senza la necessità di installare pacchetti aggiuntivi ma impiegando solamente le funzioni statistiche e grafiche di base di R.  

Iniziamo esaminando il primo fattore che condiziona l'aspetto di un kernel density plot, il fattore di forma (kernel) impiegando la concentrazione degli eritrociti (globuli rossi) nel sangue rilevata in 202 atleti australiani riportata nella colonna/variabile rcc della tabella ais inclusa nel pacchetto DAAG. Accertatevi di avere installato il pacchetto o in alternativa procedete come indicato in [3]. 

Copiate e incollate questo script nella Console di R e premete ↵ Invio:

# KERNEL DENSITY PLOT effetto del fattore di forma (kernel)

#

library(DAAG) # carica il pacchetto DAAG che include il set di dati ais

windows() # apre e inizializza una nuova finestra grafica

par(mfrow=c(2,2)) # suddivide la finestra in quattro quadranti,  uno per grafico

#

kdpg <-density(ais$rcc, kernel=c("gaussian"), bw=0.05, adjust=1) # calcola il kernel density plot 

plot(kdpg, frame.plot=FALSE, xlim=c(3,8), col="blue", lwd=1, main="kernel=gaussian", xlab="Eritrociti in 10^12/L", ylab="Stima kernel della densità di probabilità") # traccia il kernel density plot 

#

kdpe <-density(ais$rcc, kernel=c("epanechnikov"), bw=0.05, adjust=1) # calcola il kernel density plot 

plot(kdpe, frame.plot=FALSE, xlim=c(3,8), col="blue", lwd=1, main="kernel=epanechnikov", xlab="Eritrociti in 10^12/L", ylab="Stima kernel della densità di probabilità") # traccia il kernel density plot 

#

kdpr <-density(ais$rcc, kernel=c("rectangular"), bw=0.05, adjust=1) # calcola il kernel density plot 

plot(kdpr, frame.plot=FALSE, xlim=c(3,8), col="blue", lwd=1, main="kernel=rectangular", xlab="Eritrociti in 10^12/L", ylab="Stima kernel della densità di probabilità") # traccia il kernel density plot 

#

kdpt <-density(ais$rcc, kernel=c("triangular"), bw=0.05, adjust=1) # calcola il kernel density plot 

plot(kdpt, frame.plot=FALSE, xlim=c(3,8), col="blue", lwd=1, main="kernel=triangular", xlab="Eritrociti in 10^12/L", ylab="Stima kernel della densità di probabilità") # traccia il kernel density plot 

#

L'utilizzo del fattore di forma (kernel) è illustrato con quattro grafici, ciascuno generato impiegando:

1) la funzione density() con questi argomenti:

→ ais$rcc per indicare la variabile da analizzare;

→ kernel="..." per indicare il fattore di forma da applicare;

→ bw=... per indicare l'ampiezza di banda da applicare;

→ adjust=... per indicare il fattore di moltiplicazione da applicare a bw, cosa che ne facilita la personalizzazione;

2) la funzione plot() che ha come argomenti:

→ l'oggetto generato alla riga precedente dalla funzione density() e contenente i dati necessari per tracciare il grafico ;

→ frame.plot=FALSE per togliere la cornice rettangolare che delimita il grafico;

→ xlim=c(... , ...) per indicare valore minimo e massimo da applicare all'asse delle ascisse;

→ col="..." per applicare un colore specifico alla traccia del grafico (valore di default = "black");

→ lwd=... per indicare lo spessore della linea impiegata per rappresentare il grafico (valore di default = 1);

→ main="...", xlab="..." e ylab="..." per specificare titolo, etichetta dell'asse delle x ed etichetta dell'asse delle y.

Nella funzione density() il fattore di forma viene specificato con l'argomento kernel="..." che può assumere i seguenti valori: "gaussian", "epanechnikov", "rectangular", "triangular", "biweight", "cosine", "optcosine" [4].

Questo significa che l'algoritmo che genera il kernel density plot sostituisce ciascun valore della variabile in esame con una gobba che nei nostri quattro casi ha la forma di una curva gaussiana ("gaussian"), di un vertice di parabola con convessità diretta verso l'alto ("epanechnikov"), di un rettangolo ("rectangular"), di un triangolo ("triangular"), forme che risultano evidenti nei due valori isolati e singoli visibili sulla destra in ciascun grafico.

Nota bene: lasciando fissa l'ampiezza di banda bw=0.05 e variando il fattore di forma si vede che la somma di 202 rettangoli porta a un risultato più irregolare e frastagliato della somma di altrettante curve gaussiane, vertici di parabola o triangoli, ma in ogni caso il kernel density plot risultante mostra una distribuzione bimodale dei dati [5].  

Vediamo ora il secondo fattore che condiziona l'aspetto di un kernel density plot, e cioè l'ampiezza di banda (bandwidth). Copiate e incollate questo script nella Console di R e premete ↵ Invio:

# KERNEL DENSITY PLOT effetto dell'ampiezza di banda (bandwidth)

#

library(DAAG) # carica il pacchetto DAAG che include il set di dati ais

windows() # apre e inizializza una nuova finestra grafica

par(mfrow=c(2,2)) # suddivide la finestra in quattro quadranti,  uno per grafico

#

kdp50 <-density(ais$rcc, kernel=c("gaussian"), bw=0.05, adjust=1) # calcola il kernel density plot 

plot(kdp50, frame.plot=FALSE, xlim=c(3,8), col="blue", lwd=1, main="bandwidth=0.05", xlab="Eritrociti in 10^12/L", ylab="Stima kernel della densità di probabilità") # traccia il kernel density plot 

#

kdp75 <-density(ais$rcc, kernel=c("gaussian"), bw=0.075, adjust=1) # calcola il kernel density plot 

plot(kdp75, frame.plot=FALSE, xlim=c(3,8), col="blue", lwd=1, main="bandwidth=0.075", xlab="Eritrociti in 10^12/L", ylab="Stima kernel della densità di probabilità") # traccia il kernel density plot 

#

kdp100 <-density(ais$rcc, kernel=c("gaussian"), bw=0.1, adjust=1) # calcola il kernel density plot 

plot(kdp100, frame.plot=FALSE, xlim=c(3,8), col="blue", lwd=1, main="bandwidth=0.1", xlab="Eritrociti in 10^12/L", ylab="Stima kernel della densità di probabilità") # traccia il kernel density plot 

#

kdp125 <-density(ais$rcc, kernel=c("gaussian"), bw=0.125, adjust=1) # calcola il kernel density plot 

plot(kdp125, frame.plot=FALSE, xlim=c(3,8), col="blue", lwd=1, main="bandwidth=0.125", xlab="Eritrociti in 10^12/L", ylab="Stima kernel della densità di probabilità") # traccia il kernel density plot 

#

Nello script precedente abbiamo variato il fattore di forma lasciando fissa l'ampiezza di banda bw=0.05. Qui invece mantenendo fisso il fattore di forma kernel="gaussian" abbiamo fatto variare l'ampiezza di banda, portandola dal valore iniziale bw=0.05 prima a bw=0.075 poi a bw=0.1 e infine a bw=0.125. 

Risulta così documentato l'effetto dell'ampiezza di banda sul grafico: aumentandola aumenta lo smussamento del kernel density plot risultante, mentre diminuendola il grafico si adatta viepiù ai singoli punti e diventa perciò irregolare.

Nota bene: nel nostro caso il valore bw=0.05 porta ad un eccessivo adattamento della curva ai singoli punti (overfitting) rendendola irregolare mentre il valore bw=0.125 inizia a mostrare una minor separazione tra le due sottoclassi e quindi ci sta spostando verso un adattamento della curva ai dati troppo ridotto (underfitting) a causa dell'eccessivo fattore di smussamento. A questo punto la scelta tra un bw=0.075 che conserva un maggior dettaglio e un bw=0.1 che comunque mantiene l'informazione rilevante diventa una scelta personale, che dipende anche da quello che si intende comunicare e discutere in merito ai dati.

Dopo averne brevemente illustrato razionale e tecniche soggiacenti, possiamo passare a qualche esempio pratico di applicazione e personalizzazione dei kernel density plot realizzati impiegando solamente le funzioni statistiche e grafiche di base di R, aggiungendo al termine un esempio che impiega le funzioni contenute nel pacchetto aggiuntivo sm che consentono di realizzare facilmente kernel density plot multipli sovrapposti.

Gli esempi sono stati riportati in script separati, in modo che quelli che più interessano possano essere salvati indipendentemente l'uno dall'altro ed essere successivamente ripresi per analizzare i propri dati semplicemente eliminando la prima riga di codice e sostituendo ais$rcc con il nome della propria variabile.

Copiate e incollate questo primo script nella Console di R e premete ↵ Invio:

# KERNEL DENSITY PLOT SEMPLICE

#

library(DAAG) # carica il pacchetto DAAG che include il set di dati ais

windows() # apre e inizializza una nuova finestra grafica

#

plot(density(ais$rcc)) # traccia il kernel density plot

#

Il kernel density plot più semplice che è possibile realizzare con R prevede, dopo avere caricato i dati e aperto e inizializzato una nuova finestra grafica, una sola riga di codice nella quale:

→ la funzione density(), che ha come unico argomento il nome della variabile da analizzare ais$rcc, genera la stima kernel della densità di probabilità necessaria per tracciare il grafico;

→ la funzione plot() impiega i dati generati dalla funzione density() per tracciare il grafico.

Se digitate

str(density(ais$rcc))

potete visualizzare la struttura dei dati generati dalla funzione density() 

> str(density(ais$rcc))

List of 7

 $ x        : num [1:512] 3.37 3.38 3.39 3.39 3.4 ...

 $ y        : num [1:512] 0.000214 0.000251 0.000297 0.000349 0.000409 ...

 $ bw       : num 0.143

 $ n        : int 202

 $ call     : language density.default(x = ais$rcc)

 $ data.name: chr "ais$rcc"

 $ has.na   : logi FALSE

 - attr(*, "class")= chr "density"

mentre se digitate 

density(ais$rcc)$x

e

density(ais$rcc)$y

potete visualizzare rispettivamente in dettaglio tutti i 512 valori delle x e gli altrettanti valori delle y generati dalla funzione density() che la funzione plot() impiega per tracciare il grafico.

Nota bene: lo script contiene una soluzione minima che può essere utile quando si voglia effettuare una rapida analisi preliminare dei dati, ma non solo, perché la rappresentazione ottenuta con gli argomenti di default previsti da R fornisce, oltre al numero N dei dati, un'informazione importante: l'ampiezza di banda (bandwidth) impiegata per tracciare il grafico, che la funzione density() calcola automaticamente di volta in volta.

Per realizzare un kernel density plot personalizzato ora copiate e incollate questo  script nella Console di R e premete ↵ Invio:

# KERNEL DENSITY PLOT PERSONALIZZATO

#

library(DAAG) # carica il pacchetto DAAG che include il set di dati ais

windows() # apre e inizializza una nuova finestra grafica

#

kdp <-density(ais$rcc, kernel=c("gaussian"), bw=0.1, adjust=1) # calcola il kernel density plot 

plot(kdp, frame.plot=FALSE, xlim=c(3,8), ylim=c(0,1), main="Kernel density plot (kernel=gaussian, bw=0.1)", xlab="Eritrociti in 10^12/L", ylab="Stima kernel della densità di probabilità") # traccia il kernel density plot 

polygon(kdp, col="green", border="blue", lwd=1) # colora il kernel density plot

#

Qui la funzione density() è stata arricchita con gli argomenti:

→ ais$rcc per indicare la variabile da analizzare;

→ kernel="gaussian" per indicare che il fattore di forma da applicare è la curva gaussiana;

→ bw=0.1 per indicare l'ampiezza di banda da applicare;

→ adjust=1 per indicare il fattore di moltiplicazione da applicare a bw, cosa che in questo caso lascia bw immodificato.

La funzione plot() è stata personalizzata con gli argomenti già commentati nel primo script.

La novità è la funzione polygon() che consente di colorare il grafico specificando:

→ kdp per la variabile da analizzare;

→ col="green" per il colore di riempimento del kernel density plot;

→ border="blue" per il colore della linea che delimita il kernel density plot;

→ lwd=1 per lo spessore della linea che delimita il kernel density plot.

Digitate help(polygon) nella Console di R se volete approfondire la conoscenza di questa funzione.

A questo punto vediamo come fare una cosa che prima o poi viene in mente a tutti: sovrapporre un kernel density plot e l'istogramma corrispondente. Copiate e incollate questo script nella Console di R e premete ↵ Invio:

# KERNEL DENSITY PLOT CON L'ISTOGRAMMA CORRISPONDENTE e ordinate separate

#

library(DAAG) # carica il pacchetto DAAG che include il set di dati ais

windows() # apre e inizializza una nuova finestra grafica

par(mar=c(5.1,4.1,4.1,5.1)) # aumenta il margine destro da 2.1 a 5.1

#

hist(ais$rcc, xlim=c(3,8), ylim=c(0,40), freq=TRUE, col="chartreuse", breaks="FD", main="Kernel density plot con istogramma", xlab="Eritrociti in 10^12/L", ylab="Frequenza (numero dei casi)") # traccia l'istogramma 

par(new=TRUE, ann=FALSE) # consente la sovrapposizione del grafico successivo

#

kdp <-density(ais$rcc, kernel=c("gaussian"), bw=0.1, adjust=1) # calcola il kernel density plot 

plot(kdp, yaxt="n", xlim=c(3,8), ylim=c(0,1), col="blue", lwd=1) # traccia il kernel density plot

axis(4, col.ticks="blue", col.axis="blue") # riporta l'asse delle y del kernel density plot sulla destra

mtext("Stima kernel della densità di probabilità", side=4, line=3, cex=1, col="blue") # aggiunge la legenda all'asse delle y sulla destra

#

par(mar=c(5.1,4.1,4.1,2.1)) # ripristina i margini di default

#

Vogliamo sovrapporre istogramma e kernel density plot riportando due assi delle ordinate separati e con scale diverse:

→ l'asse delle y sulla sinistra riferito all'istogramma e che riporta la frequenza espressa come numero dei casi;

→ l'asse delle y sulla destra riferito al kernel density plot e che riporta la stima kernel della densità di probabilità.

Ecco le principali osservazioni sugli argomenti e sulle funzioni impiegate:

→ come prima cosa con la funzione par() il margine sulla destra del grafico (il quarto valore dell'argomento mar) viene aumentato dal valore di default 2.1 a 5.1 per creare lo spazio necessario per la legenda dell'asse delle y del kernel density plot che verrà posizionato appunto sulla destra;

→ nella funzione hist() viene riportato in ordinate il numero dei casi ponendo freq=TRUE;

→ il numero delle classi breaks viene calcolato con la regola di Freedman-Diaconis (FD);

→ la funzione par(new=TRUE, ann=FALSE) specifica che il grafico successivo verrà sovrapposto (new=TRUE) ma che non devono essere annotati automaticamente titolo e legende (ann=FALSE), perchè questo lo faremo con le righe successive;

→ con la funzione plot() il kernel density plot viene sovrapposto all'istogramma specificando con yaxt="n" che non deve essere riportato l'asse delle y che si sovrapporrebbe a quello dell'istogramma e che invece noi aggiungeremo ora sulla destra;

→ con la funzione axis() riportiamo l'asse delle y del kernel density plot nella posizione 4 cioè sulla destra (nella grafica di base di R i margini sono 1=margine inferiore, 2=margine sinistro, 3=margine superiore, 4=margine destro) e riportiamo in colore "blue" sia le tacche sia i valori della scala;

→ con la funzione mtext() aggiungiamo sulla destra (side=4) la legenda dell'asse delle y del kernel density plot, posizioniamo l'etichetta "Stima kernel della densità di probabilità" sulla line=3 (la linea 0 contiene le tacche, la 1 contiene i valori numerici, la 2 è volutamente lasciata vuota), riportiamo i caratteri nella dimensione standard (cex=1) e in colore blu (col="blue");

→ con la funzione par() ripristiniamo i margini di default, di fatto riportiamo quello di destra (il quarto valore dell'argomento mar) da 5.1 a 2.1.

Se preferite impiegare un'unica scala sia per l'istogramma sia per il kernel density plot, esprimendo per entrambi le ordinate in termini di densità di probabilità, copiate e incollate nella Console di R questo script e premete ↵ Invio:

# KERNEL DENSITY PLOT CON L'ISTOGRAMMA CORRISPONDENTE e ordinate comuni

#

library(DAAG) # carica il pacchetto DAAG che include il set di dati ais

windows() # apre e inizializza una nuova finestra grafica

#

hist(ais$rcc, xlim=c(3,8), ylim=c(0,1), freq=FALSE, breaks="FD", col="chartreuse", main="Kernel density plot con istogramma", xlab="Eritrociti in 10^12/L", ylab="Densità di probabilità") # traccia l'istogramma

par(new=TRUE, ann=FALSE) # consente la sovrapposizione del grafico successivo

#

kdp <-density(ais$rcc, kernel=c("gaussian"), bw=0.1, adjust=1) # calcola il kernel density plot 

plot(kdp, xlim=c(3,8), ylim=c(0,1), col="blue", lwd=1) # sovrappone all'istogramma il kernel density plot

#

Come si vede sono ora sufficienti (tolte le prime due) quattro righe di codice:

→ la prima per tracciare con la funzione hist() l'istogramma, per gli argomenti del quale si rimanda al bisogno al post Istogrammi;

→ la seconda per consentire con la funzione par() la sovrapposizione all'istogramma di un nuovo grafico;

→ la terza per calcolare con la funzione density() il kernel density plot;

→ la quarta per tracciare con la funzione plot() il kernel density plot calcolato con la funzione density().

Ma potete anche semplificare ulteriormente il codice impiegando questa sintetica variante:

# KERNEL DENSITY PLOT CON L'ISTOGRAMMA CORRISPONDENTE e ordinate comuni

#

library(DAAG) # carica il pacchetto DAAG che include il set di dati ais

windows() # apre e inizializza una nuova finestra grafica

#

hist(ais$rcc, xlim=c(3,8), ylim=c(0,1), freq=FALSE, breaks="FD", col="chartreuse", main="Kernel density plot con istogramma", xlab="Eritrociti in 10^12/L", ylab="Densità di probabilità") # traccia l'istogramma

#

lines(density(ais$rcc), xlim=c(3,8), ylim=c(0,1), col="blue", lwd=1) # sovrappone all'istogramma il kernel density plot

#

Come si vede sono ora sufficienti:

→ la funzione hist() per tracciare l'istogramma;

→ la funzione lines() per sovrapporre all'istogramma il kernel density plot calcolato con density(ais$rcc).

In questo caso il grafico è un po' più spartano e il kernel density plot risulta lievemente diverso dai due precedenti in quanto per semplificare al massimo il codice non viene imposto l'argomento bw=0.1 ma viene lasciato che sia R a calcolarne il valore. Se digitate

density(ais$rcc)$bw

vi viene mostrato il valore di banda passante calcolato automaticamente da R

> density(ais$rcc)$bw

[1] 0.1425658

e potete notare che in effetti è superiore a 0.1 cosa che comporta un maggior smussamento del kernel density plot, che vedete più arrotondato e con una minore separazione dei picchi rispetto ai due precedenti script. 

Nota bene: il kernel density plot proprio per il modo nel quale è calcolato presenta una distribuzione bimodale che ben si adatta a quella che si osserva nell'istogramma.

Facciamo un ulteriore passo avanti con il codice che sovrappone due kernel density plot in modo semplice e conciso. Copiate e incollate questo script nella Console di R e premete ↵ Invio:

# DUE KERNEL DENSITY PLOT SOVRAPPOSTI 

#

library(DAAG) # carica il pacchetto DAAG che include il set di dati ais

windows() # apre e inizializza una nuova finestra grafica

#

f <-as.matrix(subset(ais, sex=="f", select=c(rcc))) # sottoinsieme dei dati con sesso=f

m <-as.matrix(subset(ais, sex=="m", select=c(rcc))) # sottoinsieme dei dati con sesso=m

#

plot(density(f), xlim=c(3,7), ylim=c(0,1.5), col="red", lwd=1, main="Kernel density plot sovrapposti", xlab="Eritrociti in 10^12/L", ylab="Stima kernel della densità di probabilità") # traccia il kernel density plot del sesso=f

#

lines(density(m), xlim=c(3,7), ylim=c(0,1.5), col="green", lwd=1) # traccia il kernel density plot del sesso=m

#

legend(3, 1.5, legend=c("f","m"), col=c("red", "green"), pt.cex=2, pch=15) # aggiunge la legenda

#

Per la funzione subset() dobbiamo partire dalla tabella ais, se nella Console di R digitate 

str(ais)

vi comparirà la struttura della tabella

> str(ais)

'data.frame':   202 obs. of  13 variables:

 $ rcc   : num  3.96 4.41 4.14 4.11 4.45 4.1 4.31 4.42 4.3 4.51 ...

 $ wcc   : num  7.5 8.3 5 5.3 6.8 4.4 5.3 5.7 8.9 4.4 ...

 $ hc    : num  37.5 38.2 36.4 37.3 41.5 37.4 39.6 39.9 41.1 41.6 ...

 $ hg    : num  12.3 12.7 11.6 12.6 14 12.5 12.8 13.2 13.5 12.7 ...

 $ ferr  : num  60 68 21 69 29 42 73 44 41 44 ...

 $ bmi   : num  20.6 20.7 21.9 21.9 19 ...

 $ ssf   : num  109.1 102.8 104.6 126.4 80.3 ...

 $ pcBfat: num  19.8 21.3 19.9 23.7 17.6 ...

 $ lbm   : num  63.3 58.5 55.4 57.2 53.2 ...

 $ ht    : num  196 190 178 185 185 ...

 $ wt    : num  78.9 74.4 69.1 74.9 64.6 63.7 75.2 62.3 66.5 62.9 ...

 $ sex   : Factor w/ 2 levels "f","m": 1 1 1 1 1 1 1 1 1 1 ...

 $ sport : Factor w/ 10 levels "B_Ball","Field",..: 1 1 1 1 1 1 1 1 1 1 ...

che come si vede contiene 13 variabili. Le prime 11 sono variabili numeriche (num), le ultime due sono fattori (Factor), sono cioè variabili qualitative opportunamente codificate, la variabile sex per indicare il sesso impiegando 2 livelli di codifica ("f" e "m"), la variabile sport per indicare lo sport praticato, impiegando 10 livelli di codifica.

Con la funzione subset() sono estratti da ais i sottoinsiemi di dati f e m impiegando i seguenti argomenti:

→ ais per indicare la fonte dei dati da estrarre;

→ sex="f" ovvero sex="m" per estrarre da ais rispettivamente i dati corrispondenti agli atleti di sesso femminile e agli atleti di sesso machile;

→ select=(rcc) per estrarre da ais solamente i valori della variabile rcc che sono quelli che ci interessano per tracciare i kernel density plot.

I due kernel density plot sono tracciati con due righe di codice che impiegano:

→ la funzione plot() per tracciare il primo kernel density plot;

→ la funzione lines() per sovrapporre il secondo.

 

Da notare a proposito della funzione lines() che:

→ riporta sullo stesso grafico e unisce tra loro con segmenti di retta i punti con le coordinate x,y del kernel density plot calcolate dalla funzione density();

→ riporta come scala degli assi gli stessi valori xlim=c(3,7) e ylim=c(0,1.5) impiegati per il primo kernel density plot;

→ traccia i segmenti di retta che descrivono il kernel density plot con il colore col="green" e con una linea di spessore lwd=1.

La funzione legend() [6] identifica i dati rappresentati impiegando i seguenti argomenti:

→ 3 per l'ascissa cui allineare il lato sinistro della legenda; 

→ 1.5 per l'ordinata cui allineare l'angolo superiore della legenda;

→ legend=c(...) per l'elenco delle etichette che descrivono i dati:

→ col=c(...) per l'elenco dei colori corrispondenti ai dati

→ pt.cex=2 per indicare il fattore di ingrandimento dei simboli;

→ pch=15 per indicare il tipo di simbolo da impiegare (15=un quadrato).

Vediamo ora come sovrapporre due kernel density plot colorati e renderli visibili nella loro interezza, senza che uno nasconda parzialmente l'altro, utilizzando la trasparenza del colore. Copiate e incollate lo script nella Console di R e premete ↵ Invio:

# DUE KERNEL DENSITY PLOT SOVRAPPOSTI COLORATI

#

library(DAAG) # carica il pacchetto DAAG che include il set di dati ais

windows() # apre e inizializza una nuova finestra grafica

#

f <-as.matrix(subset(ais, sex=="f", select=c(rcc))) # sottoinsieme dei dati con sesso=f

m <-as.matrix(subset(ais, sex=="m", select=c(rcc))) # sottoinsieme dei dati con sesso=m

#

plot(density(f), xlim=c(3,7), ylim=c(0,1.5), main="Kernel density plot sovrapposti colorati", xlab="Eritrociti in 10^12/L", ylab="Stima kernel della densità di probabilità") # traccia il kernel density plot del sesso=f

polygon(density(f), col=rgb(255,0,0,127,maxColorValue=255), border="red", lwd=1) # colora il kernel density plot del sesso=f

#

lines(density(m), xlim=c(3,7), ylim=c(0,1.5)) # traccia il kernel density plot del sesso=m

polygon(density(m), col=rgb(0,255,0,127,maxColorValue=255), border="green", lwd=1) # colora il kernel density plot del sesso=m

#

legend(3, 1.5, legend=c("f","m"), col=c(rgb(255,0,0,127,maxColorValue=255), rgb(0,255,0,127,maxColorValue=255)), pt.cex=2, pch=15) # aggiunge la legenda

#

La funzione subset() impiegata per estrarre da ais i sottoinsiemi di dati f e m è stata illustrata nello script precedente.

Per sovrapporre i due kernel density plot colorati sono state impiegate le funzioni plot(), density(), polygon() e lines() anch'esse già illustrate in precedenza.

Nelle funzioni poligon() e legend() il colore è stato specificato con la funzione rgb(n,n,n,n) nella quale:

→ il primo numero indica il rosso (red);

→ il secondo il verde (green);

→ il terzo il blu (blu);

→ il quarto il grado di trasparenza del colore impiegato.

Ponendo maxColorValue=255 questi quattro valori possono essere espressi in una scala compresa tra 0 e 255. Così ad esempio i valori di rgb uguali a 255,0,0,127 del primo dei due kernel density plot stanno per rosso al 100%, no verde, no blu, trasparenza del colore uguale al 50%. Digitale help(rgb) per la documentazione completa della funzione.

Si possono anche realizzare due kernel density plot contrapposti. Copiate e incollate lo script nella Console di R e premete ↵ Invio:

# DUE KERNEL DENSITY PLOT CONTRAPPOSTI

#

library(DAAG) # carica il pacchetto DAAG che include il set di dati ais

windows() # apre e inizializza una nuova finestra grafica

#

f <-as.matrix(subset(ais, sex=="f", select=c(rcc))) # estrae il sottoinsieme di ais con sesso=f

m <-as.matrix(subset(ais, sex=="m", select=c(rcc))) # estrae il sottoinsieme di ais con sesso=m

par(mfrow=c(2,1)) # grafici organizzati in due righe e una colonna

#

par(mar=c(0,5,3,5)) # adatta i margini al grafico superiore

plot(density(f), frame.plot=FALSE, xaxt="n", yaxt="n", xlim=c(3,7), ylim=c(0,1.5), main="Kernel density plot contrapposti", xlab="", ylab="") # traccia il kernel density plot del sesso=f

axis(4) # riporta l'asse delle y del kernel density plot sulla destra

mtext("Densità di probabilità", side=4, line=3) # aggiunge la legenda all'asse delle y sulla destra

polygon(density(f), col=rgb(1,0,0,0.5), border="red", lwd=1) # colora il kernel density plot

legend("topleft", legend=c("f"), col=c(rgb(1,0,0,0.5)), pt.cex=2, pch=15) # riporta la legenda

#

par(mar=c(5,5,0,5)) # adatta i margini al grafico inferiore

plot(density(m), frame.plot=FALSE, xlim=c(3,7), ylim=c(1.5,0), main="", xlab="Eritrociti in 10^12/L", ylab="Densità di probabilità") # traccia il kernel density plot del sesso=m

polygon(density(m), col=rgb(0,1,0,0.5), border="green", lwd=1) # colora il kernel density plot

legend("right", legend=c("m"), col=c(rgb(0,1,0,0.5)), pt.cex=2, pch=15) # riporta la legenda

#

Si tratta di una estetica alternativa a quella dello script precedente, che può piacere o non piacere, e che viene realizzata in modo molto semplice con le funzioni:

→ plot() per tracciare il grafico;

→ polygon() per colorare il grafico;

→ legend() per aggiungere la legenda.

Per realizzare il grafico queste tre funzioni sono integrate con i seguenti accorgimenti:

→ par(mfrow=c(2,1)) per disporre i grafici in due righe per una colonna quindi in definitiva uno sopra l'altro;

→ par(mar=c(0,5,3,5)) per eliminare il margine inferiore e aumentare il margine destro al grafico superiore;

→ xaxt="n" per togliere l'asse delle x al grafico superiore (è sufficiente l'asse delle ascisse riportato al grafico inferiore);

→ yaxt="n" per togliere al grafico superiore l'asse delle y che per default verrebbe riportato sulla sinistra;

→ axis(4) per riportare l'asse delle y del grafico superiore sulla destra;

→ mtext() per aggiungere la legenda "Densità di probabilità" all'asse delle y di destra (side=4) sulla riga 3 (line=3);

→ par(mar=c(5,5,0,4)) per eliminare il margine superiore e aumentare il margine destro al grafico inferiore;

→ ylim=c(1.5,0) necessario per invertire nel grafico inferiore l'asse delle y rispetto al grafico superiore, nel quale abbiamo ylim=c(0,1,5).

Da notare che nelle funzioni polygon() e legend() per definire il colore viene impiegata la funzione rgb(n,n,n,n) ma, contrariamente a quanto avevamo fatto nello script precedente, questa volta non viene specificato l'argomento maxColorValue che pertanto qui assume il valore di default, che è uguale a 1. Quindi i valori n assegnati nell'ordine al rosso (red), al verde (green), al blu (blu) e al grado di trasparenza del colore, per definizione possono variare solamente tra 0 e 1. In questa scala il valore 0.5 del grado di trasparenza equivale al valore 127 che nello script precedente era espresso in una scala da 0 a 255.

Qualora si debbano sovrapporre molti kernel density plot, può risultare comodo impiegare le funzioni contenute nel pacchetto sm. Accertatevi di avere installato il pacchetto altrimenti dovete scaricarlo e installarlo dal CRAN. Ora copiate e incollate questo script nella Console di R e premete ↵ Invio:

# KERNEL DENSITY PLOT MULTIPLI SOVRAPPOSTI con il pacchetto sm

#

library(DAAG) # carica il pacchetto DAAG che include il set di dati ais

library(sm) # carica il pacchetto

windows() # apre e inizializza una nuova finestra grafica

#

sport <- c("B_Ball", "Field", "Gym", "Netball", "Row", "Swim", "T_400m", "T_Sprnt", "Tennis", "W_Polo") # vettore che riporta la codifica impiegata per la variabile/fattore sport

etichette <- c("Basketball", "Field", "Gymnastics", "Netball", "Rowing", "Swimming", "Track 400m", "Track Sprint", "Tennis", "Water Polo") # vettore con le etichette da assegnare agli sport

colori <- c("green", "red", "purple", "blue", "cyan", "grey50", "brown", "tomato", "magenta", "gold") # vettore con i colori da assegnare agli sport

#

sm.density.compare(ais$rcc, ais$sport, col=colori, xlab="Eritrociti in 10^12/L", ylab="Stima kernel di densità") # traccia i kernel density plot separati per sport

#

title(main="Concentrazione degli eritrociti per sport praticato") # aggiunge il titolo

legend(locator(1), levels(factor(ais$sport, levels=sport, labels=etichette)), fill=colori) # predispone la legenda con etichette e colori

#

# resta in attesa del click con il tasto sinistro del mouse nel punto in cui si desidera posizionare la legenda

#

Il punto di forza è che se i dati sono associati a una variabile qualitativa per la quale sono previsti N livelli di codifica, vengono generati automaticamente sia gli N kernel density plot corrispondenti sia la legenda necessaria per identificarli. Inoltre la legenda può essere posizionata nel grafico in un punto a piacere con un semplice click del mouse.

Nel nostro caso per rappresentare la concentrazione degli eritrociti rilevata nei 10 sport praticati dagli atleti australiani in altrettanti kernel density plot procediamo come segue:

→ nel vettore sport=c(...) riportiamo la codifica dei livelli impiegata nella variabile (fattore) sport;

→ nel vettore etichette=c(...) riportiamo nello stesso ordine le scritte che vogliamo impiegare per decrivere ciascuno sport;

→ nel vettore colori=c(...) riportiamo nello stesso ordine i colori che vogliamo applicare a ciacuno sport nel kernel density plot e nella legenda corrispondente;

→ con la funzione sm.density.compare() tracciamo i kernel density plot per la variabile ais$rcc suddividendoli per sport (ais$sport) impiegando i colori definiti nel vettore colori;

→ con la funzione title() riportiamo il titolo del grafico;

→ con la funzione legend() riportiamo riportiamo nella legenda, per la variabile (fattore) ais$sport, i livelli definiti nel vettore sport, con le etichette definite nel vettore etichette e i colori definiti nel vettore colori.

Una volta che sono comparsi i grafici, posizionate il puntatore del mouse dove volete che compaia la legenda e fate click con il tasto sinistro del mouse per farla comparire, terminando così lo script. Da notare che il codice riportato non consente di riposizionare dinamicamente la legenda: se non si è soddisfatti della sua posizione, è necessario rieseguire l’intero script e fare nuovamente click con il tasto sinistro del mouse nel nuovo punto in cui la si vuole posizionare.

Come al solito gli script riportati possono essere immediatamente riutilizzati sostituendo ais$rcc con il nome della propria variabile e con il semplice adattamento di titoli ed etichette degli assi [7].

__________

[1] I metodi parametrici sono basati sull'assunto che i dati abbiano una distribuzione gaussiana e utilizzano i "parametri" media (misura di posizione) e deviazione standard (misura di dispersione) nella sintesi dei dati e nell'esecuzione dei test statistici. I metodi non parametrici invece non prevedono assunti riguardo la distribuzione dei dati e per definizione non impiegano i citati "parametri" nella sintesi dei dati e nell'esecuzione dei test statistici. 

[2] Vedere il post Scale nominali, scale ordinali e scale numeriche.

[3] Vedere il post Il set di dati ais. La concentrazione degli eritrociti viene espressa in 10¹²/L o in 10⁶/µL, le due unità di misura impiegate nella pratica, che sono numericamente identiche. Nel post trovate anche come caricare i dati della tabella senza impiegare il pacchetto DAAG. 

[4] Qui vediamo i primi quattro fattori di forma, per i fattori "biweight", "cosine" e "optcosine", che forniscono risultati molto simili tra loro e praticamente identici a "gaussian", si rimanda alla documentazione di R.

[5] La distribuzione unimodale è per definizione quella della distribuzione gaussiana, la distribuzione bimodale tipicamente si osserva quando una variabile/serie di misure è composta da due sottoinsiemi aventi mode (ma anche medie) differenti, mentre una distribuzione plurimodale (o polimodale) è quella con mode multiple.

[6] Trovate altri esempi di legende nel post Aggiungere la legenda a un grafico.

[7] Un poco di lavoro in più è necessario per adattare lo script per la generazione di kernel density plot multipli sovrapposti con il pacchetto sm, nel quale ais$rcc deve essere sostituito con il nome della tabella$variabile da analizzare, ais$sport deve essere sostituito con il nome della tabella$variabile che contiene i fattori dai quali generare i kernel density plot, e infine devono essere adattati opportunamente il vettore sport che contiene i nomi dei fattori e i vettori etichette e colori.