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