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

Calcolare la media cumulativa e la mediana cumulativa

Data una sequenza di valori numerici la media cumulativa è la media, calcolata in corrispondenza di ciascun dato, di tutti i dati dal primo fino a quello corrente.

Così ad esempio in corrispondenza del dato numero 5 la media cumulativa è la media dei dati da 1 a 5, in corrispondenza del dato numero 6 la media cumulativa è la media dei dati da 1 a 6, in corrispondenza del dato numero 7 la media cumulativa è la media dei dati da 1 a 7 e così via. Lo stesso principio vale per la mediana cumulativa che a sua volta è la mediana, calcolata in corrispondenza di ciascun dato, di tutti i dati dal primo fino a quello corrente.

Per il calcolo della media cumulativa e della mediana cumulativa impieghiamo le funzioni contenute nel pacchetto cumstats, che deve essere scaricato dal CRAN.

Copiate e incollate questo script nella Console di R e premete ↵ Invio:

# CALCOLARE LA MEDIA CUMULATIVA E LA MEDIANA CUMULATIVA 

#

library(cumstats) # carica il pacchetto

rileva_x <- seq(1:50) # numero progressivo della rilevazione in ascisse

valore_y <- c(5, 9, 12, 7, 8, 13, 26, 7, 8, 19, 10, 7, 5, 36, 12, 18, 31, 17, 16, 43, 7, 12, 9, 27, 15, 18, 24, 29, 32, 23, 29, 12, 7, 48, 34, 45, 19, 12, 16, 18, 22, 31, 38, 7, 11, 4, 38, 36 ,43, 15) # valore osservato corrispondente in ordinate

mydata <- data.frame(rileva_x, valore_y) # combina i vettori in una matrice

#

mydata$mediacum <- cummean(mydata$valore_y) # calcola la media cumulativa e la aggiunge in una nuova colonna

mydata$medianacum <- cummedian(mydata$valore_y) # calcola la mediana cumulativa e la aggiunge in una nuova colonna

mydata # mostra la matrice con i dati definitivi

#

plot(mydata$rileva_x, mydata$valore_y, xlim=c(0,50), ylim=c(0,60), type="l", lty=1, lwd=1, col="black", main="Valori osservati, media cumulativa e mediana cumulativa", xlab="Rilevazioni", ylab="Valori osservati") # grafico dei dati originali

lines(mydata$rileva_x, mydata$mediacum, xlim=c(0,50), ylim=c(0,60), type="l", lty=2, lwd=1, col="red") # sovrappone grafico della media cumulativa

lines(mydata$rileva_x, mydata$medianacum, xlim=c(0,50), ylim=c(0,60), type="l", lty=2, lwd=1, col="blue") # sovrappone grafico della mediana cumulativa

#

legend(0,60, legend=c("Valori osservati nelle rilevazioni", "Media cumulativa dei valori osservati", "Mediana cumulativa dei valori osservati"), col=c("black", "red", "blue"), lty=c(1,2,2), lwd=c(1,1,1), cex=0.8) # aggiunge la legenda

#

Dopo avere caricato il pacchetto cumstats i dati impiegati come esempio sono generati in questo modo:

→ nel vettore rileva_x viene riportato il numero progressivo della rilevazione;

→ nel vettore valore_y viene riportato il valore numerico osservato in corrispondenza di ciascuna rilevazione, valore numerico sul quale calcoleremo la media cumulativa e la mediana cumulativa.

→ con la funzione dataframe() i due vettori sono combinati nella matrice mydata [1].

A questo punto alla matrice mydata sono aggiunte due colonne:

→ la colonna (vettore) mediacum nella quale è riportata la media cumulativa calcolata in corrispondenza di ogni dato con la funzione cummmean();

→ la colonna (vettore) medianacum nella quale è riportata la mediana cumulativa calcolata in corrispondenza di ogni dato con la funzione cummedian().

Per i dettagli delle funzioni impiegate digitare help(nomedellafunzione) nella Console di R. Al bisogno potete consultare il manuale del pacchetto cumstats [2].

Non resta quindi che mostrare il contenuto della matrice mydata nella versione definitiva contenente il numero della rilevazione (rileva_x), il valore osservato (valore_y) e le rispettive media cumulativa (mediacum) e mediana cumulativa (medianacum)

> mydata # mostra la matrice con i dati definitivi

   rileva_x valore_y  mediacum medianacum

1         1        5  5.000000        5.0

2         2        9  7.000000        7.0

3         3       12  8.666667        9.0

4         4        7  8.250000        8.0

5         5        8  8.200000        8.0

6         6       13  9.000000        8.5

7         7       26 11.428571        9.0

8         8        7 10.875000        8.5

9         9        8 10.555556        8.0

10       10       19 11.400000        8.5

11       11       10 11.272727        9.0

12       12        7 10.916667        8.5

13       13        5 10.461538        8.0

14       14       36 12.285714        8.5

15       15       12 12.266667        9.0

16       16       18 12.625000        9.5

17       17       31 13.705882       10.0

18       18       17 13.888889       11.0

19       19       16 14.000000       12.0

20       20       43 15.450000       12.0

21       21        7 15.047619       12.0

22       22       12 14.909091       12.0

23       23        9 14.652174       12.0

24       24       27 15.166667       12.0

25       25       15 15.160000       12.0

26       26       18 15.269231       12.0

27       27       24 15.592593       12.0

28       28       29 16.071429       12.5

29       29       32 16.620690       13.0

30       30       23 16.833333       14.0

31       31       29 17.225806       15.0

32       32       12 17.062500       14.0

33       33        7 16.757576       13.0

34       34       48 17.676471       14.0

35       35       34 18.142857       15.0

36       36       45 18.888889       15.5

37       37       19 18.891892       16.0

38       38       12 18.710526       15.5

39       39       16 18.641026       16.0

40       40       18 18.625000       16.0

41       41       22 18.707317       16.0

42       42       31 19.000000       16.5

43       43       38 19.441860       17.0

44       44        7 19.159091       16.5

45       45       11 18.977778       16.0

46       46        4 18.652174       16.0

47       47       38 19.063830       16.0

48       48       36 19.416667       16.5

49       49       43 19.897959       17.0

50       50       15 19.800000       16.5

che sono impiegati per la successiva rappresentazione grafica

che viene realizzata con tre sole funzioni:

→ la funzione plot() che traccia il grafico del valore osservato mydata$valore_y in corrispondenza di ciascuna rilevazione mydata$rileva_x impiegando una linea (type="l") continua (lty=1) di larghezza unitaria (lwd=1) e di colore nero (col="black");

→ la funzione lines() che viene impiegata una prima volta per sovrapporre il grafico della media cumulativa con una linea tratteggiata (lty=2) di colore rosso (col="red") e una seconda volta per sovrapporre il grafico della mediana cumulativa con una linea tratteggiata (lty=2) di colore blu (col="blue") 

→ la funzione legend() che riporta la legenda esplocativa dei grafici [3].

Ovviamente lo script può essere riadattato per rappresentare dati differenti ma anche per riportare solamente l'uno o l'altro dei grafici.

----------

[1] Parliamo di array o vettore nel caso di dati numerici monodimensionali, disposti su una sola riga, 

8

6

11

7

di matrice nel caso di dati numerici disposti su più righe e più colonne

8	9	15	14
6	7	18	12
11	8	17	13
7	4	19	17

e di tabella nei casi in cui il contenuto, disposto su più righe e più colonne, è rappresentato oltre che da dati numerici, anche da testo e/o operatori logici

M	7	9	VERO
F	3	12	VERO
F	5	10	FALSO

[2] Vedere la voce Reference manual del pacchetto cumstats: Cumulative Descriptive Statistics. URL consultato il 02/01/2023.

[3] Per un approfondimento sull'impiego della funzione legend() rimando al post Aggiungere la legenda a un grafico.