Sensibilità, specificità e valore predittivo

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

Nel post Curve ROC e valore soglia abbiamo visto, dati i risultati di un test per la diagnosi di una malattia X, come impiegare il pacchetto pROC per tracciare la curva ROC [1] e individuare il valore soglia che consente di classificare correttamente il maggior numero possibile di soggetti sani (VN) e di classificare correttamente il maggior numero possibile di soggetti malati (VP), minimizzando sia la possibilità di sbagliare classificando dei soggetti malati come sani (FN) sia la possibilità di sbagliare classificando dei soggetti sani come malati (FP). Dato il valore soglia individuato, risultano conseguentemente fissate la sensibilità e la specificità del test.

Ora impieghiamo il pacchetto epiR per calcolare le probabilità che rispondono alle due domande e cioè:

→ il valore predittivo di un test positivo, ovvero la probabilità di essere malato per un soggetto che presenta un test positivo, che risponde alla prima domanda

→ il valore predittivo di un test negativo, ovvero la probabilità di essere sano per un soggetto che presenta un test negativo, che risponde alla seconda domanda.

Si stabilisce se il test è positivo o negativo impiegando il valore soglia individuato, quindi si calcola il valore predittivo del test impiegando:

→ la sensibilità del test, cioè la probabilità che il test sia positivo in un soggetto malato;

→ la specificità del test, cioè la probabilità che il test sia negativo in un soggetto sano;

→ la prevalenza, cioè la probabilità della malattia.

La sensibilità e la specificità sono caratteristiche intrinseche al test, dipendono dal valore soglia individuato e sono indipendenti dalla prevalenza, che è invece una caratteristica della popolazione alla quale il test viene applicato.

Notare che la prevalenza, cioè la probabilità della malattia nella popolazione alla quale il test viene applicato, è la probabilità a priori della malattia, ovvero la probabilità per un dato paziente di essere malato (o sano) prima di avere eseguito il test. Il teorema di Bayes impiega la probabilità a priori (nel nostro caso la prevalenza della malattia), insieme a sensibilità e specificità del test, per calcolare la probabilità a posteriori, ovvero la probabilità per un dato paziente di essere malato (o sano) dopo avere eseguito il test (il valore predittivo del test positivo o del test negativo).

Notare inoltre che nella (iii) del post Teorema di Bayes e diagnosi medica:

→ P(T+|M+)/P(T+) è una costante e rappresenta l'evidenza fornita dai fatti (nel nostro caso il risultato del test per la diagnosi della malattia)

→ P(M+) è l'informazione/probabilità a priori

→ P(M+|T+) è l'informazione/probabilità a posteriori

Pertanto il teorema può essere riscritto come

probabilità a posteriori = evidenza · probabilità a priori

espressione che sintetizza l'essenza del teorema di Bayes: che consente di aggiornare l'informazione in nostro possesso, la probabilità a priori, sulla base dell'evidenza fornita dai fatti.

Mentre i calcoli eseguiti impiegando i valori di probabilità sono riportati nell'appendice al termine del post, vediamo qui i calcoli eseguiti impiegando il numero di casi osservati, organizzati in questo modo

	Malattia +	Malattia -
Test +	VP	FP
Test -	FN	VN

e definiti come [2]:

→ veri positivi (VP), i soggetti malati con il test positivo;

→ falsi positivi (FP), i soggetti sani con il test positivo;

→ falsi negativi (FN), i soggetti malati con il test negativo;

→ veri negativi (VN), i soggetti sani con il test negativo.

Supponiamo ora (primo esempio) che il test per una data malattia sia stato eseguito in 800 soggetti sani, dipendenti di un Ospedale, impiegati come controlli, e in 800 malati, che sono stati raccolti in numero così elevato in quanto nell'Ospedale esiste un Reparto che si occupa specificamente di questa malattia. Questa è la tabella con i risultati del test

	Malattia +	Malattia -
Test +	700 (VP)	44 (FP)
Test -	100 (FN)	756 (VN)

I calcoli sono effettuati mediante il pacchetto aggiuntivo epiR [3] che deve essere preventivamente scaricato dal CRAN. Copiate questo script, incollatelo nella Console di R e premete ↵ Invio:

# SENSIBILITA', SPECIFICITA' E VALORE PREDITTIVO

# 

#

library(epiR) # carica il pacchetto

cells <- c(700,44,100,756) # genera l'array cells con i valori numerici contenuti nelle celle

rnames <- c("Test +", "Test -") # genera l'array rnames con i nomi delle righe

cnames <- c("Malattia +", "Malattia -") # genera l'array cnames con i nomi delle colonne

mymatrix <- matrix(cells, nrow=2, ncol=2, byrow=TRUE, dimnames=list(rnames, cnames)) # costruisce la matrice dati mymatrix a partire dagli array cells, rnames, cnames

mymatrix # mostra la matrice dati mymatrix

mystat <- epi.tests(mymatrix, conf.level=0.95) # calcola le statistiche

summary(mystat) # mostra le statistiche

#

Con la funzione c() [4] sono generati gli array che mediante l'operatore di assegnamento <- vengono denominati cells, rnames, cnames e che ora contengono rispettivamente i dati (cells), i nomi delle righe (rnames) e i nomi delle colonne (cnames).

Quindi mediante la funzione matrix() [5] gli array sono combinati nella matrice mymatrix di 2 righe (nrow=2) per due colonne (ncol=2) inserendo i valori per riga (byrow=TRUE) e quindi da sinistra a destra e dall'alto in basso ottenendo la matrice desiderata alla quale con l'ultimo argomento (dimnames=list(rnames, cnames)) sono assegnati i nomi alle righe e alle colonne.

Questo è il contenuto dell'oggetto mymatrix riportato da R con il comando mymatrix:

> mymatrix # mostra la matrice dati mymatrix

       Malattia + Malattia -

Test +        700         44

Test -        100        756

Impiegando la funzione epi.test() sono infine calcolate sui dati contenuti in mymatrix le statistiche, ciascuna con l'intervallo di confidenza al 95% (conf.level=0.95), e i risultati sono riportati con la funzione summary():

> summary(mystat) # mostra le statistiche

   statistic          est       lower        upper

1         ap   0.46500000  0.44033136   0.48979728

2         tp   0.50000000  0.47520704   0.52479296

3         se   0.87500000  0.85006839   0.89712668

4         sp   0.94500000  0.92686534   0.95975609

5    diag.ac   0.91000000  0.89490175   0.92357301

6    diag.or 120.27272727 83.14513765 173.97925284

7       nndx   1.21951220  1.16702078   1.28711106

8     youden   0.82000000  0.77693373   0.85688277

9     pv.pos   0.94086022  0.92141771   0.95670377

10    pv.neg   0.88317757  0.85974563   0.90393102

11    lr.pos  15.90909091 11.92293523  21.22792488

12    lr.neg   0.13227513  0.11003359   0.15901245

13    p.rout   0.53500000  0.51020272   0.55966864

14     p.rin   0.46500000  0.44033136   0.48979728

15    p.tpdn   0.05500000  0.04024391   0.07313466

16    p.tndp   0.12500000  0.10287332   0.14993161

17    p.dntp   0.05913978  0.04329623   0.07858229

18    p.dptn   0.11682243  0.09606898   0.14025437

Le prime 12 statistiche calcolate sono quelle che ci interessano, riportiamo (tra parentesi) il significato della sigla, le formule con cui sono calcolate e il risultato numerico ottenuto (est) con il limite inferiore (lower) e il limite superiore (upper) del relativo intervallo di confidenza al 95% (le statistiche da 13 a 18 riguardano l'outcome cioè il risultato, e non le consideriamo).

ap (prevalenza apparente, soggetti con il test positivo)

(VP+FP)/(VP+FP+FN+VN) = 0.4650000 (0.4403314 - 0.4897973)

tp (prevalenza reale, soggetti con la malattia)

(VP+FN)/(VP+FP+FN+VN) = 0.5000000 (0.4752070 - 0.5247930)

se (sensibilità, positività nei malati)

VP/(VP+FN) = 0.8750000 (0.8500684 - 0.8971267)

sp (specificità, negatività nei sani)

VN/(VN+FP) = 0.9450000 (0.9268653 - 0.9597561)

diag.ac (accuratezza diagnostica)

(VP+VN)/(VP+FP+FN+VN) = 0.9100000 (0.8949017 - 0.9235730)

diag.or (odd ratio)

rapporto LR+/LR- = 120.2727273 (83.1451376 - 173.9792528)

nndx (numero necessario per la diagnosi)

1/(sensibilità - (1 - specificità)) = 1.2195122 (1.1670208 - 1.2871111)

youden (indice di Youden)

sensibilità + specificità - 1 = 0.8200000 (0.7769337 - 0.8568828)

pv.pos (valore predittivo di un test positivo)

VP/(VP+FP) = 0.9408602 (0.9214177 - 0.9567038)

pv.neg (valore predittivo di un test negativo)

VN/(VN+FN) = 0.8831776 (0.8597456 - 0.9039310)

lr.pos (rapporto di verosimiglianza LR+ per un test positivo) [6]

(VP/(VP+FN))/(FP/(FP+VN)) = 15.9090909 (11.9229352 - 21.2279249)

lr.neg (rapporto di verosimiglianza LR- per un test negativo) [6]

(FN/(VP+FN))/(VN/(FP+VN)) = 0.1322751 (0.1100336 - 0.1590125)

La conclusione è che il valore predittivo di un test positivo, ovvero la probabilità di essere malato per un soggetto al quale il test è risultato positivo, è uguale al 94.1% e il valore predittivo di un test negativo, ovvero la probabilità di essere sano per un soggetto al quale il test è risultato negativo, è uguale all'88.3%.

Nel caso appena esaminato la prevalenza della malattia è del 50%, un valore che è dipeso dalle modalità di raccolta dei dati: per valutare adeguatamente il test è stato raccolto il maggior numero possibile di malati, un fatto che non riflette la reale diffusione della malattia. Ora supponiamo che il test venga applicato per valutare lo stato di salute di un individuo della popolazione generale, nella quale la reale diffusione della malattia, la sua prevalenza, è del 2%. I dati (secondo esempio) sono questi:

	Malattia +	Malattia -
Test +	14 (VP)	44 (FP)
Test -	2 (FN)	756 (VN)

Copiate questo script, identico al precedente tranne che per il fatto che contiene i nuovi dati, incollatelo nella Console di R e premete ↵ Invio:

# SENSIBILITA', SPECIFICITA' E VALORE PREDITTIVO

# 

#

library(epiR) # carica il pacchetto

cells <- c(14,44,2,756) # genera l'array cells con i valori numerici contenuti nelle celle

rnames <- c("Test +", "Test -") # genera l'array rnames con i nomi delle righe

cnames <- c("Malattia +", "Malattia -") # genera l'array cnames con i nomi delle colonne

mymatrix <- matrix(cells, nrow=2, ncol=2, byrow=TRUE, dimnames=list(rnames, cnames)) # costruisce la matrice dati mymatrix a partire dagli array cells, rnames, cnames

mymatrix # mostra la matrice dati mymatrix

mystat <- epi.tests(mymatrix, conf.level=0.95) # calcola le statistiche

summary(mystat) # mostra le statistiche

#

Questi sono i dati contenuti nell'oggetto mymatrix

> mymatrix # mostra la matrice dati mymatrix

            Malattia + Malattia -

Test +         14         44

Test -          2        756

e questi sono, riportati con la funzione summary(), i risultati calcolati con la funzione epi.test(), ciascuno con l'intervallo di confidenza al 95% (conf.level=0.95)

> summary(mystat) # mostra le statistiche

   statistic           est         lower         upper

1         ap   0.071078431  0.0544119039   0.090919329

2         tp   0.019607843  0.0112480876   0.031647030

3         se   0.875000000  0.6165237632   0.984486396

4         sp   0.945000000  0.9268653430   0.959756088

5    diag.ac   0.943627451  0.9255197857   0.958436211

6    diag.or 120.272727273 26.5044169172 545.778047894

7       nndx   1.219512195  1.0590499968   1.840301892

8     youden   0.820000000  0.5433891062   0.944242484

9     pv.pos   0.241379310  0.1387012000   0.371689716

10    pv.neg   0.997361478  0.9905015414   0.999680303

11    lr.pos  15.909090909 11.3036587611  22.390907130

12    lr.neg   0.132275132  0.0361754961   0.483661941

13    p.rout   0.928921569  0.9090806715   0.945588096

14     p.rin   0.071078431  0.0544119039   0.090919329

15    p.tpdn   0.055000000  0.0402439118   0.073134657

16    p.tndp   0.125000000  0.0155136038   0.383476237

17    p.dntp   0.758620690  0.6283102838   0.861298800

18    p.dptn   0.002638522  0.0003196972   0.009498459

Si nota immediatamente che quando la prevalenza della malattia diminuisce dal 50% al 2% (precisamente diventa 1.96%)

                 est       lower        upper

tp        0.5000000   0.4752070    0.5247930 [primo esempio]

tp       0.01960784  0.01124809   0.03164703 [secondo esempio]

il valore predittivo del test positivo diminuisce passando dal 94.1% al 24.1%

                 est       lower        upper

pv.pos    0.9408602   0.9214177    0.9567038  [primo esempio]

pv.pos    0.24137931  0.13870120   0.37168972 [secondo esempio]

e il valore predittivo del test negativo aumenta passando da 88.3% a 99.7%

                 est       lower        upper

pv.neg    0.8831776   0.8597456    0.9039310  [primo esempio]

pv.neg    0.99736148  0.99050154   0.99968030 [secondo esempio]

Mentre sensibilità e specificità non cambiano al cambiare della popolazione alla quale un test diagnostico viene applicato, in quanto sono caratteristiche intrinseche al test, la prevalenza della malattia non dipende dal test, ma è una caratteristica della popolazione alla quale il test viene applicato, e può cambiare notevolmente. Per esempio la prevalenza della malaria in Italia è diversa da quella che si ha nei paesi dell'Africa equatoriale; la prevalenza della malattia reumatica nei pazienti che accedono all'ambulatorio di un medico generico è differente dalla prevalenza della malattia reumatica nei pazienti ricoverati in un Reparto ospedaliero specializzato nella cura di questa malattia; la prevalenza del cancro al polmone è differente nei soggetti fumatori rispetto ai soggetti non fumatori; e così via. E questo cambia il valore predittivo del test in quanto al diminuire della prevalenza della malattia, cioè della probabilità a priori di malattia, diminuisce il valore predittivo del test positivo (la probabilità a posteriori di malattia in un soggetto cui il test è risultato positivo), e aumenta il valore predittivo del test negativo (la probabilità a posteriori di non-malattia in un soggetto cui il test è risultato negativo) [7].

Questo è un punto cruciale, e rende conto del perché in molte situazioni di bassa prevalenza un test sia utile, e venga impiegato, soprattutto per escludere una malattia. Mentre il teorema di Bayes mostra (e non solo in campo medico) come e quanto l'incertezza delle nostre conclusioni può essere ridotta grazie all'informazione/evidenza fornita dai fatti, e R fornisce opportuni strumenti di calcolo.

----------

[1] Abbiamo anche visto in un altro post come sovrapporre due Curve ROC sovrapposte.

[2] Gerhardt W, Keller H. Evaluation of test data from clinical studies. Scand J Clin Lab Invest, 46, 1986 (supplement 181).

https://www.tandfonline.com/doi/pdf/10.1080/00365518709168461

[3] Package ‘epiR’.

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

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

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

[6] Deeks JJ, Altman DG. Diagnostic tests 4: likelihood ratios. BMJ 2004;329;168-169.

https://www.bmj.com/content/329/7458/168.full

[7] Per i temi trattati nei post sul teorema di Bayes vedere Altman, Federspil, Galen, Gerhardt, Gigerenzer, Gill, Motterlini e Nordenström nella pagina Bibliografia.

APPENDICE

Si vuole dimostrare che i risultati ottenuti sopra con i conteggi dei casi sono identici a quelli ottenuti impiegando il teorema di Bayes con i valori di probabilità e le espressioni (iv) per il valore predittivo del test positivo e (v) per il valore predittivo del test negativo (vedi post Teorema di Bayes e diagnosi medica).

Primo esempio riportato sopra: abbiamo sensibilità del test 0.875, specificità del test 0.945 e prevalenza della malattia 0.50.

Valore predittivo del test positivo

                              P(T+|M+) · P(M+)

P(M+|T+) = ———————————————————— (iv)

                   P(T+|M+) · P(M+) + P(T+|M-) · P(M-)

che si legge come

                                        sensibilità · prevalenza

P(M+|T+) = ——————————————————————————————

                   sensibilità · prevalenza + (1 – specificità) · (1 - prevalenza)

quindi il valore predittivo del test positivo è uguale a

                              0.875 · 0.50

P(M+|T+) = ———————————————— = 0.9408602 

                   (0.875 · 0.50) + (0.055 · 0.50)

Valore predittivo del test negativo

                               P(T-|M-) · P(M-)

P(M-|T-) = ——————————————————— (v)

                 P(T-|M-) · P(M-) + P(T-|M+) · P(M+)

che si legge come

                                   specificità · (1 – prevalenza)

P(M-|T-) = ——————————————————————————————

                 specificità · (1 – prevalenza) + (1 – sensibilità) · prevalenza

quindi il valore predittivo del test negativo è uguale a

                            0.945 · 0.50

P(M-|T-) = ———————————————— = 0.8831776 

                 (0.945 · 0.50) + (0.125 · 0.50)

Secondo esempio riportato sopra: abbiamo sensibilità del test 0.875, specificità del test 0.945 e prevalenza della malattia 0.01969784.

Valore predittivo del test positivo

                                     0.875 · 0.01960784

P(M+|T+) = ———————————————————————— = 0.24137928 

                   (0.875 · 0.01960784) + (0.055 · 0.98039216)

Valore predittivo del test negativo

                                    0.945 · 0.98039216

P(M-|T-) = ———————————————————————— = 0.99736148 

                 (0.945 · 0.98039216) + (0.125 · 0.01960784)

I risultati del valore predittivo ottenuti impiegando i conteggi dei casi sono identici a quelli ottenuti impiegando il teorema di Bayes con i valori di probabilità, come volevasi dimostrare.

----------

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.