Istogrammi e kernel density plot

Gli istogrammi sono lo strumento tradizionale per la rappresentazione grafica dei dati di una distribuzione univariata. 

I kernel density plot, cioè i grafici basati sulla stima kernel di densità, possono essere utilizzati in alternativa al tradizionale istogramma, del quale rappresentano l'evoluzione. La stima kernel di densità è il metodo non parametrico impiegato per realizzare i kernel density plot. Invece di raccogliere le osservazioni in barre come negli istogrammi, lo stimatore kernel di densità colloca in corrispondenza di ogni osservazione una piccola ”gobba”, determinata da un fattore di forma (kernel), quindi somma tutte le gobbe e genera una curva smussata impiegando un parametro di smussamento detto ampiezza di banda (bandwidth). 

Qui riporto uno script che vuole giusto introdurre il tema, mentre trovate il necessario approfondimento, che include vari script con le rappresentazioni più frequentemente impiegate nella pratica, nei post:

→ Istogrammi 

→ Kernel density plot

→ Kernel density plot sovrapposti

→ Istogrammi affiancati, sovrapposti e contrapposti

Vediamo quindi una breve sintesi della rappresentazione di istogramma e kernel density plot di una variabile numerica continua, impiegando come esempio il set di dati ais.  Accertatevi di avere installato il pacchetto DAAG che contiene il set di dati, o in alternativa procedete come indicato [1].

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

# ISTOGRAMMI E KERNEL DENSITY PLOT

#

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)) # predispone la suddivisione della finestra in quattro quadranti, uno per grafico

#

hist(ais$rcc, xlim=c(3,7), main = "Istogramma semplice", xlab = "Eritrociti in 10^12/L", ylab = "Frequenza") # traccia un istogramma semplice dei dati

#

hist(ais$rcc, xlim=c(3,7), breaks=30, col="red", main = "Istogramma colorato", xlab = "Eritrociti in 10^12/L", ylab = "Frequenza") # traccia un istogramma colorato

#

hist(ais$rcc, xlim=c(3,7), freq=FALSE, breaks="FD", col="red", main = "Istogramma con curva gaussiana", xlab = "Eritrociti in 10^12/L", ylab = "Densità") # traccia un istogramma colorato

x <- seq(min(ais$rcc), max(ais$rcc), length.out=40) # curva gaussiana teorica valori in ascisse

y <- dnorm(x, mean=mean(ais$rcc), sd=sd(ais$rcc)) # curva gaussiana teorica valori in ordinate

lines(x, y, col="blue", lwd=2) # sovrappone all'istogramma la curva gaussiana

#

plot(density(ais$rcc), xlim=c(3,7), main = "Kernel density plot", xlab="Eritrociti in 10^12/L", ylab = "Stima kernel di densità") # traccia il kernel density plot

polygon(density(ais$rcc), col="gray82", border="black") # colora il kernel density plot

#

Lo script traccia tre differenti istogrammi e il kernel density plot per la variabile rcc, la concentrazione dei globuli rossi del sangue (eritrociti) [2].

Per favorire il confronto tra le diverse rappresentazioni, la finestra grafica aperta con la funzione windows() viene divisa con par(mfrow=c(2,2)) in quattro quadranti, uno per grafico, che sono riempiti nell'ordine da sinistra verso destra e dall'alto verso il basso. Ma con un semplice adattamento, e cioè eliminando la riga

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

e aggiungendo prima di ciascuno dei quattro grafici la riga

windows() # apre e inizializza una nuova finestra grafica

è possibile rappresentare i grafici separatamente.

Il primo dei tre istogrammi è il più semplice che si possa realizzare, specifica come argomenti solamente i dati dei quali tracciare l'istogramma (ais$rcc), limite inferiore e superiore dell'asse delle x (xlim=c(3,7)) per rendere la rappresentazione omogenea con gli altri istogrammi, più il titolo (main = "...") e le etichette dell'asse delle x (xlab = "...") e dell'asse delle y (ylab = "..."). Per tutti gli altri argomenti sono impiegati i valori di default della funzione hist() [3].

Il secondo istogramma prevede in aggiunta rispetto al primo solamente l'argomento che specifica il numero di classi (breaks=30), ovviamente adattabile secondo le esigenze [4], e l'argomento che specifica il colore dell'istogramma (col="red").

Il terzo istogramma viene realizzato impiegando l'argomento breaks="FD" con il quale il numero delle classi viene calcolato impiegando la regola di Freedman-Diaconis [5], che consente di minimizzare la differenza tra l'area dell'istogramma e l'area della gaussiana corrispondente. Ponendo in ordinate, invece della frequenza, la densità (di probabilità) con l'argomento freq=FALSE è possibile sovrapporre all'istogramma la curva gaussiana teorica corrispondente, avente cioè la media mean=mean(ais$rcc) e la deviazione standard sd=sd(ais$rcc) della variabile indagata. Da notare come in questo caso la curva gaussiana mal si adatta alla effettiva distribuzione dei dati, che sembra essere bimodale.

Il quarto grafico è il kernel density plot. Per tracciarlo viene impiegata la funzione plot() il cui argomento principale non è più costituito come nel caso degli istogrammi dalla variabile ais$rcc bensì dalla sua densità kernel calcolata con la funzione density(). La successiva funzione polygon() colora il kernel density plot in una delle numerose possibili tonalità di colore grigio ("gray82") e ne traccia il profilo in colore nero ("black"). Da notare che il kernel density plot, contrariamente all'istogramma semplice in alto a sinistra, rende immediatamente evidente la distribuzione bimodale dei dati.

Vediamo ora come realizzare istogramma e kernel density plot sovrapposti con il pacchetto ggplot2, che deve essere scaricato e installato dal CRAN).

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

# ISTOGRAMMA E KERNEL DENSITY PLOT SOVRAPPOSTI

#

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

library(ggplot2) # carica il pacchetto per la grafica

#

windows() # apre e inizializza una nuova finestra grafica

#

ggplot(data.frame(ais$rcc), aes(x = ais$rcc)) + geom_histogram(aes(y = after_stat(density)), fill = 'red', alpha = 0.5, bins = 19) + geom_density(colour = 'blue') + xlab(expression(bold('Eritrociti in 10^12/L'))) + ylab(expression(bold('Densità'))) + xlim(3,7) + ylim(0,1) # traccia istogramma e kernel density plot sovrapposti

#

Come vedete nello script il grafico viene realizzato mediante la funzione ggplot() [6] concatenando (+) più funzioni in una sola riga di codice:

Ma potete anche semplificare il codice impiegando questa variante - basata su due sole funzioni, hist() e lines() - che non richiede l'installazione di un pacchetto aggiuntivo per la grafica:

# KERNEL DENSITY PLOT CON L'ISTOGRAMMA CORRISPONDENTE

#

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,7), ylim=c(0,1), freq=FALSE, breaks="FD", col="red", 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,7), ylim=c(0,1), col="blue", lwd=1) # sovrappone all'istogramma il kernel density plot

#

E questo è il risultato:

L'istogramma e il kernel density plot sono un ottimo esempio di come senza una adeguata rappresentazione grafica potrebbe passare inosservato che nel set di dati ais la concentrazione degli eritrociti nel sangue segue una distribuzione bimodale, dovuta alla coesistenza di due sottoinsiemi di dati aventi differente moda (e media) e correlati al sesso.

Ricordo anche la possibilità di realizzare kernel density plot sovrapposti [6], separati per ciascuno dei fattori impiegati per la classificazione dei casi (che nel set di dati ais sono sesso e sport praticato).

Come al solito il riutilizzo del codice di entrambi gli script è molto semplice in quanto richiede solamente di sostituire alla variabile ais$rcc la variabile che volete analizzare, adattando poi opportunamente titolo e legende del grafico.

----------

[1] Vedere il post Il set di dati ais nel quale trovate anche come caricare i dati della tabella senza impiegare il pacchetto DAAG.

[2] Concentrazione espressa in migliaia di miliardi di globuli rossi per litro di sangue (10¹²/L), il che è numericamente identico ad esprimerla in milioni di globuli rossi per microlitro (milionesimo di litro) di sangue (10⁶/µL), come può essere riportata in alternativa.

[3] Per gli altri argomenti che è possibile impiegare con la funzione hist() digitare help(hist) nella Console di R. Da notare che a partire da R versione 4.0 le barre dell'istogramma sono di default riempite in colore grigio.

[4] In ogni caso il numero delle classi viene poi arrotondato automaticamente da R, per la regola utilizzata vedere la funzione pretty().

[5] Vedere Wikipedia. Freedman–Diaconis rule. URL consultato il 05/01/2023.

[6] Per i dettagli relativi alle funzioni e agli argomenti impiegati vedere sul CRAN il manuale del pacchetto: Package 'ggplot2'. URL aggiornato il 14/11/2022: Package ‘ggplot2’