Curve ROC sovrapposte

Le curve ROC sono già state illustrate [1]. Ora aggiungiamo un breve script che consente di riportare sullo stesso grafico le curve ROC di due diversi test per la diagnosi di una malattia X, qui denominati test 1 e test 2, onde confrontare anche visivamente la loro capacità di discriminare tra sani e malati, espressa numericamente come AUC (Area Under the Curve, area sotto la curva ROC).

Per proseguire è necessario:

→ effettuare il download del file test_1.csv e del file test_2.csv

→ salvare entrambi i file nella cartella C:\Rdati\

Trovate i file, con le istruzioni per scaricare questi e gli altri file di dati impiegati nei post, alla pagina Dati.

Volendo impiegare le funzioni del pacchetto pROC i file di dati devono contenere i risultati del/i test strutturati in due variabili, una variabile quantitativa, qui denominata valore, che per ciascun caso riporta il risultato del test, e una variabile qualitativa, in R definita “fattore”, e qui denominata classe, che contiene la classe di appartenenza (qui vengono impiegati i valori previsti di default nel pacchetto e cioè 0 per indicare un soggetto sano/controllo e 1 per indicare un soggetto malato/caso) come vedete ad esempio aprendo con un editor di testo il file test_1.csv che contiene i risultati del test 1:

valore;classe

19;0

22;0

22;1

24;1

24;1

26;0

Quindi la coppia di valori 19;0 indica che in un soggetto sano/controllo (0) il risultato fornito dal test era 19; la coppia di valori 22;0 indica che in un soggetto sano/controllo (0) il risultato fornito dal test era 22; la coppia di valori 22;1 indica che in un soggetto malato/caso (1) il risultato fornito dal test era 22; e così via.

Per effettuare i calcoli sono impiegate le funzioni del pacchetto aggiuntivo pROC [1] che deve essere scaricato, se non ancora fatto in precedenza, dal CRAN. Copiate questo script, incollatelo nella Console di R e premete ↵ Invio:

# CONFRONTO TRA DUE CURVE ROC

#

# importa i dati dei due test da mettere a confronto

test_1 <- read.table("c:/Rdati/test_1.csv", header=TRUE, sep=";")

test_2 <- read.table("c:/Rdati/test_2.csv", header=TRUE, sep=";")

#

library(pROC)

#

# traccia la curva ROC del primo test

roc(response=test_1$classe, predictor=test_1$valore, smooth=FALSE, auc=TRUE, ci=FALSE, plot=TRUE, identity=FALSE, main="Confronto tra due curve ROC", xlab="Specificità", ylab="Sensibilità")

# calcola valore soglia, sensibilità e specificità del primo test

ci.thresholds(response=test_1$classe, predictor=test_1$valore, thresholds="best", conf.level=0.95, boot.n=100)

#

# sovrappone la curva ROC del secondo test

roc(response=test_2$classe, predictor=test_2$valore, smooth=FALSE, auc=TRUE, ci=FALSE, plot=TRUE, add=TRUE, col="red", lty=4)

# calcola valore soglia, sensibilità e specificità del secondo test

ci.thresholds(response=test_2$classe, predictor=test_2$valore, thresholds="best", conf.level=0.95, boot.n=100)

#

# aggiunge una legenda al grafico

legend(0.35, 0.8, legend=c(paste("Test 1, AUC =", round(auc(response=test_1$classe, predictor=test_1$valore), digits=3), sep=" "), paste("Test 2, AUC =", round(auc(response=test_2$classe, predictor=test_2$valore), digits=3), sep=" ")), col=c("black", "red"), lty=c(1,4), cex=0.8) # aggiunge al grafico la legenda

#

Le prime due righe di codice provvedono a importare i dati. Con la terza viene caricato il pacchetto pROC.

Con la quarta e la quinta riga di codice le funzioni roc() e ci.treshold() sono impiegate rispettivamente per tracciare il grafico e calcolare le statistiche del primo test. Entrambe le funzioni le trovate illustrate in dettaglio in [1], nel manuale di riferimento del pacchetto [1] e anche digitando help(roc) e help(ci.treshold) nella Console di R.

Nelle due righe successive la funzione roc() viene impiegata con l'argomento add=TRUE per sovrapporre la seconda curva ROC alla prima, quindi con la funzione ci.treshold() sono calcolate le statistiche del secondo test.

Infine con l'ultima riga di codice al grafico viene aggiunta una legenda impiegando la funzione legend() i cui argomenti sono:

→ 0.3 per indicare il valore in ascisse al quale posizionare l'angolo superiore sinistro del box della legenda;

→ 0.75 per indicare il valore in ordinate al quale posizionare l'angolo superiore sinistro del box della legenda;

→ legend che riporta con la funzione c() le etichette che descrivono il primo e il secondo test;

→ col che riprende con la funzione c() i colori assegnati alle due curve ROC;

→ lty che riprende lo stile delle linee delle due curve ROC;

→ cex che dimensiona i caratteri della legenda.

Da notare che la funzione legend() è stata condensata in una sola riga di codice sfruttando una caratteristica fondamentale della programmazione a oggetti, cioè il fatto che una funzione può impiegare come valore per un argomento il risultato di una funzione, che a sua volta può impiegare come valore per un argomento il risultato di una funzione, e così via, come avviene in questo caso, nel quale:

→ la funzione legend(), che costruisce la legenda,

→ impiega come terzo argomento legend, che riporta il contenuto della legenda

→ che a sua volta impiega il risultato della funzione c(), che fornisce i valori all'argomento legend

→ impiegando come primo argomento la funzione paste()

→ che impiega come secondo argomento il risultato della funzione round(), che arrotonda i risultati

→ impiegando come argomento il risultato della funzione auc(), che calcola il valore dell'area sotto la curva AUC impiegando come argomenti i dati del file test_1.csv

ripetendo infine gli ultimi tre passi anche per i dati del file test_2.csv. Anche se va a scapito della leggibilità del codice, che di norma è preferibile spezzare su più righe introducendo opportune variabili accessorie, questo è comunque un esempio della possibilità di concisione offerta da R.

Ecco il grafico ottenuto mediante lo script:

La conclusione è che per la diagnosi della malattia X è preferibile impiegare il test 1, che con una AUC = 0.963 fornisce un'informazione diagnostica superiore a quella del test 2 (AUC = 0.902). E questo è espresso anche dal fatto che il test 1 ha sensibilità e specificità entrambe superiori a quelle del test 2.

Una analisi separata di ciascuno dei due test diagnostici qui messi a confronto può essere effettuata mediante lo script riportato in [1] semplicemente cambiando nella prima riga di codice il nome del file, sostituendo quindi bayes.csv prima con test_1.csv e poi con test_2.csv.

----------

[1] Vedere il post Curve ROC e valore soglia.

[1] Per la funzione roc() e la funzione ci.tresholds() vedere il manuale di riferimento del pacchetto Package ‘pROC’.

https://cran.r-project.org/web/packages/pROC/pROC.pdf

Coefficienti di correlazione parametrici e non parametrici

Come sottolineato altrove, l'assunzione implicita di un rapporto di causa-effetto in seguito ad un coefficiente di correlazione significativo è una delle trappole mentali più ricorrenti e più subdole della statistica: perché non esiste una connessione logica tra correlazione e causazione [1].

Nonostante da questo derivi l'invito ad impiegare la correlazione con parsimonia, con prudenza e con spirito critico nel trarne le conclusioni, risulta difficile ignorare l'esistenza in R di due strumenti.

Il primo sono i classici coefficienti di correlazione e precisamente:

→ il coefficiente di correlazione (lineare) r di Pearson, il classico test parametrico;

→ il coefficiente di correlazione per ranghi ρ (rho) di Spearman (test non parametrico);

→ il coefficiente di correlazione τ (tau) di Kendall (test non parametrico).

Il secondo sono i correlogrammi [2], una rappresentazione un po' particolare ma interessante della correlazione tra variabili, utile come complemento grafico ai coefficienti di correlazione di cui sopra. 

Va ricordato che il più ampiamente usato (e abbondantemente abusato) coefficiente di correlazione r fornisce una misura del grado di correlazione lineare ed è pertanto valido solamente se le due variabili a confronto hanno una relazione di tipo lineare, mentre i coefficienti di correlazione non parametrici forniscono una misura della correlazione valida anche anche nel caso di relazioni non lineari.

In questo primo script viene impiegato il set di dati galton [3], nel quale sono riportate le altezze (in pollici) di 928 coppie padre-figlio. Sono richiesti i pacchetti aggiuntivi psychTools e Hmisc, che devono essere preventivamente scaricati e installati dal CRAN.

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

# COEFFICIENTI DI CORRELAZIONE parametrici e non parametrici

#

library(psychTools) # carica il pacchetto che include il set di dati galton

str(galton) # mostra la struttura dei dati

#

# calcola i coefficienti di correlazione con i metodi di pearson (il classico r), spearman, kendall

#

round(cor(galton, use="complete.obs", method="pearson"), digits=2) # r di Pearson

round(cor(galton, use="complete.obs", method="spearman"), digits=2) # rho di Spearman

round(cor(galton, use="complete.obs", method="kendall"), digits=2) # tau di Kendall

#

# calcola i coefficienti di correlazione con i livelli di significatività

#

library(Hmisc) # carica il pacchetto

rcorr(as.matrix(galton), type="pearson") # type può essere solo “pearson” (il classico r) o “spearman”

rcorr(as.matrix(galton), type="spearman")

#

Questi sono i risultati ottenuti con i tre metodi per il calcolo della correlazione:

> round(cor(galton, use="complete.obs", method="pearson"), digits=2) # r di Pearson 

       parent child

parent   1.00  0.46

child    0.46  1.00

> round(cor(galton, use="complete.obs", method="spearman"), digits=2) # rho di Spearman 

       parent child

parent   1.00  0.43

child    0.43  1.00

> round(cor(galton, use="complete.obs", method="kendall"), digits=2) # tau di Kendall 

       parent child

parent   1.00  0.33

child    0.33  1.00

La funzione cor() viene impiegata per calcolare i coefficienti di correlazione con gli argomenti:

→ galton che indica i dati da elaborare;

→ use=”complete.obs” che comporta l'eliminazione automatica dei casi nei quali manchi l'uno o l'altro dato, argomento che qui non servirebbe impiegare, ma che viene riportato per completezza;

→ method= che di volta in volta riporta i valori "pearson", "spearman", "kendall", per calcolare le rispettive statistiche.

Da notare che la funzione cor() è stata annidata nella funzione round() che con l'argomento digits=2 provvede ad arrotondare a due cifre decimali i risultati: e che i coefficienti di correlazione non parametrici forniscono, come atteso, valori inferiori a quelli del classico coefficiente di correlazione parametrico r di Pearson (sono più conservativi).

Infine lo script impiega la funzione rcorr() del pacchetto Hmisc per calcolare la significatività del coefficiente r di Pearson e del coefficiente ρ (rho) di Spearman. La probabilità p di osservare per caso i valori di r e di rho qui ottenuti è molto bassa e, come potrete vedere se eseguite lo script, viene riportata come uguale a 0 (zero). Quindi entrambi questi coefficienti di correlazione risultano essere significativi.

In questo secondo script sono impiegati i dati contenuti nella tabella ais del pacchetto DAAG - accertatevi di avere installato il pacchetto o in alternativa procedete come indicato in [4].

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

# COEFFICIENTI DI CORRELAZIONE parametrici e non parametrici

#

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

str(ais) # struttura di ais

mydata <- ais[c(1,2,3,4,5,6,7,8,9,10,11)] # salva sole le colonne con le variabili numeriche in mydata

str(mydata) # mostra la struttura di mydata

#

# method può essere pearson (il classico r), spearman, kendall

#

round(cor(mydata, use="complete.obs", method="pearson"), digits=2) # r di Pearson

round(cor(mydata, use="complete.obs", method="spearman"), digits=2) # rho di Spearman

round(cor(mydata, use="complete.obs", method="kendall"), digits=2) # tau di Kendall

#

# calcola i coefficienti di correlazione con i livelli di significatività

#

library(Hmisc) # carica il pacchetto

rcorr(as.matrix(mydata), type="pearson") # type può essere solo “pearson” (il classico r) o “spearman”

rcorr(as.matrix(mydata), type="spearman")

#

# approfondimenti grafici sulla correlazione lineare

#

library(car) # carica il pacchetto

windows() # apre una finestra grafica

scatterplotMatrix(~rcc+wcc+hc+hg+ferr+bmi+ssf+pcBfat+lbm+ht+wt, regLine = list(method=lm, lty=1, lwd=2, col="red"), smooth=FALSE, diagonal=list(method="density", bw="nrd0", adjust=1, kernel="gaussian", na.rm=TRUE), col = "black", main="Matrice di dispersione", data=ais) # traccia il grafico di dispersione xy che incrocia tutte le variabili

#

Nella parte iniziale dello script dopo avere caricato il pacchetto DAAG con la funzione str(ais) viene mostrata la struttura dell'oggetto/set di dati, che contiene 13 variabili, 11 numeriche e 2 non numeriche. Dato che siamo interessati alle sole variabili numeriche, queste, che sono comprese nelle colonne da 1 a 11 della tabella (o dataframe o dataset) ais, le salviamo in un nuovo oggetto denominato mydata che è quello sul quale lavoreremo. Infine per la conferma dell'avvenuta operazione di selezione delle colonne viene mostrata la struttura del nuovo oggetto con str(mydata).

Il resto del codice è praticamente identico allo script precedente. I risultati sono presentati sotto forma di una matrice che contiene i valori dei coefficienti di correlazione che la funzione cor() calcola per tutte le coppie di variabili. Questi valori sono simmetrici rispetto alla diagonale, nella quale è riportato sempre un valore uguale a 1 per i coefficienti di correlazione di ciascuna variabile con se stessa. I valori evidenziati (manualmente dal sottoscritto) sono quelli che nel successivo grafico mostrano una relazione tra le variabili che più si avvicina ad una retta.

> round(cor(mydata, use="complete.obs", method="pearson"), digits=2) # r di Pearson

         rcc  wcc    hc    hg  ferr  bmi   ssf pcBfat   lbm    ht   wt

rcc     1.00 0.15  0.92  0.89  0.25 0.30 -0.40  -0.49  0.55  0.36 0.40

wcc     0.15 1.00  0.15  0.13  0.13 0.18  0.14   0.11  0.10  0.08 0.16

hc      0.92 0.15  1.00  0.95  0.26 0.32 -0.45  -0.53  0.58  0.37 0.42

hg      0.89 0.13  0.95  1.00  0.31 0.38 -0.44  -0.53  0.61  0.35 0.46

ferr    0.25 0.13  0.26  0.31  1.00 0.30 -0.11  -0.18  0.32  0.12 0.27

bmi     0.30 0.18  0.32  0.38  0.30 1.00  0.32   0.19  0.71  0.34 0.85

ssf    -0.40 0.14 -0.45 -0.44 -0.11 0.32  1.00   0.96 -0.21 -0.07 0.15

pcBfat -0.49 0.11 -0.53 -0.53 -0.18 0.19  0.96   1.00 -0.36 -0.19 0.00

lbm     0.55 0.10  0.58  0.61  0.32 0.71 -0.21  -0.36  1.00  0.80 0.93

ht      0.36 0.08  0.37  0.35  0.12 0.34 -0.07  -0.19  0.80  1.00 0.78

wt      0.40 0.16  0.42  0.46  0.27 0.85  0.15   0.00  0.93  0.78 1.00

> round(cor(mydata, use="complete.obs", method="spearman"), digits=2) # rho di Spearman

         rcc  wcc    hc    hg  ferr  bmi   ssf pcBfat   lbm    ht    wt

rcc     1.00 0.17  0.91  0.89  0.25 0.29 -0.40  -0.50  0.59  0.40  0.42

wcc     0.17 1.00  0.17  0.14  0.05 0.23  0.15   0.14  0.11  0.04  0.17

hc      0.91 0.17  1.00  0.95  0.25 0.32 -0.44  -0.54  0.63  0.43  0.45

hg      0.89 0.14  0.95  1.00  0.31 0.37 -0.43  -0.54  0.67  0.41  0.49

ferr    0.25 0.05  0.25  0.31  1.00 0.30 -0.11  -0.19  0.34  0.17  0.30

bmi     0.29 0.23  0.32  0.37  0.30 1.00  0.29   0.15  0.70  0.35  0.84

ssf    -0.40 0.15 -0.44 -0.43 -0.11 0.29  1.00   0.96 -0.22 -0.10  0.13

pcBfat -0.50 0.14 -0.54 -0.54 -0.19 0.15  0.96   1.00 -0.41 -0.25 -0.05

lbm     0.59 0.11  0.63  0.67  0.34 0.70 -0.22  -0.41  1.00  0.80  0.91

ht      0.40 0.04  0.43  0.41  0.17 0.35 -0.10  -0.25  0.80  1.00  0.78

wt      0.42 0.17  0.45  0.49  0.30 0.84  0.13  -0.05  0.91  0.78  1.00

> round(cor(mydata, use="complete.obs", method="kendall"), digits=2) # tau di Kendall

         rcc  wcc    hc    hg  ferr  bmi   ssf pcBfat   lbm    ht    wt

rcc     1.00 0.11  0.75  0.71  0.17 0.19 -0.26  -0.33  0.41  0.27  0.28

wcc     0.11 1.00  0.11  0.09  0.03 0.15  0.10   0.10  0.08  0.03  0.12

hc      0.75 0.11  1.00  0.83  0.17 0.22 -0.29  -0.35  0.44  0.29  0.30

hg      0.71 0.09  0.83  1.00  0.22 0.25 -0.29  -0.36  0.46  0.28  0.33

ferr    0.17 0.03  0.17  0.22  1.00 0.20 -0.07  -0.13  0.22  0.11  0.20

bmi     0.19 0.15  0.22  0.25  0.20 1.00  0.20   0.11  0.52  0.24  0.65

ssf    -0.26 0.10 -0.29 -0.29 -0.07 0.20  1.00   0.82 -0.13 -0.07  0.09

pcBfat -0.33 0.10 -0.35 -0.36 -0.13 0.11  0.82   1.00 -0.24 -0.16 -0.01

lbm     0.41 0.08  0.44  0.46  0.22 0.52 -0.13  -0.24  1.00  0.61  0.77

ht      0.27 0.03  0.29  0.28  0.11 0.24 -0.07  -0.16  0.61  1.00  0.59

wt      0.28 0.12  0.30  0.33  0.20 0.65  0.09  -0.01  0.77  0.59  1.00

L'ultima parte dello script realizza un grafico che fornisce la necessaria integrazione grafica all'analisi numerica.

Impiegando la funzione scatterplotMatrix() viene realizzato con una sola riga di codice un grafico che riporta sulla diagonale per ciascuna variabile il kernek density plot, e per ciascuna coppia di variabili il grafico di dispersione (xy) con la relativa retta di regressione.

I dettagli delle funzioni qui impiegate, tutte molto semplici, li trovate come al solito digitando help(nomedellafunzione) nella Console di R.

----------

[1] Vedere il post Correlazione e causazione.

[2] Vedere il post Correlogrammi.

[3] Vedere il post Il set di dati Galton.

[4] Vedere il post Il set di dati ais nel quale trovate anche come caricare i dati della tabella senza impiegare il pacchetto DAAG. 

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 rappresentare 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 piccola differenza (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 cioè letteralmente la "Deviazione Mediana Assoluta 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).