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: Freedman–Diaconis rule.

https://en.wikipedia.org/wiki/Freedman%E2%80%93Diaconis_rule

Statistiche elementari parametriche

Siamo al terzo passo delle valutazioni che è sempre necessario effettuare nell'ambito di un percorso logico che prevede, per il calcolo delle statistiche elementari di una singola variabile (analisi univariata):

→ analisi esplorativa dei dati;

→ esecuzione dei test di normalità (gaussianità) per valutare se i dati seguono una distribuzione gaussiana;

→ calcolo delle statistiche elementari parametriche (media, deviazione standard, varianza, quantili parametrici) se i dati seguono una distribuzione gaussiana;

→ calcolo delle statistiche elementari non parametriche (mediana, deviazione assoluta mediana o MAD, quartili e altri quantili non parametrici) se i dati non seguono una distribuzione gaussiana.

Si tratta del passo conclusivo per la variabile altezza (espressa in cm) rilevata in 202 atleti australiani, in quanto con i test di normalità abbiamo stabilito che la variabile è distribuita in modo gaussiano: questa era la sintesi grafica dei risultati.

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

# STATISTICHE ELEMENTARI PARAMETRICHE 

#

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

mydata <- as.matrix(ais[c(10)]) # ci interessa la sola colonna con l'altezza

#

# statistiche elementari parametriche 

#

min(mydata) # valore minimo

mean(mydata) # media 

max(mydata) # valore massimo

range(mydata) # range 

max(mydata)-min(mydata) # campo di variazione

var(mydata) # varianza

sd(mydata) # deviazione standard

(sd(mydata)/mean(mydata))*100 # coefficiente di variazione CV %

#

# confronta media e mediana

#

media <- mean(mydata) # calcola la media

mediana <- median(mydata) # calcola la mediana

data.frame(media, mediana) # le mette a confronto

#

# confronta deviazione standard e MAD

#

ds <- sd(mydata) # calcola la deviazione standard

mad <- mad(mydata) # calcola la Median Absolute Deviation (about the median) o MAD

data.frame(ds, mad) # le mette a confronto

#

# confronta quantili parametrici e non parametrici

#

qpar <- round(qnorm(c(seq(0.025, 0.975, 0.025)), mean=mean(mydata), sd=sd(mydata)), digits=2) # calcola i quantili parametrici

qnon <- round(quantile(mydata, probs=seq(0.025, 0.975, 0.025)), digits=2) # calcola i quantili non parametrici

dper <- abs(round(100*(qpar-qnon)/((qpar+qnon)/2), digits=2)) # calcola la differenza in percentuale

data.frame(qpar, qnon, dper) # li mette a confronto

#

Non avrebbe senso ripetere l'analisi esplorativa dei dati, i test di normalità e la sintesi grafica dei risultati - per questi si rimanda agli script riportati nel percorso logico indicato all'inizio. Per cui qui ci limitiamo a calcolare le usuali statistiche elementari parametriche, riportando per confronto i risultati delle statistiche elementari non parametriche corrispondenti. 

Dopo avere caricato il pacchetto DAAG [1] con la tabella ais, la variabile che ci interessa, contenuta nella colonna 10 della tabella ais viene assegnata (<-) all'oggetto mydata, sul quale sono calcolate le statistiche elementari parametriche con le funzioni:

→ min() per riportare il valore minimo osservato;

→ mean() per calcolare la media;

→ max() per riportare il valore massimo osservato;

→ range() per calcolare il range, che qui è correttamente;

→ max(mydata)-min(mydata) per calcolare il campo di variazione, talora indicato (impropriamente) come "range";

→ var() per calcolare la varianza,

→ sd() per calcolare la deviazione standard;

→ e con l'espressione (sd(mydata)/mean(mydata))*100 che calcola il coefficiente di variazione percentuale.

> min(mydata) # valore minimo

[1] 148.9

> mean(mydata) # media 

[1] 180.104

> max(mydata) # valore massimo

[1] 209.4

> range(mydata) # range 

[1] 148.9 209.4

> max(mydata)-min(mydata) # campo di variazione

[1] 60.5

> var(mydata) # varianza

         ht

ht 94.76038

> sd(mydata) # deviazione standard

[1] 9.734494

> (sd(mydata)/mean(mydata))*100 # coefficiente di variazione CV %

[1] 5.404931

Sono poi messe a confronto la media con la mediana

> media <- mean(mydata) # calcola la media

> mediana <- median(mydata) # calcola la mediana

> data.frame(media, mediana) # le mette a confronto

    media mediana

1 180.104   179.7

la deviazione standard con la MAD [2]

> ds <- sd(mydata) # calcola la deviazione standard

> mad <- mad(mydata) # calcola la Median Absolute Deviation (about the median) o MAD

> data.frame(ds, mad) # le mette a confronto

        ds     mad

1 9.734494 8.96973

e i quantili parametrici con i quantili non parametrici, riportando nell'ultima colonna la differenza percentuale tra i due

> qpar <- round(qnorm(c(seq(0.025, 0.975, 0.025)), mean=mean(mydata), sd=sd(mydata)), digits=2) # calcola i quantili parametrici

> qnon <- round(quantile(mydata, probs=seq(0.025, 0.975, 0.025)), digits=2) # calcola i quantili non parametrici

> dper <- abs(round(100*(qpar-qnon)/((qpar+qnon)/2), digits=2)) # calcola la differenza in percentuale

> data.frame(qpar, qnon, dper) # li mette a confronto

        qpar   qnon dper

2.5%  161.02 158.98 1.28

5%    164.09 163.96 0.08

7.5%  166.09 167.35 0.76

10%   167.63 169.17 0.91

12.5% 168.91 170.36 0.85

15%   170.01 171.40 0.81

17.5% 171.01 172.22 0.71

20%   171.91 172.76 0.49

22.5% 172.75 173.52 0.44

25%   173.54 174.00 0.26

27.5% 174.29 174.18 0.06

30%   175.00 175.00 0.00

32.5% 175.69 175.73 0.02

35%   176.35 176.07 0.16

37.5% 177.00 177.30 0.17

40%   177.64 177.94 0.17

42.5% 178.26 178.44 0.10

45%   178.88 178.99 0.06

47.5% 179.49 179.60 0.06

50%   180.10 179.70 0.22

52.5% 180.71 180.10 0.34

55%   181.33 180.36 0.54

57.5% 181.94 181.00 0.52

60%   182.57 182.66 0.05

62.5% 183.21 183.06 0.08

65%   183.85 183.90 0.03

67.5% 184.52 184.67 0.08

70%   185.21 185.17 0.02

72.5% 185.92 185.60 0.17

75%   186.67 186.18 0.26

77.5% 187.46 187.28 0.10

80%   188.30 188.62 0.17

82.5% 189.20 189.18 0.01

85%   190.19 190.67 0.25

87.5% 191.30 191.94 0.33

90%   192.58 192.79 0.11

92.5% 194.12 193.86 0.13

95%   196.12 195.17 0.49

97.5% 199.18 197.48 0.86

Le differenze trascurabili tra statistiche elementari parametriche e statistiche elementari non parametriche riflettono il fatto che nel caso di una distribuzione perfettamente gaussiana [3] le due coincidono - è quindi assolutamente lecito impiegare le seconde al posto delle prime. Mentre il contrario è sbagliato - nel caso di una distribuzione non gaussiana le statistiche elementari parametriche non devono essere impiegate.

Lo script può essere riutilizzato molto semplicemente, assegnando (<-) all'oggetto mydata i propri dati.

----------

[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] La Median Absolute Deviation about median o MAD ovvero la deviazione assoluta mediana (dalla mediana) è l'equivalente non parametrico della deviazione standard. Vedere: Rousseeuw PJ, Croux C. Alternatives to the Median Absolute Deviation. Journal of the American Statistical Association 88 (424), 1273-1283, 1993.

https://www.jstor.org/stable/2291267?seq=1#page_scan_tab_contents

[3] Un esempio di dati distribuiti in modo perfettamente gaussiano è riportato nel post La distribuzione gaussiana.

Test di normalità (gaussianità)

Siamo al secondo passo delle valutazioni preliminari che è sempre necessario effettuare nell'ambito di un percorso logico che prevede, per il calcolo delle statistiche elementari di una singola variabile (analisi univariata), i seguenti passi:

→ analisi esplorativa dei dati;

→ esecuzione dei test di normalità (gaussianità) per valutare se i dati seguono una distribuzione gaussiana [1];

→ calcolo delle statistiche elementari parametriche (media, deviazione standard, varianza, quantili parametrici) se i dati seguono una distribuzione gaussiana [2];

→ calcolo delle statistiche elementari non parametriche (mediana, deviazione assoluta mediana o MAD, quartili e altri quantili non parametrici) se i dati non seguono una distribuzione gaussiana.

I dati sono contenuti nella tabella ais del pacchetto DAAG - accertatevi di avere installato il pacchetto o in alternativa procedete come indicato in [3]. Dovete scaricare e installare anche i pacchetti moments e nortest. Trovate la loro documentazione completa, incluso il manuale di riferimento, sul repository della documentazione di R [4].

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

# TEST DI NORMALITA' (GAUSSIANITA')

#

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

library(moments) # carica il pacchetto per asimmetria e curtosi

library(nortest) # carica il pacchetto per ulteriori test di normalità

#

mydata <- as.matrix(ais[c(5)]) # ci interessa la sola colonna con la ferritina

#

options(scipen=999) # esprime i numeri in formato fisso

#

agostino.test(mydata) # test di D'Agostino per il coefficiente di asimmetria

anscombe.test(mydata) # test di Anscombe-Glynn per il coefficiente di curtosi

#

ad.test(mydata) # test di Anderson-Darling per la normalità

cvm.test(mydata) # test di Cramer-von Mises per la normalità

pearson.test(mydata) # test chi-quadrato di Pearson per la normalità

sf.test(mydata) # test di Shapiro-Francia per la normalità

#

options(scipen=0) # ripristina la notazione scientifica

#

Dopo avere installato i tre pacchetti necessari, selezionato la colonna 5 contenente i dati che ci interessano, imposto di riportare i numeri in formato fisso con scipen=999, sono eseguiti il test di asimmetria, il test di curtosi e altri quattro test di normalità per la variabile ferritina [5]:

> agostino.test(mydata) # test di D'Agostino per il coefficiente di asimmetria

        D'Agostino skewness test

data:  mydata

skew = 1.2806, z = 6.1376, p-value = 0.0000000008377

alternative hypothesis: data have a skewness

> anscombe.test(mydata) # test di Anscombe-Glynn per il coefficiente di curtosi

        Anscombe-Glynn kurtosis test

data:  mydata

kurt = 4.4202, z = 2.9449, p-value = 0.003231

alternative hypothesis: kurtosis is not equal to 3

> ad.test(mydata) # test di Anderson-Darling per la normalità

        Anderson-Darling normality test

data:  mydata

A = 6.097, p-value = 0.000000000000004904

> cvm.test(mydata) # test di Cramer-von Mises per la normalità

        Cramer-von Mises normality test

data:  mydata

W = 0.94473, p-value = 0.000000002495

> pearson.test(mydata) # test chi-quadrato di Pearson per la normalità

        Pearson chi-square normality test

data:  mydata

P = 63.772, p-value = 0.00000002531

> sf.test(mydata) # test di Shapiro-Francia per la normalità

        Shapiro-Francia normality test

data:  mydata

W = 0.89138, p-value = 0.000000001251

Tutti i test di normalità concordano - con valori di p sempre di molto inferiori al valore soglia tradizionalmente impiegato di p=0.05 - nell'attribuire ai valori della ferritina una distribuzione non gaussiana [6].

Se nello script sostituite ais[c(5)] con ais[c(10)] potete eseguire i test di normalità sulla variabile altezza (espressa in centimetri, cm), ottenendo questi risultati:

> agostino.test(mydata) # test di D'Agostino per il coefficiente di asimmetria

        D'Agostino skewness test

data:  mydata

skew = -0.1993, z = -1.1856, p-value = 0.2358

alternative hypothesis: data have a skewness

> anscombe.test(mydata) # test di Anscombe-Glynn per il coefficiente di curtosi

        Anscombe-Glynn kurtosis test

data:  mydata

kurt = 3.5281, z = 1.5412, p-value = 0.1233

alternative hypothesis: kurtosis is not equal to 3

> ad.test(mydata) # test di Anderson-Darling per la normalità

        Anderson-Darling normality test

data:  mydata

A = 0.4462, p-value = 0.2793

> cvm.test(mydata) # test di Cramer-von Mises per la normalità

        Cramer-von Mises normality test

data:  mydata

W = 0.058974, p-value = 0.3886

> pearson.test(mydata) # test chi-quadrato di Pearson per la normalità

        Pearson chi-square normality test

data:  mydata

P = 17.653, p-value = 0.223

> sf.test(mydata) # test di Shapiro-Francia per la normalità

        Shapiro-Francia normality test

data:  mydata

W = 0.98919, p-value = 0.1174

I risultati, con valori sempre abbondantemente superiori al p=0.05, ci dicono che questa volta siamo di fronte a dati distribuiti in modo gaussiano.

A conferma dei valori numerici forniti dai test di normalità aggiungiamo ora quattro grafici, partendo da quello della ferritina. Copiate questo script, incollatelo nella Console di R e premete ↵ Invio:

# VALUTAZIONE GRAFICA DELLA NORMALITA' (GAUSSIANITA')

#

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

#

mydata <- as.matrix(ais[c(5)]) # ci interessa la sola colonna con la ferritina

#

windows() # apre una nuova finestra

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

#

# istogramma e kernel density plot

#

hist(mydata, xlim=c(0, 250), ylim=c(0, 0.020), freq=FALSE, breaks=25, 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=1000, from=0, to=250), xlim=c(0,250), ylim=c(0, 0.020), 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=1000) # 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

#

Con la funzione par(mfrow=c(2,2)) la finestra grafica viene suddivisa in quattro quadranti, che verranno riempiti da sinistra a destra e dall'alto verso il basso dai quattro grafici generati con i successivi blocchi di codice, migliorando la sintesi dei risultati.

Per le funzioni qui riportate [7] sono previsti molti altri possibili argomenti, per i quali si rimanda alla documentazione di R. Quelli che vi capiterà di utilizzare più frequentemente, a parte ovviamente la variabile da analizzare, che è sempre il primo argomento in ogni funzione, sono:

→ main, il titolo del grafico, che compare in alto;

→ xlab, l'etichetta dell'asse delle ascisse x;

→ ylab, l'etichetta dell'asse delle ordinate y;

→ xlim, che indica limite inferiore e limite superiore della scala da applicare all'asse delle ascisse x;

→ ylim, che indica limite inferiore e limite superiore della scala da applicare all'asse delle ordinate y;

→ xaxp, che specifica il limite inferiore, il limite superiore e il numero di intervalli da applicare alla scala dell'asse delle ascisse x;

→ yaxp, che specifica il limite inferiore, il limite superiore e il numero di intervalli da applicare alla scala dell'asse delle ordinate y.

Questi sono i quattro grafici risultanti, con la distribuzione dei dati campionari (in colore nero) che si discosta dalla distribuzione gaussiana teorica prevista (in colore rosso) confermando anche visivamente che i valori di concentrazione delle ferritina nei 202 atleti australiani non sono distribuiti in modo gaussiano.

Se nello script sostituite ais[c(5)] con ais[c(10)] e aggiustate le scale degli assi potete visualizzare il grafico dei valori della variabile altezza.

Questi sono i quattro grafici risultanti, che mostrano come la distribuzione dei dati campionari (in colore nero) segue la distribuzione gaussiana teorica prevista (in colore rosso), illustrando graficamente quanto avevamo già visto con i test di normalità: i valori di altezza nei 202 atleti australiani sono distribuiti in modo gaussiano.

Conclusione:

→ nel caso della ferritina i dati non sono distribuiti in modo gaussiano, pertanto dovremo impiegare metodi non parametrici e calcoleremo la mediana, la deviazione assoluta mediana o MAD e i quantili non parametrici;

→ nel caso dell'altezza i dati sono distribuiti in modo gaussiano, pertanto potremo impiegare metodi parametrici e calcoleremo la media, la deviazione standard e i quantili parametrici.

----------

[1] Ghasemi A, Zahediasl S. Normality Tests for Statistical Analysis: A Guide for Non-Statisticians. Int J Endocrinol Metab. 2012;10(2);486–489.

https://pmc.ncbi.nlm.nih.gov/articles/PMC3693611/

[2] Per il quadro tipico di una distribuzione gaussiana vedere il post: La distribuzione gaussiana

[3] Vedere i dati ematologici e i dati biometrici rilevati in 202 atleti australiani nel post Il set di dati ais nel quale trovate anche come caricare i dati della tabella senza impiegare il pacchetto DAAG. 

[4] Available CRAN Packages By Name.

https://cran.r-project.org/web/packages/available_packages_by_name.html

[5] 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. La concentrazione è espressa in µg/L.

[6] Se si vogliono avere i risultati di questi test di normalità sotto forma di una tabella compatta, vedere il post Tabulare una serie di test di normalità (gaussianità)

[7] Per la documentazione delle funzioni digitare help(nomedellafunzione) nella Console di R. Per gli argomenti xlog, ylog, xlim, ylim, xaxp, yaxp vedere anche la funzione par() con help(par).