Analizzare graficamente la distribuzione di una variabile

Questa volta l'obiettivo è di fornire uno strumento pratico, immediatamente riadattabile a nuovi dati, che condensa in una sola immagine quattro grafici che illustrano la distribuzione di una variabile: una sintesi grafica che rappresenta un utile complemento ai test di normalità [1].

Come dati impieghiamo la concentrazione delle ferritina nel sangue rilevata in 202 atleti australiani riportata nella variabile ferr della tabella ais inclusa nel pacchetto DAAG. Potete scaricare il pacchetto dal CRAN oppure acquisire la tabella seguendo le indicazioni alternative riportate in [2].

Incollate questo script nella Console di R e premete ↵ Invio:

# ANALIZZARE GRAFICAMENTE LA DISTRIBUZIONE DI UNA VARIABILE

#

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

mydata <- ais$ferr # salva in mydata i valori della variabile ferritina

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

#

# istogramma e kernel density plot

#

hist(mydata, xlim=c(0, 250), ylim=c(0, 0.012), freq=FALSE, breaks="FD", main="Istogramma\ne kernel density plot", xlab="Ferritina in µg/L", ylab="Densità di probabilità") # traccia l'istogramma dei dati

par(new=TRUE, ann=FALSE) # consente di sovrapporre il grafico successivo

plot(density(mydata, n=1024, from=0, to=250), xlim=c(0, 250), ylim=c(0, 0.012), yaxt="n", col="black") # sovrappone il kernel density plot della distribuzione campionaria

#

# distribuzione campionaria vs. gaussiana teorica

#

plot(density(mydata, n=1024, from=0, to=250), xlim=c(0, 250), ylim=c(0, 0.012)) # traccia il kernel density plot della distribuzione campionaria

par(new=TRUE, ann=FALSE) # consente di sovrapporre il grafico successivo

x <- seq(0, 250, length.out=1024) # calcola i valori in ascisse della gaussiana teorica

y <- dnorm(x, mean=mean(mydata), sd=sd(mydata)) # calcola i valori in ordinate della gaussiana teorica

plot(x, y, xlim=c(0, 250), ylim=c(0, 0.012), yaxt="n", col="red", type="l") # sovrappone la distribuzione gaussiana teorica in colore rosso

title(main="Distribuzione campionaria\nvs. gaussiana teorica", xlab="Ferritina in µg/L", ylab="Densità di probabilità") # aggiunge titolo e legende degli assi

#

# distribuzione cumulativa campionaria vs. distribuzione cumulativa teorica

#

par(new=FALSE, ann=TRUE) # per mostrare titolo e legende degli assi

plot(ecdf(scale(mydata)), main="Cumulativa campionaria\nvs. cumulativa teorica", xlab="Deviata normale standardizzata z", ylab="Frequenza cumulativa", xlog = FALSE, ylog = FALSE, xlim=c(-4, 4), ylim=c(0, 1), xaxp=c(-4, 4, 5), yaxp=c(0, 1, 5)) # traccia il grafico della distribuzione cumulativa campionaria

par(new=TRUE, ann=FALSE) # consente di sovrapporre il grafico successivo

plot(ecdf(rnorm(10000, mean=0, sd=1)), col="red", xlog=FALSE, ylog=FALSE, xlim=c(-4, 4), ylim=c(0, 1), xaxp=c(-4, 4, 5), yaxp=c(0, 1, 5)) # sovrappone la distribuzione cumulativa teorica in colore rosso

#

# quantili campionari vs. quantili teorici

#

par(new=FALSE, ann=TRUE) # per mostrare titolo e legende degli assi

qqnorm(scale(mydata), main="Quantili campionari\nvs. quantili teorici", xlab="Quantili teorici", ylab="Quantili campionari", xlim=c(-4, 4), ylim=c(-4, 4)) # traccia il grafico della distribuzione dei quantili campionari

abline (0, 1, col="red") # sovrappone la distribuzione dei quantili teorica in colore rosso

#

La ferritina è una proteina che, pur con alcune limitazioni, fornisce una misura dei depositi di ferro presenti nell'organismo, necessari per garantire una adeguata produzione di emoglobina. L'analisi grafica dei dati consente di confermare graficamente, per la concentrazione della ferritina nel sangue, una distribuzione non gaussiana.

Il codice dello script è illustrato dai commenti via via inseriti, ma vediamo brevemente il significato di questa rappresentazione.

Il primo grafico presenta la distribuzione dei dati sotto forma del classico istogramma [3] impiegando la funzione hist(), calcolando il numero di classi (breaks="FD") mediante la regola di Freedman-Diaconis [4], e ponendo in ordinate la densità di probabilità (freq=FALSE).

Questo consente di sovrapporre all'istogramma mediante la funzione plot(), con la stessa scala delle ordinate, il corrispondente kernel density plot, la cui densità di probabilità è calcolata mediante la funzione density() in corrispondenza di 1024 valori (n=1024) compresi tra 0 (from=0) e 250 (to=250).

Con il secondo grafico entriamo nel vivo del confronto tra la distribuzione dei dati osservata, o distribuzione empirica o distribuzione campionaria che dir si voglia, e quella che i dati dovrebbero avere se fossero distribuiti secondo una gaussiana. 

Questo viene fatto riportando prima il kernel density plot calcolato come al punto precedente, e riportando per confronto (in colore rosso) la distribuzione gaussiana teorica corrispondente, ottenuta calcolando con la funzione dnorm() la densità di probabilità (y) - quella che avremmo se dati con la media mean=mean(mydata) e la deviazione standard sd=sd(mydata) dei dati osservati fossero distribuiti in modo gaussiano - in corrispondenza di 1024 valori (length.out=1024) della x ottenuti con la funzione seq().

Il terzo grafico riporta la distribuzione cumulativa campionaria calcolata con la funzione ecdf() sui dati (mydata) del campione, e riporta per confronto (in colore rosso) la corrispondente distribuzione cumulativa teorica, calcolata sempre con la funzione ecdf() su un campione di 10 000 dati distribuiti in modo gaussiano con media 0 e deviazione standard 1, generati impiegando la funzione rnorm() con gli argomenti 10 000, mean=0, sd=1.

Il quarto grafico impiega la funzione qqplot() per calcolare i quantili della distribuzione campionaria in funzione dei quantili teorici, e ne confronta l'andamento con quello previsto (la linea retta in colore rosso) per dati che seguono una distribuzione gaussiana. Concettualmente, si tratta della trasformazione lineare del grafico della distribuzione cumulativa di cui al punto precedente. 

In definitiva vediamo come i grafici, realizzati in modo semplice, forniscono una indicazione concisa di come e di quanto una distribuzione campionaria si discosta da una distribuzione gaussiana. E, come già detto all'inizio, offrono una sintesi che rappresenta un indispensabile complemento ai test statistici di normalità (gaussianità) [1].

Aggiungo che, nonostante siano ovviamente tecnicamente ineccepibili, distribuzione cumulativa e quantili tendono a comprimere, dal punto di vista grafico, la differenza tra la distribuzione campionaria e la distribuzione attesa nel caso di dati distribuiti in modo gaussiano. Tale differenza per certi versi appare più evidente nel semplice confronto tra la distribuzione osservata, rappresentata sotto forma di kernel density plot e la corrispondente distribuzione gaussiana (grafico in alto a destra).

I dettagli delle funzioni e degli argomenti impiegati sono riportati nei post che trattano le rispettive rappresentazioni grafiche (vedere alla pagina Indice) e sono come sempre accessibili digitando help(nomedellafunzione) nella Console di R.

Si fa infine notare la comparsa nel codice R di un \n che impone un “a capo” in alcune delle righe di testo.

Per riutilizzare questo script per analizzare altri dati è sufficiente:

→ nella prima riga di codice sostituire ais con il nome che volete assegnare ai vostri dati e sostituire C:/Rdati/ais.csv con nome e posizione del file dal quale importare i dati, adeguando al bisogno il separatore dei campi (sep=";") e il separatore dei decimali (dec=",") a quelli da voi impiegati;

→ essendo (ad esempio) dati il nome dei vostri dati e var il nome della variabile che volete analizzare, sostituire nella seconda riga di codice ais$ferr con dati$var;

→ adattare opportunamente xlim, ylim, titoli e legende degli assi.

Vale infine la pena di ricordare che, impiegando l'utilità riportata nel post Salvare i grafici di R in un file, potete anche salvare le immagini realizzate con R, come quella qui riportata, sotto forma di file .bmp, .jpeg, .png, .pdf, .ps per poterle stampare, archiviare, inserire in una pubblicazione, in un post o in un sito web.

----------

[1] Vedere ad esempio la sintesi dei test statistici riportata nel post Tabulare una serie di test di normalità (gaussianità).

[2] Vedere nel post Il set di dati ais come acquisire i dati della tabella ais anche senza istallare il pacchetto DAAG. La concentrazione della ferritina è espressa in microgrammi per litro di siero (µL). 

[3] Da notare che a partire da R versione 4.0 le barre dell'istogramma di default sono riempite in colore grigio. Per una trattazione dettagliata, che include anche vari script con rappresentazoni utili nella pratica, vedere il post Istogrammi.

[4] Vedere: Wikipedia, the free encyclopedia. Freedman–Diaconis rule. URL consultato il 12/01/2023.