venerdì 14 giugno 2019

Curve ROC sovrapposte

Le curve ROC sono illustrate nel post Curve ROC e valore soglia, al quale si rimanda. 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\

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 (vedere anche il post Importazionedei dati da un file .csv). 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 sono illustrate in dettaglio nel post Curve ROC e valore soglia, nel manuale di riferimento del pacchetto [1] oppure 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:









fare click per ingrandire










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 per la diagnosi della malattia X qui messi a confronto può essere effettuata mediante lo script riportato nel post Curve ROC e valore soglia 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

venerdì 10 maggio 2019

Curve ROC e valore soglia

Nel post Teorema di Bayes e diagnosi medica abbiamo visto che è possibile impiegare il teorema per rispondere alle domande: quale probabilità ha di essere affetto da una data malattia un soggetto cui un test per quella malattia è risultato positivo? quale probabilità ha di essere sano un soggetto cui il test è risultato negativo?

Per rispondere a queste due domande impieghiamo il pacchetto pROC per calcolare il valore soglia che meglio discrimina tra sani e malati, e per ricavare la sensibilitàla specificità del test corrispondenti. Per farlo sono necessari:
→ i risultati del test in questione nei sani e nei 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 nella Fig. 2 per un test ideale, nel quale le distribuzioni sono completamente separate, e l'area sotto la curva AUC (Area Under the Curve) è uguale a 1






Fig. 2 - fare click per ingrandire









nella Fig. 3 per un test inutile nel quale le distribuzioni sono completamente sovrapposte, non è possibile discriminare, impiegando i risultati del test, tra malati e sani, e l'AUC è uguale a 0.5







Fig. 3 - fare click per ingrandire








e nella Fig. 4 per un test 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.







Fig. 4 - fare click per ingrandire








Quindi la curva ROC fornisce una misura della informazione fornita dal test che va da 0.5 (equivalente a decidere se un soggetto è sano o malato mediante il lancio di una moneta) a 1 (il test fornisce una risposta certa). 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\

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
#

Nota: per importare ed elaborare i vostri dati, dopo averli strutturati come nel file bayes.csv dovete solo modificare opportunamente il nome e la posizione del vostro file di dati nella prima riga di codice. 

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 al post Aggiungere la legenda a un grafico.

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 con la curva ROC del test [4]










fare click per ingrandire











quindi riporta il valore dell'AUC 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, insieme ai valori di specificità e di sensibilità 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 al bisogno 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 al bisogno 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 calcolato è 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 [5]. Per questo valore soglia 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, la sensibilità e la specificità corrispondenti non dipendono né dal numero di sani né dal numero di malati né dal loro rapporto numerico. Semplicemente il loro valore sarà tanto meno incerto (con un intervallo di confidenza ristretto) quanto più numericamente consistente sarà la casistica.

Nel post Sensibilità, specificità, valore predittivo vediamo 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 o, in alternativa, possono essere calcolati impiegando la sensibilità del test, la specificità del test e la prevalenza della malattia mediante la (iv) e la (v) riportate nel post Teorema di Bayes e diagnosi medica. Rispondendo così finalmente alle domande riportate all'inizio.


----------

[1] Per una sintesi sul tema delle curve ROC vedere Altman nella pagina della 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. URL consultato il 10/05/2019: https://goo.gl/hLC9BB

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

[4] 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à.

[5] 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.