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