mercoledì 18 settembre 2019

Efficienza e accuratezza diagnostica

Nel post Teorema di Bayes e diagnosi medica abbiamo visto che è possibile impiegare il teorema di Bayes 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?

E abbiamo visto come un modo semplice per farlo si basi sul numero di casi osservati, organizzati come in questa tabella


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

nella quale sono definiti come:
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.

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 ricavare il valore soglia (Fig. 1d) 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). Fissato il valore soglia, risultano determinate la sensibilità e la specificità del test.

Nel post Sensibilità, specificità e valore predittivo abbiamo impiegato il pacchetto epiR per calcolare le probabilità che rispondono alle due domande: 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, e 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.

In questo post, che conclude la saga bayesiana, impieghiamo R per completare quanto visto finora con una analisi numerica e una analisi grafica dettagliate [2] dei dati del file bayes.csv [3]. Nel primo passaggio, impiegando il pacchetto ROCR, sono calcolati, per ciascuno dei possibili valori soglia compresi tra il valore minimo e il valore massimo osservati, il numero di VP, FP, FN e VN. Nel secondo passaggio, impiegando semplicemente le funzioni base di R, sono calcolati e rappresentati graficamente, per ciascuno dei possibili valori soglia compresi tra il valore minimo e il valore massimo osservati, sensibilità, specificità, efficienza diagnostica, valore predittivo del test positivo, valore predittivo del test negativo e accuratezza diagnostica.

Per proseguire è necessario:
effettuare il download del file bayes.csv 
salvare il file nella cartella C:\Rdati\

Inoltre il pacchetto aggiuntivo ROCR [4] deve essere scaricato dal CRAN. Copiate questo script, incollatelo nella Console di R e premete ↵ Invio:

# EFFICIENZA E ACCURATEZZA DIAGNOSTICA
#
library(ROCR) # carica il pacchetto
mydata <- read.table("c:/Rdati/bayes.csv", header=TRUE, sep=";") # importa i dati
attach(mydata) # consente di impiegare direttamente i nomi delle variabili
#
# calcola sensibilità, specificità, valore predittivo
#
pred <- prediction(valore, classe) # genera un oggetto che contiene i dati pre-elaborati
#
cut <- sort(unlist(pred@cutoffs), decreasing=FALSE) # valori soglia ordinati in ordine crescente
vp <- sort(unlist(pred@tp), decreasing=TRUE) # veri positivi per ciascun valore soglia ordinati in ordine decrescente
fp <- sort(unlist(pred@fp), decreasing=TRUE) # falsi positivi per ciascun valore soglia ordinati in ordine decrescente
fn <- sort(unlist(pred@fn), decreasing=FALSE) # falsi negativi per ciascun valore soglia ordinati in ordine crescente
vn <- sort(unlist(pred@tn), decreasing=FALSE) # veri negativi per ciascun valore soglia ordinati in ordine crescente
#
tabpred <- data.frame(cut, vp, fp, fn, vn) # genera la tabella dei valori calcolati
names(tabpred) <- c("cutoff", "VP", "FP", "FN", "VN") # assegna i nomi alle variabili/colonne
tabpred # mostra la tabella
#
sens <- vp/(vp+fn) # calcola sensibilità per ciascun valore soglia
spec <- vn/(vn+fp) # calcola specificità per ciascun valore soglia
effi <- sens*spec # calcola efficienza diagnostica per ciascun valore soglia
vptp <- vp/(vp+fp) # calcola valore predittivo del test positivo per ciascun valore soglia
vptn <- vn/(vn+fn) # calcola valore predittivo del test negativo per ciascun valore soglia
accu <- (vp+vn)/(vp+fp+fn+vn) # calcola accuratezza diagnostica per ciascun valore soglia
#
tabprob <- data.frame(cut, sens, spec, effi, vptp, vptn, accu) # genera la tabella dei valori calcolati
names(tabprob) <- c("cutoff", "sensibilità", "specificità", "efficienza", "vptp", "vptn", "accuratezza") # assegna i nomi alle variabili/colonne
tabprob # mostra la tabella
#
paste("Massima efficienza diagnostica ", round(max(effi), digits=3), " per un valore soglia di ", cut[which.max(effi)], sep="") # riporta la massima efficienza diagnostica e il corrispondente valore soglia
paste("Massima accuratezza diagnostica ", round(max(accu), digits=3), " per un valore soglia di ", cut[which.max(accu)], sep="") # riporta la massima accuratezza diagnostica e il corrispondente valore soglia
#
# traccia istogrammi con curve di sensibilità, specificità ed efficienza diagnostica
#
sani <- subset(valore, classe==0) # sottoinsieme con i valori dei soggetti sani per tracciare gli istogrammi
malati <- subset(valore, classe==1) # sottoinsieme con i valori dei soggetti malati per tracciare gli istogrammi
#
windows() # apre una nuova finestra
par(mar=c(5,4,4,5)+.1) # margini per fare spazio alla scala dell'asse delle y sulla destra
#
hist(malati, xlim=c(0,250), ylim=c(0,150), xaxp = c(0, 250, 5), yaxp = c(0, 150, 3), breaks=seq(from=0, to=250, by=5), border="red", col="red", density = 20, angle = 45, main="Sensibilità, specificità, efficienza diagnostica", xlab="Risultato del test", ylab="Frequenza", freq=TRUE, plot=TRUE) # traccia istogramma malati
hist(sani, xlim=c(0,250), breaks=seq(from=0, to=250, by=5), border="green4", col="lightgreen", density = 20, angle = -45, freq=TRUE, add=TRUE, plot=TRUE) # sovrappone istogramma sani
#
par(new=TRUE) # sovrappone curva sensibilità
plot(cut, sens, xlim=c(0,250), ylim=c(0,1), yaxp=c(0,1,5), type="l", lty=1, lwd=2, pch= 20, cex=0.5, col="blue", xaxt="n", yaxt="n", xlab="", ylab="") # traccia la curva
axis(4) # riporta un secondo asse delle y sulla destra per i valori di probabilità
mtext("Probabilità p", side=4, line=3) # aggiunge la legenda all'asse delle y di destra
#
par(new=TRUE) # sovrappone curva specificità
plot(cut, spec, xlim=c(0,250), ylim=c(0,1), yaxp=c(0,1,5), type="l", lty=1, lwd=2, pch=20, cex=0.5, col="mediumorchid", xaxt="n", yaxt="n", xlab="", ylab="") # traccia la curva
#
par(new=TRUE) # sovrappone curva efficienza diagnostica
plot(cut, effi, xlim=c(0,250), ylim=c(0,1), yaxp=c(0,1,5), type="l", lty=3, lwd=2, pch=20, cex=0.5, col="orange", xaxt="n", yaxt="n", xlab="", ylab="") # traccia la curva
#
legend(140, 0.97, legend=c("Sani", "Malati", "Sensibilità", "Specificità", "Efficienza diagnostica"), col=c("green", "red", "blue", "mediumorchid", "orange"), lty=c(1,1,1,1,3), lwd=c(1,1,2,2,2), cex=0.8) # aggiunge al grafico la legenda
#
# traccia istogrammi con curve di valore predittivo ed accuratezza diagnostica
#
windows() # apre una nuova finestra
par(mar=c(5,4,4,5)+.1) # margini per fare spazio alla scala dell'asse delle y sulla destra
#
hist(malati, xlim=c(0,250), ylim=c(0,150), xaxp = c(0, 250, 5), yaxp = c(0, 150, 3), breaks=seq(from=0, to=250, by=5), border="red", col="red", density = 20, angle = 45, main="Valore predittivo del test, accuratezza diagnostica", xlab="Risultato del test", ylab="Frequenza", freq=TRUE, plot=TRUE) # traccia istogramma malati
hist(sani, xlim=c(0,250), breaks=seq(from=0, to=250, by=5), border="green4", col="lightgreen", density = 20, angle = -45, freq=TRUE, add=TRUE, plot=TRUE) # sovrappone istogramma sani
#
par(new=TRUE) # sovrappone curva valore predittivo del test positivo
plot(cut, vptp, xlim=c(0,250), ylim=c(0,1), yaxp=c(0,1,5), type="l", lty=1, lwd=2, pch= 20, cex=0.5, col="blue", xaxt="n", yaxt="n", xlab="",ylab="") # traccia la curva
axis(4) # riporta l'asse delle y sulla destra
mtext("Probabilità p", side=4, line=3) # aggiunge legenda asse y
#
par(new=TRUE) # sovrappone curva valore predittivo del test negativo
plot(cut, vptn, xlim=c(0,250), ylim=c(0,1), yaxp=c(0,1,5), type="l", lty=1, lwd=2, pch=20, cex=0.5, col="mediumorchid", xaxt="n", yaxt="n", xlab="", ylab="") # traccia la curva
#
par(new=TRUE) # sovrappone curva accuratezza diagnostica
plot(cut, accu, xlim=c(0,250), ylim=c(0,1), yaxp=c(0,1,5), type="l", lty=3, lwd=2, pch=20, cex=0.5, col="orange", xaxt="n", yaxt="n", xlab="", ylab="") # traccia la curva
#
legend(135, 0.98, legend=c("Sani", "Malati", "Valore predittivo del test +", "Valore predittivo del test -", "Accuratezza diagnostica"), col=c("green", "red", "blue", "mediumorchid", "orange"), lty=c(1,1,1,1,3), lwd=c(1,1,2,2,2), cex=0.8) # aggiunge al grafico la legenda
#

Importare ed elaborare i vostri dati con questo script risulterà semplicissimo: dovete solamente strutturare i dati come nel file bayes.csv [5] quindi modificare opportunamente il nome e la posizione del file nella prima riga di codice.

Con le prime tre righe di codice viene caricato il pacchetto, sono importati i dati, e viene reso possibile con la funzione attach() impiegare direttamente i nomi delle variabili.

La riga di codice cruciale è la quarta, che salva nell'oggetto pred i dati generati dalla funzione prediction() del pacchetto ROCR e precisamente il vettore che contiene tutti i possibili valori soglia (pred@cutoffs) e i vettori che riportano per ciascun valore soglia:
→ i veri positivi (pred@tp)
→ i falsi positivi (pred@fp)
→ i falsi negativi (pred@fn)
→ i veri negativi (pred@tn).

Questi cinque vettori sono quindi combinati nella tabella tabpred che così vi apparirà:

    cutoff  VP  FP  FN  VN
1       19 843 853   0   0
2       22 843 852   0   1
3       26 843 851   0   2
4       28 843 850   0   3
          ..........
70      94 669   4 174 849
71      95 666   4 177 849
72      96 662   4 181 849
73      97 655   3 188 850
          .........
200    235   6   0 837 853
201    237   4   0 839 853
202    242   2   0 841 853
203    Inf   0   0 843 853

Pertanto fissando ad esempio il valore soglia (cutoff) tra sani e malati a 19 avremo 843 VP, 853 FP, 0 FN e 0 VN, fissando il valore soglia a 94 avremo 669 VP, 4 FP, 174 FN e 849 VN, e così via.

Nelle righe successive a partire dai vettori vp, fp, fn e vn sono calcolati per ciascuno dei valori soglia (cutoff) individuati:
→ sensibilità = VP/(VP+FN)
→ specificità = VN/(VN+FP)
→ efficienza diagnostica = sensibilità · specificità
valore predittivo del test positivo = VP/(VP+FP)
valore predittivo del test negativo = VN/(VN+FN)
→ accuratezza diagnostica = (VP+VN)/(VP+FP+FN+VN)
quindi i vettori sens, spec, effi, vptp, vptn, accu che contengono i dati calcolati sono combinati nella tabella tabprob:

    cutoff sensibilità specificità  efficienza      vptp      vptn accuratezza
1       19 1.000000000 0.000000000 0.000000000 0.4970519       NaN   0.4970519
2       22 1.000000000 0.001172333 0.001172333 0.4973451 1.0000000   0.4976415
3       26 1.000000000 0.002344666 0.002344666 0.4976387 1.0000000   0.4982311
4       28 1.000000000 0.003516999 0.003516999 0.4979327 1.0000000   0.4988208
                                    .........
70      94 0.793594306 0.995310668 0.789872879 0.9940565 0.8299120   0.8950472
71      95 0.790035587 0.995310668 0.786330848 0.9940299 0.8274854   0.8932783
72      96 0.785290629 0.995310668 0.781608140 0.9939940 0.8242718   0.8909198
73      97 0.776986951 0.996483001 0.774254289 0.9954407 0.8188825   0.8873821
                                    .........
200    235 0.007117438 1.000000000 0.007117438 1.0000000 0.5047337   0.5064858
201    237 0.004744958 1.000000000 0.004744958 1.0000000 0.5041371   0.5053066
202    242 0.002372479 1.000000000 0.002372479 1.0000000 0.5035419   0.5041274

203    Inf 0.000000000 1.000000000 0.000000000       NaN 0.5029481   0.5029481

Da questa tabella sono ricavati il valore soglia che corrisponde alla massima efficienza diagnostica

[1] "Massima efficienza diagnostica 0.858 per un valore soglia di 75"

e il valore soglia, identico al precedente, che corrisponde alla massima accuratezza diagnostica

[1] "Massima accuratezza diagnostica 0.927 per un valore soglia di 75"

Da notare nelle funzioni paste() che generano le due righe gli argomenti cut[which.max(effi)] e cut[which.max(accu)] nei quali la funzione which.max() restituisce la posizione, all'interno del rispettivo vettore, del valore massimo di efficienza o di accuratezza, posizione che viene impiegata per trovare nel vettore cut il valore soglia corrispondente. Scorrendo i dati della tabella è facile confermare che sia nel caso dell'efficienza diagnostica sia nel caso della accuratezza diagnostica il valore soglia “ottimale” tra sani e malati è quello che corrisponde al massimo valore trovato di efficienza o di accuratezza.

Dalla tabella sono infine generati i grafici che riportano le sei grandezze calcolate in funzione del valore soglia (cutoff). In questa parte l'unica cosa un po' particolare è la funzione subset() [6] che viene impiegata per generare, dai dati importati, due sottoinsiemi separati, uno contenente solamente i valori dei soggetti sani e l'altro contenente solamente i valori dei soggetti malati, 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).

Le funzioni e gli argomenti impiegati per costruire gli istogrammi con le distribuzioni dei risultati nel test nei sani e nei malati, e per sovrapporre ad essi le curve che rappresentano, in funzione dei valori soglia riportati in ascisse, la sensibilità, la specificità e l'efficienza diagnostica, sono illustrati sia nella documentazione di R delle funzioni hist() e plot() sia nei vari post dedicati alla rappresentazione grafica dei dati [7], ai quali si rimanda. Questo è il grafico risultante










fare click per ingrandire









Questo invece è il grafico che sovrappone agli istogrammi delle distribuzioni dei risultati nel test nei sani e nei malati le curve che rappresentano, in funzione dei valori soglia riportati in ascisse, il valore predittivo del test positivo, il valore predittivo del test negativo e l'accuratezza diagnostica









fare click per ingrandire










La conclusione è che un'analisi completa e dettagliata, sia numerica sia grafica, delle caratteristiche di un test per la diagnosi di una malattia, esempio paradigmatico di applicazione del teorema di Bayes, può essere realizzata a partire semplicemente dai conteggi di VP, FP, FN e VN effettuati per una serie opportuna di valori soglia.


----------

[1] Abbiamo anche visto in un altro post come realizzare 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). URL consultato il 18/06/2019: http://bit.ly/2IVp4ox

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

[4] Package ‘ROCR’. URL consultato il 15/06/2019: http://bit.ly/2MOmUvT

[5] Vedere il post Curve Roc e valore soglia.

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

[7] Per la funzione legend() consultare anche il post Aggiungere la legenda a un grafico.

domenica 14 luglio 2019

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, e 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.
come già illustrato nel post Teorema di Bayes e diagnosi medica.

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
                 est      lower       upper
aprev      0.4650000  0.4403314   0.4897973
tprev      0.5000000  0.4752070   0.5247930
se         0.8750000  0.8500684   0.8971267
sp         0.9450000  0.9268653   0.9597561
diag.acc   0.9100000  0.8949017   0.9235730
diag.or  120.2727273 83.1451376 173.9792528
nnd        1.2195122  1.1670208   1.2871111
youden     0.8200000  0.7769337   0.8568828
ppv        0.9408602  0.9214177   0.9567038
npv        0.8831776  0.8597456   0.9039310
plr       15.9090909 11.9229352  21.2279249
nlr        0.1322751  0.1100336   0.1590125

Riprendiamo le statistiche calcolate riportando (tra parentesi) il loro significato, 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%:

aprev (prevalenza apparente, soggetti con il test positivo)
(VP+FP)/(VP+FP+FN+VN) = 0.4650000 (0.4403314 - 0.4897973)

tprev (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.acc (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)

nnd (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)

ppv (valore predittivo di un test positivo)
VP/(VP+FP) = 0.9408602 (0.9214177 - 0.9567038)

npv (valore predittivo di un test negativo)
VN/(VN+FN) = 0.8831776 (0.8597456 - 0.9039310)

plr (rapporto di verosimiglianza LR+ per un test positivo) [6]
(VP/(VP+FN))/(FP/(FP+VN)) = 15.9090909 (11.9229352 - 21.2279249)

nlr (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)

> mystat <- epi.tests(mymatrix, conf.level=0.95) # calcola le statistiche
> summary(mystat) # mostra le statistiche
                  est       lower        upper
aprev      0.07107843  0.05441190   0.09091933
tprev      0.01960784  0.01124809   0.03164703
se         0.87500000  0.61652376   0.98448640
sp         0.94500000  0.92686534   0.95975609
diag.acc   0.94362745  0.92551979   0.95843621
diag.or  120.27272727 26.50441692 545.77804789
nnd        1.21951220  1.05905000   1.84030189
youden     0.82000000  0.54338911   0.94424248
ppv        0.24137931  0.13870120   0.37168972
npv        0.99736148  0.99050154   0.99968030
plr       15.90909091 11.30365876  22.39090713
nlr        0.13227513  0.03617550   0.48366194


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

                 est       lower        upper
tprev      0.5000000   0.4752070    0.5247930 [primo esempio]
tprev     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
ppv        0.9408602   0.9214177    0.9567038 [primo esempio]
ppv       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
npv        0.8831776   0.8597456    0.9039310 [primo esempio]
npv       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 realizzare 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). URL consultato il 18/06/2019: http://bit.ly/2IVp4ox

[3] Package ‘epiR’. URL consultato il 25/01/2019: https://goo.gl/wyuNsr

[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. URL consultato il 29/10/2018: https://goo.gl/ahpbpN

[7] Per i temi trattati nei tre post sul teorema di Bayes vedere Altman, Federspil, Galen, Gerhardt, Gigerenzer, Gill, Motterlini e Nordenström nella 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.

----------