Curve ROC e valore soglia

Nel post Teorema di Bayes e diagnosi medica abbiamo visto che è possibile impiegare il teorema per rispondere, sulla base dei risultati di un test diagnostico per una specifica malattia, alle domande: 

→ quale probabilità ha di essere malato un soggetto cui il test è risultato positivo? 

→ quale probabilità ha di essere sano un soggetto cui il test è risultato negativo?

Per farlo dobbiamo calcolare il valore soglia che meglio discrimina tra sani e malati, e ricavare la sensibilità e la specificità del test che ad esso corrispondono, impiegando:

→ i risultati del test in questione in un campione rappresentativo di sani e di malati;

→ un metodo per calcolare il valore soglia del test che meglio discrimina tra sani e malati (sensibilità e specificità sono una conseguenza di questo).

I risultati del test sono raccolti eseguendolo, in condizioni controllate, nel corso di una sperimentazione, a un congruo numero di soggetti sani e di soggetti affetti dalla malattia diagnosticata dal test. Il metodo per calcolare il valore soglia applica la strategia diagnostica riportata nella Fig. 1d

che consente di identificare malati (casi) e sani (controlli) massimizzando veri positivi (VP) e veri negativi (VN), minimizzando falsi positivi (FP) e falsi negativi (FN) e individuando quindi il miglior equilibrio tra sensibilità e specificità del test. Alcuni cenni su questa e le altre strategie applicabili sono riportati nel post Teorema di Bayes e diagnosi medica.

La relazione che intercorre tra i risultati del test ottenuti in un campione di soggetti sani, quelli ottenuti in un campione di soggetti malati e la loro curva ROC [1] è riportata in questa figura per un test diagnostico ideale, nel quale le distribuzioni dei risultati del test sono completamente separate, e l'area sotto la curva AUC (Area Under the Curve) è uguale a 1

in quest'altra figura per un test diagnostico inutile nel quale le distribuzioni dei risultati del test sono completamente sovrapposte, non è possibile discriminare, impiegando i risultati dello specifico test diagnostico, tra malati e sani, e l'AUC è uguale a 0.5

e infine in questa figura per un test diagnostico reale, nel quale è presente una parziale sovrapposizione tra i risultati del test nei sani e nei malati e l'AUC è tanto più prossima a 1 quanto meglio il test è in grado di discriminare i soggetti sani da quelli malati.

Quindi la curva ROC fornisce una misura della informazione fornita dal test, espressa come probabilità (misurata dall'AUC) che va da 0.5 (equivalente a decidere se un soggetto è sano o malato mediante il lancio di una moneta) a 1 (il test ci dice con certezza se il paziente è sano o malato). Per impiegare le funzioni del pacchetto pROC, i risultati del test raccolti durante la sperimentazione devono essere strutturati in due variabili, 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) e una variabile quantitativa, qui denominata valore, che per ciascun caso riporta il risultato del test, come appare in questi dati, i primi cinque e gli ultimi cinque del file bayes.csv che ci apprestiamo a impiegare:

classe;valore
0;19
0;22
0;26
0;28
0;29
.....
1;235
1;237
1;237
1;242
1;242

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

Per proseguire ora è necessario:

→ effettuare il download del file bayes.csv 

→ salvare il file nella cartella C:\Rdati\

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

Per effettuare i calcoli impieghiamo alcune funzioni contenute nei pacchetti aggiuntivi pROC e sm [2], che devono essere preventivamente scaricati dal CRAN. Copiate questo script, incollatelo nella Console di R e premete ↵ Invio:

# CURVE ROC E VALORE SOGLIA

#

mydata <- read.table("c:/Rdati/bayes.csv", header=TRUE, sep=";") # importa i dati

attach(mydata) # consente di impiegare direttamente i nomi delle variabili

str(mydata) # mostra la struttura dei dati della tabella importata

head(mydata, n=5) # lista dei primi 5 casi

tail(mydata, n=5) # lista degli ultimi 5 casi

#

sani <- subset(valore, classe==0) # sottoinsieme con i valori dei soggetti sani

malati <- subset(valore, classe==1) # sottoinsieme con i valori dei soggetti malati

#

windows() # apre una nuova finestra

par(mfrow=c(2,2)) # predispone la suddivisione della finestra in quattro quadranti, uno per grafico

#

# primo istogramma

#

hist(malati, xlim=c(0,250), ylim=c(0,0.04), xaxp = c(0, 250, 5), yaxp = c(0, 0.04, 4), breaks=seq(from=0, to=250, by=5), border="red", col="white", main="Istogramma", xlab="Valore osservato", ylab="Densità", freq=FALSE, plot=TRUE) # traccia istogramma malati

hist(sani, xlim=c(0,250), breaks=seq(from=0, to=250, by=5), border="green4", col="lightgreen", freq=FALSE, add=TRUE, plot=TRUE) # sovrappone istogramma sani

legend(150, 0.04, legend=c("Sani", "Malati"), col=c("green", "red"), lty=c(1,1), lwd=c(1,1), cex=0.8) # aggiunge al grafico la legenda

#

# secondo istogramma

#

hist(sani, xlim=c(0,250), ylim=c(0,0.04), xaxp = c(0, 250, 5), yaxp = c(0, 0.04, 4), breaks=seq(from=0, to=250, by=5), border="green4", col="lightgreen", main = "Istogramma", xlab = "Valore osservato", ylab = "Densità", freq=FALSE, plot=TRUE) # traccia istogramma sani

hist(malati, breaks=seq(from=0, to=250, by=5), border="red", col="white", main="", freq=FALSE, add=TRUE, plot=TRUE) # sovrappone istogramma malati

legend(150, 0.04, legend=c("Sani", "Malati"), col=c("green", "red"), lty=c(1,1), lwd=c(1,1), cex=0.8) # aggiunge al grafico la legenda

#

# kernel density plot

#

library(sm) # carica il pacchetto

sm.density.compare(valore, classe, xlim=c(0,250), col=c("green", "red"), main="Kernel density plot", xlab="Valore osservato", ylab="Stima kernel di densità") # traccia kernel density plot separati per sani e malati

title(main="Kernel density plot") # aggiunge al grafico un titolo

legend(150, 0.03, legend=c("Sani", "Malati"), col=c("green", "red"), lty=c(1,2), lwd=c(1,1), cex=0.8) # aggiunge al grafico la legenda

#

# curva ROC

#

library(pROC) # carica il pacchetto

roc(response=classe, predictor=valore, smooth=FALSE, auc=TRUE, ci=TRUE, plot=TRUE, col="orange", identity=TRUE, identity.col="grey", identity.lty=3, main="Curva ROC", xlab="Specificità", ylab="Sensibilità") # traccia la curva ROC e riporta l'AUC con gli intervalli di confidenza

#

# calcola valore soglia, sensibilità e specificità con gli intervalli di confidenza

#

ci.thresholds(response=classe, predictor=valore, thresholds="best", conf.level=0.95, boot.n=100) # calcola il miglior valore soglia tra sani e malati e la sensibilità e la specificità corrispondenti con gli intervalli di confidenza al 95%

ci.se(response=classe, predictor=valore, specificities=seq(0,1,.1), conf.level=0.95, boot.n=100) # calcola gli intervalli di confidenza al 95% della sensibilità per valori di specificità da 0 a 1 con passo .1

ci.sp(response=classe, predictor=valore, sensitivities=seq(0,1,.1), conf.level=0.95, boot.n=100) # calcola gli intervalli di confidenza al 95% della specificità per valori di sensibilità da 0 a 1 con passo .1

#

La parte dello script che precede la creazione dei grafici si spiega da sola, con i commenti inseriti in ciascuna riga. L'unica cosa un po' particolare è la funzione subset() [3] che viene impiegata per generare, dai dati importati, due sottoinsiemi separati, uno contenente solamente i valori dei soggetti sani/controlli e l'altro contenente solamente i valori dei soggetti malati/casi, impiegando come argomenti:

→ il nome della variabile (valore) dalla quale si desidera estrarre il sottoinsieme di valori;

→ il criterio di selezione dei casi da estrarre, laddove classe==0 indica di estrarre i valori per i quali nella variabile classe è riportato 0 (sano) e laddove classe==1 indica di estrarre i valori per i quali nella variabile classe è riportato 1 (malato).

Nel blocco di codice che realizza gli istogrammi sono degni di nota i seguenti argomenti:

→ xlim = c(0,250) indica limite inferiore e limite superiore della scala applicata all'asse delle ascisse;

→ ylim = c(0,0.04) indica limite inferiore e limite superiore della scala applicata all'asse delle ordinate, sulla quale è riportata la densità (di probabilità) della distribuzione dei dati;

→ xaxp = c(0, 250, 5) specifica che sull'asse delle ascisse la scala che va da 0 a 250 deve essere suddivisa in cinque intervalli, quindi sull'asse delle ascisse compariranno sei tacche con i valori 0, 50, 100, 150, 200, 250;

→ yaxp = c(0,0.04,4) specifica che sull'asse delle ordinate la scala che va da 0 e 0.04 deve essere suddivisa in quattro intervalli, quindi sull'asse delle ordinate compariranno cinque tacche con i valori 0, 0.01, 0.02, 0.03, 0.04;

→ breaks=seq(from=0, to=250, by=5) traccia l'istogramma tra 0 e 250 con classi di ampiezza 5;

→ freq=FALSE comporta di riportare sull'asse delle ordinate non il numero di casi osservati in ciascuna classe bensì la densità di probabilità;

→ add=TRUE permette che il secondo istogramma sia sovrapposto al primo.

Per la funzione legend() si rimanda a [4].

Per i kernel density plot sovrapposti l'unica cosa significativa è che sono tracciati mediante la funzione sm.density.compare() del pacchetto aggiuntivo sm, mentre la funzione legend() semplicemente adatta ai nuovi grafici gli argomenti già visti nel caso degli istogrammi.

La curva ROC viene tracciata mediante la funzione roc() che prevede come argomenti:

→ la variabile che contiene la classificazione di ciascun caso (response=caso);

→ la variabile che contiene il valore di ciascun caso (predictor=valore);

→ smooth=FALSE in quanto non vogliamo che la curva ROC sia smussata (ma potete provare a farlo mettendo TRUE);

→ auc=TRUE in quanto vogliamo che sia calcolata l'area sotto la curva AUC;

→ ci=TRUE in quanto vogliamo che sia riportato l'intervallo di confidenza al 95% della AUC;

→ il colore della curva ROC (col="orange");

→ identity=TRUE per riportare nel grafico la retta che va dall'angolo inferiore sinistro all'angolo superiore destro, identity.col per specificarne il colore, identity.lty per specificare lo stile della linea, che rappresenta la curva ROC di un test nel quale l'AUC è uguale a 0.5, i valori dei sani e quelli dei malati sono completamente sovrapposti, e non è possibile discriminare gli uni dagli altri (Fig. 2). In altre parole, se la curva ROC del test coincide con la diagonale, il test fornisce la stessa informazione che viene fornita dal lancio di una moneta;

→ i noti argomenti impiegati per aggiungere a un grafico il titolo e le legende degli assi.

La funzione roc() completa i grafici [5] con la curva ROC del test [6]

quindi riporta il valore dell'AUC (0.9734) con il suo intervallo di confidenza al 95%:

Data: valore in 853 controls (classe 0) < 843 cases (classe 1).

Area under the curve: 0.9734

95% CI: 0.9664-0.9805 (DeLong)

La funzione ci.thresholds() riporta il miglior valore soglia tra sani e malati (74.5), insieme ai valori di specificità (0.9484) e di sensibilità (0.9063) corrispondenti, con i relativi intervalli di confidenza al 95% (conf.level=0.95):

> ci.thresholds(response=classe, predictor=valore, thresholds="best", conf.level=0.95, boot.n=100) # calcola il miglior valore soglia tra sani e malati e la sensibilità e la specificità corrispondenti con gli intervalli di confidenza al 95%

95% CI (100 stratified bootstrap replicates):

 thresholds sp.low sp.median sp.high se.low se.median se.high

       74.5 0.9349    0.9484  0.9643 0.8843    0.9063  0.9242

La funzione ci.se() consente di tabulare i valori di sensibilità corrispondenti ai valori di specificità prescelti, che qui con l'argomento specificities=seq(0,1,.1) sono fatti variare da 0 a 1 con passo 0.1, con i relativi intervalli di confidenza al 95% (conf.level=0.95):

> ci.se(response=classe, predictor=valore, specificities=seq(0,1,.1), conf.level=0.95, boot.n=100) # calcola gli intervalli di confidenza al 95% della sensibilità per valori di specificità da 0 a 1 con passo .1

  sp se.low se.median se.high

 0.0 1.0000    1.0000  1.0000

 0.1 1.0000    1.0000  1.0000

 0.2 0.9941    0.9976  1.0000

 0.3 0.9882    0.9941  0.9988

 0.4 0.9844    0.9912  0.9975

 0.5 0.9754    0.9858  0.9923

 0.6 0.9709    0.9821  0.9899

 0.7 0.9567    0.9678  0.9762

 0.8 0.9480    0.9597  0.9719

 0.9 0.9138    0.9311  0.9453

 1.0 0.6494    0.6821  0.7788

La funzione ci.sp() consente di tabulare i valori di specificità corrispondenti ai valori di sensibilità prescelti, che qui con l'argomento sensitivities=seq(0,1,.1) sono fatti variare da 0 a 1 con passo 0.1, con i relativi intervalli di confidenza al 95% (conf.level=0.95):

> ci.sp(response=classe, predictor=valore, sensitivities=seq(0,1,.1), conf.level=0.95, boot.n=100) # calcola gli intervalli di confidenza al 95% della specificità per valori di sensibilità da 0 a 1 con passo .1

  se sp.low sp.median sp.high

 0.0 1.0000    1.0000  1.0000

 0.1 1.0000    1.0000  1.0000

 0.2 1.0000    1.0000  1.0000

 0.3 1.0000    1.0000  1.0000

 0.4 1.0000    1.0000  1.0000

 0.5 1.0000    1.0000  1.0000

 0.6 1.0000    1.0000  1.0000

 0.7 0.9965    0.9990  1.0000

 0.8 0.9900    0.9953  0.9999

 0.9 0.9312    0.9538  0.9753

 1.0 0.1354    0.1583  0.2446

Le funzioni ci.treshold(), ci.se(), ci.sp() impiegano un metodo di bootstrap per calcolare gli intervalli di confidenza: questo significa che i risultati generati ogni volta saranno lievemente differenti da quelli qui riportati.

L'AUC è uguale a 0.9734 (0.9664-0.9805); il miglior valore soglia tra sani e malati è 74.5 - per un valore inferiore o uguale al valore soglia il test è considerato negativo e per un valore superiore al valore soglia il test è considerato positivo [7].

Per il valore soglia 74.5 abbiamo:

→ un valore di specificità del test (mediana) uguale a 0.9484 che significa che il test è negativo nel 94.8% dei soggetti sani;

→ un valore di sensibilità del test (mediana) uguale a 0.9063 che significa che il test è positivo nel 90.6% dei soggetti malati.

Da notare che una volta fissato il valore soglia tra sani e malati, la sensibilità e la specificità ad esso corrispondenti diventano le due caratteristiche diagnostiche tipiche del test.

Infine per importare ed elaborare i vostri dati, dopo averli strutturati come nel file bayes.csv, dovete semplicemente sostituire nella prima riga di codice c:/Rdati/bayes.csv con nome  e posizione del vostro file di dati. 

Potete trovare:

→ nel post Curve ROC sovrapposte come riportare sullo stesso grafico due curve ROC per un più immediato confronto tra i risultati di due test per la diagnosi della stessa malattia;

→ nel post Sensibilità, specificità e valore predittivo come il valore predittivo del test positivo (la probabilità di essere malato per un soggetto che presenta un test positivo) e il valore predittivo del test negativo (la probabilità di essere sano per un soggetto che presenta un test negativo) possono essere calcolati impiegando il numero di casi osservati VP, FP, FN, VN;

→ nel post Teorema di Bayes e diagnosi medica come il valore predittivo del test positivo e quello del test negativo possono essere calcolati impiegando la sensibilità del test, la specificità del test e la prevalenza della malattia mediante la (iv) e la (v) lì riportate.

Rispondendo così finalmente alle due domande iniziali.

----------

[1] Per una sintesi sul tema delle curve ROC vedere Altman nella pagina Bibliografia.

[2] Trovate la documentazione completa dei pacchetti, indispensabile per impiegarne queste e le altre funzioni, nei loro manuali di riferimento che trovate alla pagina Available CRAN Packages By Name.

https://cran.r-project.org/web/packages/available_packages_by_name.html

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

[4] Vedere il post Aggiungere la legenda a un grafico.

[5] A scopo didattico i grafici sono stati compattati in un'unica figura, ma è possibile riportarli separatamente, impiegando la trasparenza dei colori per un miglior aspetto della sovrapposizione (vedere il post Istogrammi e il post Kernel density plot).

[6] Sull'asse delle ascisse del grafico di una curva ROC sono normalmente riportati i valori 1 – Specificità, valori che vanno da 0 a 1, con lo zero a sinistra e i valori in ordine crescente fino a 1 verso destra: come si confà a un grafico cartesiano. La rappresentazione realizzata con il pacchetto pROC riporta invece sull'asse delle ascisse i valori di Specificità, che, da sinistra verso destra, vanno da 1 a 0: una rappresentazione non canonica, che ha chiaramente lo scopo di evidenziare il rapporto inverso esistente tra sensibilità e specificità.

[7] L'esempio qui impiegato rappresenta la situazione più frequente, quella nella quale la malattia determina un aumento del risultato dei test, ma ovviamente è sufficiente ribaltare dati e conclusioni per riadattare il tutto a un caso nel quale la malattia determina una diminuzione del risultato del test.

Identificare i dati in un grafico xy con un click

Accade con un certa frequenza, nel corso dell'analisi esplorativa dei dati, osservando un grafico di dispersione (xy), di porsi la domanda: questo valore, che si discosta così tanto dagli altri, a quale dato corrisponde?

R fornisce con uno strumento molto efficace per identificare il numero del record/riga che corrisponde ad uno specifico punto riportato in un grafico, la funzione identify() [1], il cui utilizzo viene illustrato con questo script. Copiatelo e incollatelo nella Console di R e premete ↵ Invio:

# IDENTIFICA I DATI IN UN GRAFICO CON UN SEMPLICE CLICK

#

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

attach(ais) # consente di impiegare direttamente i nomi delle variabili

#

windows() # apre una nuova finestra

#

plot(rcc, hc, main="Identifica i punti in un grafico", xlab="Eritrociti (10^12/L)", ylab="Ematòcrito (%)", pch=1) # traccia il grafico di dispersione (xy)

#

# fare click accanto al punto da identificare, ripetere secondo necessità, premere <esc> per terminare

identify(rcc, hc, plot=TRUE, atpen=FALSE, offset=0.5, tolerance=0.25)

#

I dati impiegati sono contenuti nel set di dati ais. Dopo che con la funzione plot() [2] è stato tracciato il grafico di dispersione (xy), la funzione identify() consente di identificare nel grafico uno o più punti impiegando i seguenti argomenti:

→ rcc la variabile in ascisse;

→ hc la variabile in ordinate;

→ plot=TRUE che visualizza l'identificativo del dato sul quale si è fatto click con il mouse;

→ atpen=FALSE che posiziona l'identificativo del dato lievemente discosto dal punto al quale corrisponde (provare a cambiare l'argomento ponendo atpen=TRUE per vedere la differenza);

→ offset=0.5 che indica la distanza che l'identificativo del dato deve avere dal punto al quale corrisponde;

→ tolerance=0.25 che indica la tolleranza ammessa per la distanza definita con l'argomento offset.

Facendo click con il tasto sinistro del mouse accanto a un punto, vedrete comparire il suo identificativo. L'operazione può essere ripetuta per più punti. Di default viene emesso un suono ogniqualvolta si fa click e viene identificato un punto: per silenziare la funzione identify() è necessario aggiungere l'argomento locatorBell=FALSE.

Per uscire dal processo di identificazione e terminare lo script è sufficiente premere il tasto <esc>.

Questo è il grafico ottenuto dopo avere fatto click con il tasto sinistro del mouse sui valori che più si discostano dai rimanenti, valori identificati dalla funzione identify() come corrispondenti ai dati numero 68, 78, 161 e 166.

Come è facile intuire questa funzione grafica avanzata è molto interessante perché permette di individuare e identificare quelli che i test statistici possono avere indicato come possibili dati anomali (o dati aberranti o outliers).

Potete riutilizzare facilmente lo script sostituendo all'oggetto ais l'oggetto contenente i vostri dati, opportunamente strutturati, e modificando conseguentemente i nomi della variabili. Per una guida rapida all'importazione dei dati potete consultare i link:

→ importazione dei dati da un file .csv

→ importazione dei dati da un file .xls o .xlsx

→ gestionedei dati mancanti

----------

[1] Per la documentazione della funzione identify() digitare help(identify) nella Console di R.

[2] Per la documentazione della funzione plot() digitare help(plot) nella Console di R.

Grafici di dispersione (scatterplot)

I grafici di dispersione (scatter plot o scatterplot) sono lo strumento base per visualizzare su un sistema di coordinate cartesiane le coppie di valori x,y di due variabili numeriche nel caso di distribuzioni bivariate.

Prendiamo come esempio i dati ematologici relativi agli eritrociti (o globuli rossi), la cui funzione essenziale come noto è trasportare l'ossigeno dai polmoni a tutti gli altri organi e tessuti, e per i quali valgono in media i seguenti dati (arrotondati per semplicità) :

→ gli eritrociti hanno una concentrazione nel sangue all'incirca di 5·10¹²/L (ovvero ve ne sono all'incirca cinquemila miliardi in un litro di sangue);

→ un eritrocito ha un volume all'incirca di 90 fL (un femtolitro corrisponde a 10^-15 litri quindi un globulo rosso ha un volume all'incirca di 0.000 000 000 000 090 litri).

Da questi dati si ricava che la percentuale del volume del sangue occupata dagli eritrociti (detta valore ematòcrito, e uguale alla concentrazione degli eritrociti moltiplicata per il volume di un eritrocito) è il 45% (la restante parte, esclusi i globuli bianchi e le piastrine, che occupano un volume irrilevante, è liquida ed è rappresentata dal plasma). Ma non solo. Da questi dati si ricava anche che il valore ematòcrito dipende in modo lineare dalla concentrazione degli eritrociti. Possiamo quindi dare un senso al prossimo script, che traccia i grafici di dispersione di queste due grandezze.

I dati sono contenuti nella tabella ais del pacchetto DAAG - accertatevi di avere installato il pacchetto o in alternativa procedete come indicato in [1].

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

# GRAFICI DI DISPERSIONE (scatterplot)

#

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

str(ais) # struttura di ais

attach(ais) # consente di impiegare direttamente i nomi delle variabili

#

windows() # apre e inizializza una nuova finestra grafica

par(mfrow=c(2,2)) # predispone la suddivisione della finestra in quattro quadranti, uno per grafico

#

plot(rcc, hc, xlab="Eritrociti (10^12/L)", ylab="Ematòcrito (%)", main="Grafico di dispersione") # grafico di dispersione semplice

#

plot(rcc, hc, pch=2, col="goldenrod", xlab="Eritrociti (10^12/L)", ylab="Ematòcrito (%)", main="Cambia simbolo e colore") # cambia simbolo e colore

#

plot(rcc, hc, pch=1, col="black", abline(lm(hc~rcc), col="red", lty=1, lwd=1), xlab="Eritrociti (10^12/L)", ylab="Ematòcrito (%)", main="Retta di regressione") # grafico con retta di regressione

#

plot(rcc, hc, xlim=c(3,7), ylim=c(30,70), xaxp=c(3, 7, 4), yaxp=c(30, 70, 4), pch=20, col="black", abline(lm(hc~rcc), col="black", lty=1, lwd=1), xlab="Eritrociti (10^12/L)", ylab="Ematòcrito (%)", main="Ridimensiona gli assi") # ridimensiona gli assi

#

detach(ais) # termina l'impiego diretto dei nomi delle variabili 

#

Dopo avere caricato il pacchetto che include il set di dati ais e dopo averne mostrato la sua struttura con str(ais), viene impiegata la funzione attach(ais) al fine di consente, nelle funzioni successive, di impiegare direttamente i nomi delle variabili incluse nel set di dati.

Per semplificare la presentazione grafica, con par(mfrow=c(2,2)) la finestra viene suddivisa in quattro quadranti, uno per ciascuno dei grafici che verranno preparati.

Il primo grafico, in alto a sinistra, lascia per la maggior parte degli argomenti della funzione plot() i valori di default, specificando solamente quelli essenziali [2]:

→ rcc la variabile da riportare sull'asse delle ascisse x;

→ hc la variabile da riportare sull'asse delle ordinate y;

→ xlab="" l'etichetta da riportare per l'asse delle x;

→ ylab="" l'etichetta riportata per l'asse delle y;

→ main="" il titolo del grafico che compare in alto.

Nel secondo grafico, in alto a destra, sono aggiunti gli argomenti:

→ pch=2 che specifica per i dati l'impiego di un simbolo (triangolo) differente da quello di default (che è il cerchio) [3];

→ col="goldenrod" che specifica per il simbolo dei dati l'impiego di un colore diverso dal nero [4].

Nel terzo grafico, in basso a sinistra, viene riportata una retta mediante la funzione abline() che impiega come argomenti:

→ l'equazione della retta, che altro non è se non la retta di regressione calcolata con la funzione lm() che a sua volta impiega come argomenti hc~rcc per specificare rispettivamente la variabile in ordinate e la variabile in ascisse;

→ il colore della retta (col="red");

→ il tipo di linea da impiegare (lty=1), uno dei sei possibili con i valori da 1 a 6 [5];

→ lo spessore della linea (lwd=1) espresso come multiplo dello spessore di default (1 uguale spessore di default, 2 spessore doppio, eccetera).

Nel quarto grafico, in basso a destra, compaiono quattro nuovi argomenti, che consentono di gestire in modo personalizzabile a piacere le scale degli assi e la loro suddivisione:

→ xlim = c(3,7) indica limite inferiore e limite superiore della scala applicata all'asse delle ascisse, trattandosi di valori di concentrazione degli eritrociti espressi in 10¹²/L questo significa che l'asse delle ascisse riporterà i valori compresi tra 3·10¹²/L e 7·10¹²/L;

→ ylim = c(30,70) indica limite inferiore e limite superiore della scala applicata all'asse delle ordinate, che, trattandosi di valori di percentuale, varieranno quindi tra 30% e 70%;

→ xaxp = c(3, 7, 4) dice che sull'asse delle ascisse la scala che va da 3·10¹²/L a 7·10¹²/L deve essere suddivisa in quattro intervalli, quindi sull'asse delle ascisse compariranno cinque tacche con i valori 3, 4, 50, 6, 7;

→ yaxp = c(30,70,4) dice che sull'asse delle ordinate la scala che va da 30% e 70% deve essere suddivisa in quattro intervalli, quindi sull'asse delle ordinate compariranno cinque tacche con i valori 30, 40, 50, 60, 70.

Lo script si chiude con la funzione detach() che pone termine all'impiego diretto dei nomi delle variabili. In altre parole la variabile rcc ora diventa ais$rcc, la variabile hc ora diventa ais$hc e così via.

I grafici di dispersione riportati qui sopra sono interessanti in quanto sono stati realizzati impiegando esclusivamente le funzioni base di R, tuttavia si può fare di meglio come nell'esempio che segue, che per tracciare lo scatterplot impiega le funzioni del pacchetto ggplot2.

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

# GRAFICO DI DISPERSIONE con ggplot

# colori differenziati in base a un fattore

#

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

library(ggplot2) # carica il pacchetto per la grafica

#

windows() # apre e inizializza una nuova finestra grafica

ggplot(ais, aes(x=rcc, y=hc, color=sport)) + geom_point(size=4) + theme_minimal() # traccia il grafico

#

Ed ecco il grafico, nel quale è interessante notare come la funzione ggplot() consente di attribuire ai punti un differente colore (color=) sulla base dei valori assunti dal fattore (sport) – lo sport praticato – evidenziando così nei dati i corrispondenti sottoinsiemi. Altri esempi che illustrano le varie possibilità di realizzare scatterplot con ggplot2 sono riportati in un post al quale rimando [6].

Vediamo ora due grafici di dispersione multipli che aiutano ad illustrare le relazioni tra le variabili riportando tutti i possibili grafici di dispersione singoli che risultano incrociando i dati ematologici contenuti nelle seguenti variabili del set di dati ais: rcc (eritrociti), wcc (leucociti), hc (ematòcrito), hg (emoglobina) e ferr (ferritina). 

I grafici sono realizzati mediante i pacchetti aggiuntivi car e gclus: se non l'avete già fatto, dovete scaricarli e installarli dal CRAN.

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

# GRAFICI DI DISPERSIONE MULTIPLI 

#

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

library(car) # carica il pacchetto per la grafica

windows() # apre e inizializza una nuova finestra grafica

#

# grafici di dispersione con istogrammi delle variabili

#

scatterplotMatrix(~rcc+wcc+hc+hg+ferr, col="black", pch=20, regLine = list(method=lm, lty=1, lwd=2, col="chartreuse3"), smooth=FALSE, diagonal=list(method ="histogram", breaks="FD"), main="Matrice di dispersione con rette di regressione", data=ais) # mostra il grafico

#

L'argomento ~rcc+wcc+hc+hg+ferr definisce le variabili da rappresentare, e può essere facilmente riscritto per cambiare la rappresentazione a piacere. Ad esempio ponendo questo argomento uguale a ~rcc+wcc+hc+hg+ferr+bmi+ssf+pcBfat+lbm+ht+wt si possono rappresentare in uno stesso grafico tutte le variabili del set di dati ais.

Per la retta di regressione sono impiegati una linea continua (lty=1), uno spessore doppio (lwd=2) e un colore (col="chartreuse3") diverso dal nero, mentre l'argomento smooth=FALSE fa si che non venga aggiunta ai grafici la rappresentazione della funzione non lineare prevista di default nel pacchetto car.

Infine l'argomento diagonal=list() indica la rappresentazione della variabile riportata nella diagonale, che qui è l'istogramma (method ="histogram") con il numero di classi (breaks="FD") calcolato mediante la regola di Freedman-Diaconis. L'istogramma può essere sostituito con altri tipi di rappresentazione, l'elenco delle espressioni che è possibile impiegare in alternativa per questo argomento (provatele) è:

→ diagonal=list(method="density",  bw="nrd0", adjust=1, kernel="gaussian", na.rm=TRUE)

→ diagonal=list(method="boxplot")

→ diagonal=list(method="qqplot")

→ diagonal=list(method="oned")

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

# GRAFICI DI DISPERSIONE MULTIPLI ordinati in base al valore di r

#

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

library(gclus) # carica il pacchetto per la grafica

windows() # apre e inizializza una nuova finestra grafica

mydata <- ais[c(1,2,3,4,5)] # tabella con i dati delle variabili delle colonne 1,2,3,4,5

mydata.r <- abs(cor(mydata)) # calcola la correlazione

mydata.col <- dmat.color(mydata.r) # applica i colori

mydata.o <- order.single(mydata.r) # ordina le variabili, le meglio correlate più vicine alla diagonale

cpairs(mydata, mydata.o, panel.colors=mydata.col, gap=.5, main="Variabili meglio correlate vicine alla diagonale") # mostra il grafico

#

Viene realizzato un grafico di dispersione multiplo che incrocia ancora i dati ematologici degli atleti contenuti nelle cinque variabili rcc (eritrociti), wcc (leucociti), hc (ematòcrito), hg (emoglobina) e ferr (ferritina) ma con due novità.

La prima è che le variabili ora sono individuate non mediante il loro nome bensì indicando, con l'espressione ais[c(1,2,3,4,5)], il set di dati che le contiene e il numero della colonna corrispondente a ciascuna variabile

La seconda novità è che impiegando la funzione cpairs() [7] i singoli grafici di dispersione vengono ordinati in base al valore del coefficiente di correlazione r di Pearson, in modo che le variabili meglio correlate, quelle per le quali la retta di regressione pare essere un'approssimazione ragionevole, siano le più vicine alla diagonale ed evidenziate in colore, aumentando così la chiarezza della rappresentazione grafica.

Potete riutilizzare facilmente lo script sostituendo all'oggetto ais l'oggetto contenente i vostri dati, opportunamente strutturati, e modificando conseguentemente i nomi della variabili. 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

----------

[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] Per la documentazione sugli argomenti della funzione plot() digitare help(plot) nella Console di R, e in aggiunta consultare la documentazione sugli argomenti della funzione par() con help(par).

[3] Per i diversi simboli che è possibile impiegare per tracciare i punti vedere il post I simboli dei punti di R.

[4] Per i colori che è possibile impiegare vedere il post La tabella dei colori di R.

[5] Per le linee che è possibile impiegare vedere il post Gli stili delle linee di R.

[6] Vedere il post Grafici di dispersione (scatterplot) con ggplot.

[7] Digitate help(cpairs) nella Console di R per la documentazione della funzione cpairs() ovvero consultate il manuale di riferimento del pacchetto gclus su: Available CRAN Packages By Name.

https://cran.r-project.org/web/packages/available_packages_by_name.html