La distribuzione gaussiana

Innazitutto una premessa. Uno dei punti più delicati nella statistica di base è rappresentato dalla distribuzione gaussiana (l'altro punto delicato è rappresentato dal coefficiente di correlazione r).

Media, deviazione standard, ANOVA, test t di Student, solo per citare i principali, sono test statistici (o più semplicemente statistiche) basati sull'assunto che i dati seguano una distribuzione gaussiana.

Poiché media e deviazione standard sono i “parametri” che descrivono una distribuzione gaussiana, i metodi statistici che sono basati sull'assunto che i dati siano distribuiti in modo gaussiano sono detti metodi parametrici, in contrapposizione a quelli che non prevedono questo assunto, che sono detti metodi non parametrici [1, 2].

Applicare metodi parametrici a variabili che non sono distribuite in modo gaussiano è inappropriato e porta a conclusioni sbagliate. Vediamo un esempio impiegando i valori della concentrazione della ferritina nel sangue (espressa in µg/L) rilevati in 202 atleti australiani.

Copiate il testo che segue, incollatelo nella Console di R e premete ↵ Invio:

x <- c(60, 68, 21, 69, 29, 42, 73, 44, 41, 44, 38, 26, 30, 48, 30, 29, 43, 34, 53, 59, 43, 40, 92, 48, 77, 71, 37, 71, 73, 85, 64, 19, 39, 41, 13, 20, 59, 22, 30, 78, 21, 109, 102, 71, 64, 68, 78, 107, 39, 58, 127, 102, 86, 40, 50, 58, 33, 51, 82, 25, 86, 22, 30, 27, 60, 115, 124, 54, 29, 164, 50, 36, 40, 62, 90, 12, 36, 45, 43, 51, 22, 44, 26, 16, 58, 46, 43, 34, 41, 57, 73, 88, 182, 80, 88, 109, 69, 41, 66, 63, 34, 97, 55, 76, 93, 46, 155, 99, 35, 124, 176, 107, 177, 130, 64, 109, 125, 150, 115, 89, 93, 183, 84, 70, 118, 61, 72, 91, 58, 97, 110, 72, 36, 53, 72, 39, 61, 49, 35, 8, 106, 32, 41, 20, 44, 103, 50, 41, 101, 73, 56, 74, 58, 87, 139, 82, 87, 80, 67, 132, 43, 212, 94, 213, 122, 76, 53, 91, 36, 101, 184, 44, 66, 220, 191, 40, 189, 141, 212, 97, 53, 126, 234, 50, 94, 156, 124, 87, 97, 102, 55, 117, 52, 133, 214, 143, 65, 90, 38, 122, 233, 32)

Ora calcoliamo su questi dati la classica misura di posizione, la media, e la classica misura di dispersione dei dati, la deviazione standard (ds). Digitate o copiate e incollate nella Console di R queste due righe e premete  ↵ Invio:

mean(x) # calcola la media della ferritina

sd(x) # calcola la deviazione standard della ferritina

Questo è il risultato:

> mean(x) # calcola la media della ferritina

[1] 76.87624

> sd(x) # calcola la deviazione standard della ferritina

[1] 47.50124

Comunicate i risultati di questa ricerca a un vostro amico, il quale vi fa notare quanto segue:

→ in base alle proprietà della distribuzione gaussiana sappiamo che tra la media - 1.96 · ds e la media + 1.96 · ds si deve trovare il 95% dei dati campionari;

→ il valore della media della ferritina (arrotondato) è 76.9 µg/L;

→ la deviazione standard della ferritina è 47.5 µg/L;

→ moltiplicando 47.5 µg/L per 1.96 otteniamo 93.1 µg/L.

→ dal fatto che 76.9 - 93.1 = -16.2 µg/L e 76.9 + 93.1 = 170 µg/L dovremmo dedurre che il 95% dei valori di ferritina rilevati negli atleti australiani cadeva nell’intervallo

-16.2 ÷ 170 µg/L

e pertanto concludere che in una certa quota degli atleti la concentrazione nel sangue della ferritina aveva valori inferiori a zero.

Voi ci restate malissimo, ma il punto è chiaro. Un evidente controsenso consente di individuare un grave errore nella metodologia statistica impiegata: non si possono applicare media e deviazione standard a dati che non sono distribuiti in modo gaussiano - e che questo sia proprio quello che accade nel caso della ferritina lo vedete seguendo il percorso logico indicato in [3].

Nel caso delle statistiche elementari, l’alternativa a media e deviazione standard è rappresentata dalle statistiche non parametriche mediana e deviazione assoluta mediana (MAD), dalla distanza interquartile e dai quantili non parametrici (quartili, decili, percentili).

Ma c'è un altro aspetto interessante che riguarda le statistiche non parametriche. Prendiamo questi dati

8

6

11

7

9

14

3

la cui media è 8.3, quindi ordiniamoli in ordine crescente

3

6

7

8

9

11

14

per calcolarne la mediana che è 8 (il valore centrale nella lista dei dati ordinati in ordine crescente).

Supponiamo ora di aggiungere una nuova osservazione che si discosta molto dalle precedenti ma di non avere ragioni plausibili per scartarla

8

6

11

7

9

14

3

45

Calcoliamo di nuovo la media, che ora è diventata 12.9, quindi di nuovo ordiniamo i dati in ordine crescente

3

6

7

8

9

11

14

45

per calcolarne la mediana, che ora è diventata (8 + 9) / 2 = 8.5.

Mentre con l'inserimento del nuovo dato la media è passata da 8.3 a 12.9, la mediana è passata da 8 a 8.5. Il fatto che la mediana sia una misura di posizione più robusta nei confronti di dati estremi rispetto alla media è un altro argomento a favore dell'impiego delle statistiche non parametriche [1, 2].

L'esempio della ferritina conferma la necessità che le statistiche parametriche siano impiegate solamente dopo che una attenta valutazione dei dati ha confermato che questi sono distribuiti in modo gaussiano.

Arriviamo quindi al tema. Mentre per la teoria della distribuzione gaussiana si rimanda ovviamente ai test di statistica, vediamo alcuni possibili modi per impiegare le funzioni di R nella valutazione dell'assunto che i dati siano distribuiti in modo gaussiano.

Le funzioni sono applicate a una distribuzione gaussiana, in modo da rapprsentare il quadro di riferimento da punto di vista numerico e da punto di vista grafico, quindi nell'ordine vediamo:

→ i test di normalità dei dati;

→ il confronto tra media e mediana (in un distribuzione gaussiana devono essere uguali);

→ il confronto tra deviazione standard e MAD (in un distribuzione gaussiana devono essere uguali);

→ il confronto tra percentili parametrici e non parametrici (in un distribuzione gaussiana devono essere uguali);

→ un istogramma e un kernel density plot dei dati campionari (in un distribuzione gaussiana devono essere simmetrici);

→ la distribuzione campionaria (deve sovrapporsi alla corrispondente curva gaussiana teorica);

→ la distribuzione cumulativa dei dati campionari (deve sovrapporsi alla funzione di distribuzione cumulativa teorica);

→ i quantili campionari (devono essere allineati sulla retta della distribuzione teorica).

Cinquemila valori di una distribuzione gaussiana con media 0 e deviazione standard 1 sono generati con la funzione rnorm() [4] e salvati nel file gauss.csv con questo script, copiatelo e incollatelo nella Console di R e premete ↵ Invio:

# SALVA I DATI DI UNA DISTRIBUZIONE GAUSSIANA IN UN FILE CSV

#

mydata <- rnorm(5000, mean = 0, sd = 1) # genera cinquemila valori con media=0, ds=1

#

write.table(mydata,file="C:/Rdati/gauss.csv",sep=";", dec=",", quote=FALSE, row.names=FALSE) # esporta i dati in un file di testo .csv sostituendo il punto (.) con la virgola (,)

#

Questo è quindi il nostro campione, che ci consente di illustrare cosa ci si deve attendere dall'analisi di una variabile se questa è distribuita in modo gaussiano.

Se non l'avete già fatto, prima di continuare dovete scaricare e installare dal CRAN il pacchetto moments e il pacchetto nortest.

I prossimi script impiegano i dati  salvati nel file gauss.csv [5].

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

# TEST DI NORMALITA' (GAUSSIANITA')

#

library(moments) # carica il pacchetto

library(nortest) # carica il pacchetto

data <- read.table("C:/Rdati/gauss.csv", header=TRUE, sep=";", dec=",") # importa i dati

mydata <- unlist(data) # li converte in formato numerico

#

# asimmetria, curtosi e altri quattro test di normalità

#

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 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à

#

Come atteso tutti i test di normalità confermano per i dati una distribuzione gaussiana [3]:

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

        D'Agostino skewness test

data:  mydata

skew = -0.0050374, z = -0.1456345, p-value = 0.8842

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.03034, z = 0.48405, p-value = 0.6283

alternative hypothesis: kurtosis is not equal to 3

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

        Anderson-Darling normality test

data:  mydata

A = 0.13271, p-value = 0.9808

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

        Cramer-von Mises normality test

data:  mydata

W = 0.019132, p-value = 0.9754

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

        Pearson chi-square normality test

data:  mydata

P = 47.14, p-value = 0.8452

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

        Shapiro-Francia normality test

data:  mydata

W = 0.99982, p-value = 0.9497

In una distribuzione gaussiana statistiche parametriche e statistiche non parametriche devono fornire gli stessi risultati. Copiate quindi questo script, incollatelo nella Console di R e premete ↵ Invio:

# TEST AGGIUNTIVI DI NORMALITA' (GAUSSIANITA') 

#

data <- read.table("C:/Rdati/gauss.csv", header=TRUE, sep=";", dec=",") # importa i dati

mydata <- unlist(data) # li converte in formato numerico

#

# 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=0, sd=1), 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

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

#

A meno di una differenza insignificante (dipendente dal generatore di numeri casuali di R) media e mediana sono uguali

       media      mediana

1 0.01061802 0.0009510756

deviazione standard e MAD [6] sono uguali

        ds       mad

1 1.000571 0.9948351

come pure i quantili parametrici e i quantili non parametrici

       qnon  qpar

2.5%  -1.97 -1.96

5%    -1.64 -1.64

7.5%  -1.42 -1.44

10%   -1.25 -1.28

12.5% -1.14 -1.15

15%   -1.02 -1.04

17.5% -0.92 -0.93

20%   -0.83 -0.84

22.5% -0.75 -0.76

25%   -0.66 -0.67

27.5% -0.59 -0.60

30%   -0.51 -0.52

32.5% -0.44 -0.45

35%   -0.37 -0.39

37.5% -0.31 -0.32

40%   -0.24 -0.25

42.5% -0.18 -0.19

45%   -0.12 -0.13

47.5% -0.05 -0.06

50%    0.00  0.00

52.5%  0.07  0.06

55%    0.13  0.13

57.5%  0.20  0.19

60%    0.26  0.25

62.5%  0.33  0.32

65%    0.41  0.39

67.5%  0.47  0.45

70%    0.54  0.52

72.5%  0.61  0.60

75%    0.68  0.67

77.5%  0.76  0.76

80%    0.84  0.84

82.5%  0.94  0.93

85%    1.06  1.04

87.5%  1.18  1.15

90%    1.31  1.28

92.5%  1.45  1.44

95%    1.64  1.64

97.5%  1.95  1.96

Ai valori numerici forniti dai test di normalità e dalle statistiche aggiuntive aggiungiamo ora quattro grafici. Copiate questo script, incollatelo nella Console di R e premete ↵ Invio:

# VALUTAZIONE GRAFICA DELLA NORMALITA' (GAUSSIANITA')

#

data <- read.table("C:/Rdati/gauss.csv", header=TRUE, sep=";", dec=",") # importa i dati

mydata <- unlist(data) # li converte in formato numerico

#

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(-4, 4), ylim=c(0, 0.5), freq=FALSE, breaks=50, 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=-4, to=4), xlim=c(-4, 4), ylim=c(0, 0.5), yaxt="n", col="black") # sovrappone il kernel density plot della distribuzione campionaria

#

# distribuzione campionaria vs. gaussiana teorica

#

plot(density(mydata, n=1024, from=-4, to=4), xlim=c(-4, 4), ylim=c(0, 0.5)) # traccia il kernel density plot della distribuzione campionaria

par(new=TRUE, ann=FALSE) # consente di sovrapporre il grafico successivo

x <- seq(-4, 4, 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(-4, 4), ylim=c(0, 0.5), 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

#

Questi sono i quattro grafici risultanti, con i dati campionari (in nero) che si sovrappongono perfettamente alla distribuzione teorica (in rosso), quello che aspettiamo appunto nel caso di dati distribuiti in modo gaussiano:

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.

Per riutilizzare gli script è sufficiente importare i vosti dati nell'oggetto mydata.

Conclusione: l'analisi statistica e grafica qui effettuata fornisce un esempio didattico di quello che tipicamente ci si deve attendere dall'analisi di una variabile se questa è distribuita in modo gaussiano.

----------

[1] Siegel S, Castellan NJ. Statistica non parametrica. McGraw-Hill Libri Italia, Milano, 1992, ISBN 88-386-0620-X. 

[2] Maritz JS. Distribution-free statistical methods. Chapman and Hall, New York, 1984, ISBN 0-412-25410-7. 

[3] Il percorso logico prevede, per il calcolo delle statistiche elementari di una singola variabile, i seguenti passi:

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

[4] Digitate help(rnorm) nella Console di R per la documentazione della funzione rnorm().

[5] Nuovi valori sono generati ogni volta che viene eseguita la funzione rnorm(), quindi se la eseguissimo più volte i risultati sarebbero sempre lievemente diversi, se non predisponessimo il file gauss.csv.

[6] 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. URL consultato il 04/01/2019: https://goo.gl/4Rh53b

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