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’. URL consultato il 14/06/2019: http://bit.ly/2HBk7R6