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 i dati campionari che confermano 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 per i valori per la variabile altezza.

I grafici confermano 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. URL consultato il 05/01/2019: https://goo.gl/AkoUdn

[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. URL consultato il 05/01/2019: https://goo.gl/hLC9BB

[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] Quei stessi stessi test di normalità sono riportati, in una tabella più compatta, nel 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).