Statistiche elementari parametriche

Siamo al terzo passo delle valutazioni che è sempre necessario effettuare nell'ambito di un percorso logico che prevede, per il calcolo delle statistiche elementari di una singola variabile (analisi univariata):

→ analisi esplorativa dei dati;

→ esecuzione dei test di normalità (gaussianità) per valutare se i dati seguono una distribuzione gaussiana;

→ calcolo delle statistiche elementari parametriche (media, deviazione standard, varianza, quantili parametrici) se i dati seguono una distribuzione gaussiana;

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

Si tratta del passo conclusivo per la variabile altezza (espressa in cm) rilevata in 202 atleti australiani, in quanto con i test di normalità abbiamo stabilito che la variabile è distribuita in modo gaussiano: questa era la sintesi grafica dei risultati.

Ora copiate questo script, incollatelo nella Console di R e premete ↵ Invio:

# STATISTICHE ELEMENTARI PARAMETRICHE 

#

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

mydata <- as.matrix(ais[c(10)]) # ci interessa la sola colonna con l'altezza

#

# statistiche elementari parametriche 

#

min(mydata) # valore minimo

mean(mydata) # media 

max(mydata) # valore massimo

range(mydata) # range 

max(mydata)-min(mydata) # campo di variazione

var(mydata) # varianza

sd(mydata) # deviazione standard

(sd(mydata)/mean(mydata))*100 # coefficiente di variazione CV %

#

# confronta media e mediana

#

media <- mean(mydata) # calcola la media

mediana <- median(mydata) # calcola la mediana

data.frame(media, mediana) # le mette a confronto

#

# confronta deviazione standard e MAD

#

ds <- sd(mydata) # calcola la deviazione standard

mad <- mad(mydata) # calcola la Median Absolute Deviation (about the median) o MAD

data.frame(ds, mad) # le mette a confronto

#

# confronta quantili parametrici e non parametrici

#

qpar <- round(qnorm(c(seq(0.025, 0.975, 0.025)), mean=mean(mydata), sd=sd(mydata)), digits=2) # calcola i quantili parametrici

qnon <- round(quantile(mydata, probs=seq(0.025, 0.975, 0.025)), digits=2) # calcola i quantili non parametrici

dper <- abs(round(100*(qpar-qnon)/((qpar+qnon)/2), digits=2)) # calcola la differenza in percentuale

data.frame(qpar, qnon, dper) # li mette a confronto

#

Non avrebbe senso ripetere l'analisi esplorativa dei dati, i test di normalità e la sintesi grafica dei risultati - per questi si rimanda agli script riportati nel percorso logico indicato all'inizio. Per cui qui ci limitiamo a calcolare le usuali statistiche elementari parametriche, riportando per confronto i risultati delle statistiche elementari non parametriche corrispondenti. 

Dopo avere caricato il pacchetto DAAG [1] con la tabella ais, la variabile che ci interessa, contenuta nella colonna 10 della tabella ais viene assegnata (<-) all'oggetto mydata, sul quale sono calcolate le statistiche elementari parametriche con le funzioni:

→ min() per riportare il valore minimo osservato;

→ mean() per calcolare la media;

→ max() per riportare il valore massimo osservato;

→ range() per calcolare il range, che qui è correttamente;

→ max(mydata)-min(mydata) per calcolare il campo di variazione, talora indicato (impropriamente) come "range";

→ var() per calcolare la varianza,

→ sd() per calcolare la deviazione standard;

→ e con l'espressione (sd(mydata)/mean(mydata))*100 che calcola il coefficiente di variazione percentuale.

> min(mydata) # valore minimo

[1] 148.9

> mean(mydata) # media 

[1] 180.104

> max(mydata) # valore massimo

[1] 209.4

> range(mydata) # range 

[1] 148.9 209.4

> max(mydata)-min(mydata) # campo di variazione

[1] 60.5

> var(mydata) # varianza

         ht

ht 94.76038

> sd(mydata) # deviazione standard

[1] 9.734494

> (sd(mydata)/mean(mydata))*100 # coefficiente di variazione CV %

[1] 5.404931

Sono poi messe a confronto la media con la mediana

> media <- mean(mydata) # calcola la media

> mediana <- median(mydata) # calcola la mediana

> data.frame(media, mediana) # le mette a confronto

    media mediana

1 180.104   179.7

la deviazione standard con la MAD [2]

> ds <- sd(mydata) # calcola la deviazione standard

> mad <- mad(mydata) # calcola la Median Absolute Deviation (about the median) o MAD

> data.frame(ds, mad) # le mette a confronto

        ds     mad

1 9.734494 8.96973

e i quantili parametrici con i quantili non parametrici, riportando nell'ultima colonna la differenza percentuale tra i due

> qpar <- round(qnorm(c(seq(0.025, 0.975, 0.025)), mean=mean(mydata), sd=sd(mydata)), digits=2) # calcola i quantili parametrici

> qnon <- round(quantile(mydata, probs=seq(0.025, 0.975, 0.025)), digits=2) # calcola i quantili non parametrici

> dper <- abs(round(100*(qpar-qnon)/((qpar+qnon)/2), digits=2)) # calcola la differenza in percentuale

> data.frame(qpar, qnon, dper) # li mette a confronto

        qpar   qnon dper

2.5%  161.02 158.98 1.28

5%    164.09 163.96 0.08

7.5%  166.09 167.35 0.76

10%   167.63 169.17 0.91

12.5% 168.91 170.36 0.85

15%   170.01 171.40 0.81

17.5% 171.01 172.22 0.71

20%   171.91 172.76 0.49

22.5% 172.75 173.52 0.44

25%   173.54 174.00 0.26

27.5% 174.29 174.18 0.06

30%   175.00 175.00 0.00

32.5% 175.69 175.73 0.02

35%   176.35 176.07 0.16

37.5% 177.00 177.30 0.17

40%   177.64 177.94 0.17

42.5% 178.26 178.44 0.10

45%   178.88 178.99 0.06

47.5% 179.49 179.60 0.06

50%   180.10 179.70 0.22

52.5% 180.71 180.10 0.34

55%   181.33 180.36 0.54

57.5% 181.94 181.00 0.52

60%   182.57 182.66 0.05

62.5% 183.21 183.06 0.08

65%   183.85 183.90 0.03

67.5% 184.52 184.67 0.08

70%   185.21 185.17 0.02

72.5% 185.92 185.60 0.17

75%   186.67 186.18 0.26

77.5% 187.46 187.28 0.10

80%   188.30 188.62 0.17

82.5% 189.20 189.18 0.01

85%   190.19 190.67 0.25

87.5% 191.30 191.94 0.33

90%   192.58 192.79 0.11

92.5% 194.12 193.86 0.13

95%   196.12 195.17 0.49

97.5% 199.18 197.48 0.86

Le differenze trascurabili tra statistiche elementari parametriche e statistiche elementari non parametriche riflettono il fatto che nel caso di una distribuzione perfettamente gaussiana [3] le due coincidono - è quindi assolutamente lecito impiegare le seconde al posto delle prime. Mentre il contrario è sbagliato - nel caso di una distribuzione non gaussiana le statistiche elementari parametriche non devono essere impiegate.

Lo script può essere riutilizzato molto semplicemente, assegnando (<-) all'oggetto mydata i propri dati.

----------

[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] La Median Absolute Deviation about median o MAD ovvero la deviazione assoluta mediana (dalla mediana) è l'equivalente non parametrico della deviazione standard. Vedere: Rousseeuw PJ, Croux C. Alternatives to the Median Absolute Deviation. Journal of the American Statistical Association 88 (424), 1273-1283, 1993.

https://www.jstor.org/stable/2291267?seq=1#page_scan_tab_contents

[3] Un esempio di dati distribuiti in modo perfettamente gaussiano è riportato nel post La distribuzione gaussiana.

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.