venerdì 4 gennaio 2019

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 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 NameURL consultato il 04/01/2019: https://goo.gl/hLC9BB

[2] D'Agostino RB. Transformation to Normality of the Null Distribution of g1. Biometrika 1970;57(3);679-681. URL consultato il 05/01/2019: https://goo.gl/xMVUEx

[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 Wikipedia, the free encyclopedia. Freedman–Diaconis rule. URL consultato il 22/02/2023.

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


Nessun commento:

Posta un commento