Tabulare una serie di test di normalità (gaussianità)

Quando con R si esegue una serie di test di normalità (gaussianità) su una variabile, può essere fastidioso consultare e confrontare tra loro risultati che nella Console di R sono riportati separatamente e lontani l'uno dall'altro. Personalmente trovo che gli stessi risultati compattati in una tabella e incolonnati siano molto meglio fruibili. Ma non solo. I valori di probabilità p incolonnati sono poi ancor meglio confrontabili tra loro se riportati in formato fisso anziché in formato esponenziale.    

Vediamo quindi come incolonnare e riportare in formato fisso i risultati dei più comuni test di gaussianità con lo script che segue. Accertatevi di avere installati i pacchetti moments e nortest. Altrimenti dovete scaricarli e installarli dal CRAN. Accertatevi inoltre di avere installato il pacchetto DAAG che contiene i dati impiegati nello script, o in alternativa procedete come indicato [1].

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

# TABULARE UNA SERIE DI TEST DI NORMALITA' (GAUSSIANITA')

#

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

mydata <- ais[,c(5)] # salva in mydata la colonna 5 con la variabile ferritina

#

# asimmetria, curtosi e quattro test di gaussianità

#

library(moments) # carica il pacchetto

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

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

library(nortest) # carica il pacchetto

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

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

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

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

#

# sintetizza in una tabella i test con i risultati e i valori di p

#

stat <- c(names(agostino.test(mydata)$statistic[2]), names(anscombe.test(mydata)$statistic[2]), names(ad.test(mydata)$statistic), names(cvm.test(mydata)$statistic), names(pearson.test(mydata)$statistic), names(sf.test(mydata)$statistic)) # array con le denominazioni delle statistiche calcolate

#

ris <- c(agostino.test(mydata)$statistic[2], anscombe.test(mydata)$statistic[2], ad.test(mydata)$statistic, cvm.test(mydata)$statistic, pearson.test(mydata)$statistic, sf.test(mydata)$statistic) # array con i risultati delle statistiche calcolate

#

pval <- c(agostino.test(mydata)$p.value, anscombe.test(mydata)$p.value, ad.test(mydata)$p.value, cvm.test(mydata)$p.value, pearson.test(mydata)$p.value, sf.test(mydata)$p.value) # array con i valori di p

#

mydataset <- data.frame(stat, ris, pval) # combina gli array in una tabella

#

rownames(mydataset) <- c(agostino.test(mydata)$method, anscombe.test(mydata)$method, ad.test(mydata)$method, cvm.test(mydata)$method, pearson.test(mydata)$method, sf.test(mydata)$method) # aggiunge i nomi delle righe

colnames(mydataset) <- c("Statistica", "Risultato", "Valore p") # aggiunge i nomi delle colonne

mydataset # mostra la tabella

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

mydataset # mostra nuovamente la tabella

options(scipen=0) # ripristina la notazione scientifica/esponenziale

#

Il test di asimmetria di D'Agostino agostino.test() valuta quanto e come la simmetria della distribuzione campionaria (quella effettivamente osservata) si discosta da quella della distribuzione gaussiana teorica, che è per definizione perfettamente simmetrica. Il test di curtosi di Anscombe-Glynn anscombe.test() valuta se la distribuzione campionaria è troppo appiattita (bassa e larga) o troppo appuntita (alta e stretta) rispetto alla distribuzione gaussiana teorica, in ciascun punto della quale vale un preciso rapporto tra distanza dal centro della distribuzione e altezza della curva.

Mentre i due test precedenti valutano separatamente asimmetria e curtosi, gli altri quattro test qui impiegati, il test di Anderson-Darling ad.test(), il test di Cramer-von Mises cvm.test(), il test chi-quadrato di Pearson pearson.test() e il test di Shapiro-Francia sf.test(), sono test globali di gaussianità (forniscono una valutazione globale del grado di scostamento della distribuzione osservata dalla distribuzione gaussiana teorica, ma non del tipo di scostamento).

Dopo avere calcolato e riportato separatamente i risultati dei sei test, viene generata la tabella che mostra per ciascun test la denominazione, il simbolo della statistica calcolata, il risultato numerico di detta statistica e il corrispondente valore di probabilità p procedendo in questo modo:

→ per ciascun test viene riportato nell'array (o se preferite vettore) stat il simbolo della statistica calcolata names();

→ per ciascun test viene riportato nell'array ris il risultato numerico della statistica calcolata $statistic;

→ per ciascun test viene riportato nell'array pval il valore di probabilità p del risultato numerico $p.value;

→ mediante la funzione data.frame() gli array stat, ris e pval sono combinati nella tabella mydataset;
→ mediante la funzione rownames() a ciascuna riga della tabella mydataset viene assegnata la denominazione del test corrispondente $method;

→ mediante la funzione colnames() sono assegnati i nomi alle colonne della tabella mydataset che contengono, nell'ordine, il simbolo stat della statistica (Statistica), il risultato ris della statistica (Risultato) e il valore pval di probabilità p (Valore p) [2].

La tabella viene quindi mostrata nella Console di R prima con i valori numerici espressi in notazione scientifica

                                  Statistica  Risultato     Valore p

D'Agostino skewness test                   z  6.1376120 8.377114e-10

Anscombe-Glynn kurtosis test               z  2.9448982 3.230609e-03

Anderson-Darling normality test            A  6.0970351 4.903750e-15

Cramer-von Mises normality test            W  0.9447263 2.495032e-09

Pearson chi-square normality test          P 63.7722772 2.530969e-08

Shapiro-Francia normality test             W  0.8913765 1.251172e-09

e poi con i valori numerici espressi in notazione fissa mediante options(scipen=999)

                                  Statistica  Risultato               Valore p

D'Agostino skewness test                   z  6.1376120 0.00000000083771142715

Anscombe-Glynn kurtosis test               z  2.9448982 0.00323060945215406870

Anderson-Darling normality test            A  6.0970351 0.00000000000000490375

Cramer-von Mises normality test            W  0.9447263 0.00000000249503195167

Pearson chi-square normality test          P 63.7722772 0.00000002530968668914

Shapiro-Francia normality test             W  0.8913765 0.00000000125117173937

L'incolonnamento e il formato fisso rendono sicuramente più chiara la sintesi e la valutazione comparativa dei risultati numerici. La probabilità p è la probabilità di osservare per caso il risultato numerico della statistica (z, A, W, P secondo i casi) che misura lo scostamento della distribuzione osservata dalla distribuzione gaussiana teorica. Se tale probabilità è sufficientemente piccola, si conclude che il risultato ottenuto non è imputabile al caso e che pertanto la distribuzione dei dati si discosta significativamente dalla distribuzione gaussiana. Tutti e sei i test confermano che la ferritina non segue una distribuzione gaussiana: infatti abbiamo un valore p < 0.05 - cioè un p inferiore al valore del 5% tradizionalmente assunto come soglia per la significatività - in tutti i casi, e addirittura di molti ordini di grandezza inferiore al 5% in cinque casi su sei. Potete avere una conferma visiva di queste conclusioni con le rappresentazioni grafiche, sempre altamente consigliate come complemento all'analisi statistica, che potete realizzare con lo script riportato nel post Analizzare graficamente la distribuzione di una variabile. 

Come al solito lo script è realizzato in modo da essere facilmente riutilizzabile, essendo sufficiente a questo scopo:

→ sostituire la prima riga di codice con il codice per importare i vostri dati [3];

→ sostituire nella seconda riga di codice ais[,c(5)] con il nome della variabile da analizzare contenuta nei dati che avete importato.

----------

[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] Potete al bisogno approfondire la documentazione delle funzioni qui impiegate digitando help(nomedellafunzione) nella Console di R.

[3] Per una guida rapida all'importazione dei dati potete consultare i post:

→ importazione dei dati da un file .csv

→ importazione dei dati da un file .xls o .xlsx

→ gestione dei dati mancanti

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).

Valutare asimmetria e curtosi

I due test di normalità (gaussianità) più classici e tradizionali sono:

→ il test di asimmetria (test di D'Agostino) che valuta il grado di scostamento della simmetria della distribuzione campionaria (quella effettivamente osservata) da quella della distribuzione gaussiana teorica, che è per definizione perfettamente simmetrica;

→ il test di curtosi (test di Anscombe-Glynn) che ci dice se la distribuzione campionaria è troppo appiattita (bassa e larga) o troppo appuntita (alta e stretta) rispetto alla distribuzione gaussiana teorica, in ciascun punto della quale vale un preciso rapporto tra larghezza e altezza.

Prima di eseguire lo script è necessario scaricare e installare il pacchetto aggiuntivo moments. Trovate il manuale di riferimento, che include i riferimenti bibliografici, sul repository della documentazione di R [1]. Per il test di D'Agostino vedere anche [2].

I dati sono contenuti nella tabella ais del pacchetto DAAG - accertatevi di avere installato il pacchetto o in alternativa procedete come indicato in [3].

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

#  TEST DI ASIMMETRIA E TEST DI CURTOSI

#

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

library(moments) # carica il pacchetto per l'analisi statistica

#

skewness(ais$ferr) # coefficiente di asimmetria

agostino.test(ais$ferr) # test di D'Agostino per il coefficiente di asimmetria

kurtosis(ais$ferr) # coefficiente di curtosi

anscombe.test(ais$ferr) # test di Anscombe-Glynn per il coefficiente di curtosi

#

Dopo avere caricato il pacchetto DAAG che contiene i dati nella tabella ais  e il pacchetto moments che contiene le funzioni che ci interessano, nella terza riga di codice con la funzione skewness() viene calcolato il coefficiente di asimmetria della variabile ais$ferr (concentrazione della ferritina nel siero espressa in µg/L) la cui significatività viene valutata nella riga successiva con la funzione agostino.test().

Il tutto viene ripetuto per il coefficiente di curtosi, calcolato con la funzione kurtosis() e la cui significatività viene poi valutata nella riga successiva mediante la funzione anscombe.test().

Il coefficiente di asimmetria è significativo, con un valore osservato di asimmetria skew=1.2806 che dista z=6.1376 deviazioni standard dal valore 0 previsto per una distribuzione gaussiana, e una probabilità di osservare per caso tale distanza molto bassa (p-value=8.377e-10):

        D'Agostino skewness test

data:  ais$ferr

skew = 1.2806, z = 6.1376, p-value = 8.377e-10

alternative hypothesis: data have a skewness

Il coefficiente di curtosi è anch'esso significativo, con un valore osservato di curtosi kurt=4.4202 che dista z=2.9449 deviazioni standard dal valore 3 previsto per una distribuzione gaussiana, e una probabilità di osservare per caso tale distanza bassa (p-value=0.003231):

        Anscombe-Glynn kurtosis test

data:  ais$ferr

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

alternative hypothesis: kurtosis is not equal to 3

I criteri per valutare la significatività statistica sono ovviamente i soliti:

→ si  considera il risultato non significativo quando la probabilità di osservare per caso il valore (effettivamente ottenuto nel caso specifico) della statistica è "elevato";

→ si  considera il risultato significativo quando la probabilità di osservare per caso il valore (effettivamente ottenuto nel caso specifico) della statistica è "sufficientemente basso";

→ come valore soglia tra "elevato" e "sufficientemente basso" viene tradizionalmente impiegato il valore di probabilità p=0.05 (ma nulla vieta, al fine di ridurre ulteriormente la possibilità di trarre dal test statistico conclusioni errate, di impiegare un valore soglia inferiore, come per esempio p=0.01 e quindi considerare significativo un risultato statistico solamente quando la probabilità di osservarlo per caso è dell'1% o meno).

Integriamo ora l'analisi statistica con una rappresentazione grafica dei dati che riporta, sovrapposti, l'istogramma e il kernel density plot della distribuzione campionaria e per confronto la distribuzione gaussiana teorica, quella che gli stessi dati avrebbero se fossero distribuiti in modo gaussiano.

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

# DISTRIBUZIONE CAMPIONARIA vs. GAUSSIANA TEORICA

#

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$ferr, xlim=c(0,250), ylim=c(0,50), freq=TRUE, breaks="FD", density=20, angle=45, border="black", main="Distribuzione campionaria vs. gaussiana teorica", xlab="Ferritina in µg/L", ylab="Frequenza (numero dei casi)") # traccia l'istogramma 

#

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

plot(density(ais$ferr), yaxt="n", xlim=c(0,250), ylim=c(0,0.012), col="blue", lwd=1) # sovrappone all'istogramma il kernel density plot

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

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(new=TRUE, ann=FALSE) # consente la sovrapposizione del grafico successivo

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

y <- dnorm(x, mean=mean(ais$ferr), sd=sd(ais$ferr)) # 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", lty=2, lwd=2) # sovrappone la distribuzione gaussiana teorica in colore rosso

#

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

#

Questo è il grafico risultante:

La conclusione dei test statistici trova una conferma empirica nella rappresentazione grafica dei dati: la concentrazione della ferritina nel siero, rispetto alla corrispondente distribuzione gaussiana teorica (in tratteggio di colore rosso) – cioè rispetto a quella che i dati dovrebbero avere se fossero distribuiti in modo gaussiano – è eccessivamente appuntita (leptocurtica) e asimmetrica (asimmetria positiva) pertanto la distribuzione dei dati non è gaussiana.

Nel codice riportato si possono notare alcune cose:

→ il numero delle classi breaks="FD" è calcolato con la regola di Friedman-Diaconis [4];

→ le barre dell'istogramma sono profilate in nero (border="black") riempite di colore (il colore grigio di default) con un tratteggio che copre il 20% della superficie (density=20) ed è inclinato a 45 gradi (angle=45);

→ la scala delle ordinate di sinistra riporta, per l'istogramma, il numero dei casi osservati in ciascuna classe;

→ la scala delle ordinate di destra riporta, per il kernel density plot, la stima kernel della densità di probabilità [5];

→ per la distribuzione gaussiana teorica vale la stessa scala di densità di probabilità riportata sulla destra;

→ le ordinate della distribuzione gaussiana corrispondente ai dati sono calcolate con la funzione dnorm() impiegando la media (mean(ais$ferr)) e la deviazione standard (sd(ais$ferr)) campionarie della ferritina;

→ la distribuzione gaussiana teorica, quella che i dati avrebbero avere, è riportata con una linea tratteggiata (lty=2) di spessore aumentato (lwd=2) in colore rosso (col="red").

Nel post Analizzare graficamente la distribuzione di una variabile potete trovare, riunite in un'unica vista, quattro grafici per una rappresentazione più completa delle differenze tra la vostra distribuzione campionaria e la distribuzione gaussiana teorica.

I test di asimmetria e di curtosi sono molto importanti, ma più importante ancora è che facciano parte di un approccio complessivo e organico al calcolo delle statistiche elementari di una variabile che, si ricorda, prevede di:

→ effettuare una analisi esplorativa dei dati

→ effettuare uno o più test di normalità (gaussianità).

→ impiegare le statistiche elementari parametriche se i dati sono distribuiti in modo gaussiano;

→ impiegare le statistiche elementari non parametriche se i dati non sono distribuiti in modo gaussiano.

----------

[1] Available CRAN Packages By Name.

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

[2] D'Agostino RB. Transformation to Normality of the Null Distribution of g1. Biometrika 1970;57(3);679-681.

https://www.jstor.org/stable/2334794#page_scan_tab_contents

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

[4] Vedere: Freedman–Diaconis rule.

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

[5] Per i dettagli della rappresentazione grafica vedere il post Istogrammi e il post Kernel density plot.