Inserire più grafici nella stessa immagine

Vediamo come sia possibile disporre in vari modi più grafici all'interno di una stessa immagine impiegando una sola riga di codice.

Per farlo ci serviamo delle funzioni par(), layout() e matrix(), impiegando come esempi alcuni grafici realizzati con la funzione boxplot(), per i cui dettagli si rimanda al post Grafici a scatola con i baffi (boxplot) e alla documentazione di R [1].

Accertatevi di avere installato il pacchetto DAAG che contiene il set di dati ais. Altrimenti dovete scaricarlo e installarlo dal CRAN. Tutte le funzioni qui impiegate fanno invece parte del pacchetto base di funzioni installato con R.

Copiate lo script che segue, incollatelo nella Console di R e premete ↵ Invio:

# INSERIRE PIU' GRAFICI NELLA STESSA IMMAGINE

#

library(DAAG) # carica il pacchetto DAAG incluso il set di dati ais

attach(ais) # consente di impiegare direttamente i nomi delle variabili

#

par(mfrow=c(2,2)) # (1) quattro grafici in due righe per due colonne

boxplot(wt~sex, data=ais, xlab="Sesso", ylab="Peso corporeo in kg", notch=TRUE, col="green") # boxplot per sesso del peso corporeo

boxplot(ht~sex, data=ais, xlab="Sesso", ylab="Altezza in cm", notch=TRUE, col="green") # traccia i boxplot per sesso del peso corporeo

boxplot(bmi~sex, data=ais, xlab="Sesso", ylab="BMI in kg/m2", notch=TRUE, col="green") # boxplot per sesso della percentuale di grasso corporeo

boxplot(pcBfat~sex, data=ais, xlab="Sesso", ylab="Grasso corporeo in %", notch=TRUE, col="yellow") # boxplot per sesso della percentuale di grasso corporeo

#

windows() # apre una nuova finestra

par(mfrow=c(1,2)) # (2) due grafici in una riga e due colonne:

boxplot(bmi~sex, data=ais, xlab="Sesso", ylab="BMI in kg/m2", notch=TRUE, col="green") # boxplot per sesso della percentuale di grasso corporeo

boxplot(pcBfat~sex, data=ais, xlab="Sesso", ylab="Grasso corporeo in %", notch=TRUE, col="yellow") # boxplot per sesso della percentuale di grasso corporeo

#

windows() # apre una nuova finestra

layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE)) # (3) tre grafici, il primo occupa riga 1 / colonne 1 e 2, il secondo riga 2 / colonna 1, il terzo riga 2 / colonna 2

boxplot(pcBfat~sport, horizontal=FALSE, boxwex = 0.4, cex.axis = 0.8, las = 2, data=ais, main="Percentuale di grasso corporeo nei vari sport", xlab="Sport praticato", ylab="Grasso corporeo in %", notch=FALSE, col="yellow") # traccia i boxplot per sport praticato

boxplot(wt~sex, data=ais, xlab="Sesso", ylab="Peso corporeo in kg", notch=TRUE, col="green") # boxplot per sesso del peso corporeo

boxplot(ht~sex, data=ais, xlab="Sesso", ylab="Altezza in cm", notch=TRUE, col="green") # traccia i boxplot per sesso del peso corporeo

#

windows() # apre una nuova finestra

layout(matrix(c(1,3,2,3), 2, 2, byrow = TRUE)) # (4) tre grafici, il primo occupa riga 1 / colonna 1, il secondo riga 2 / colonna 1, il terzo colonna 2 / righe 1 e 2

boxplot(wt~sex, data=ais, xlab="Sesso", ylab="Peso corporeo in kg", notch=TRUE, col="green") # boxplot per sesso del peso corporeo

boxplot(ht~sex, data=ais, xlab="Sesso", ylab="Altezza in cm", notch=TRUE, col="green") # traccia i boxplot per sesso del peso corporeo

boxplot(pcBfat~sex, data=ais, xlab="Sesso", ylab="Grasso corporeo in %", notch=TRUE, col="yellow") # boxplot per sesso della percentuale di grasso corporeo

#

Nel primo blocco di codice con par(mfrow=c(2,2)) si predispone la suddivisione della finestra grafica in 2 righe per 2 colonne quindi in quattro quadranti

1	2
3	4

nei quali sono inseriti, procedendo da sinistra verso destra e dall'altro verso il basso, i quattro grafici realizzati con la funzione boxplot() con questo risultato

fare click sull'immagine per ingrandirla










Nel secondo blocco di codice con par(mfrow=c(1,2)) si predispone la suddivisione della finestra in 1 riga e 2 colonne

1

2

per realizzare due grafici affiancati, il primo sulla sinistra e il secondo sulla destra, con questo risultato

fare click sull'immagine per ingrandirla










Al bisogno potete impiegare par(mfrow=c(2,1)) per predisporre la suddivisione della finestra in 2 righe e 1 colonna

1

2

(immagine non riportata) ma anche aumentare il numero delle righe e/o delle colonne, prestando ovviamente sempre molta attenzione alla leggibilità e alla fruibilità del risultato finale.

Nel terzo blocco di codice con layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE)) vediamo come organizzare i grafici all'interno dell'immagine in maniera più articolata. Viene infatti predisposta la suddivisione della finestra in 2 righe per 2 colonne (, 2, 2,) procedendo per riga (byrow = TRUE), ma specificando con il primo argomento c(1,1,2,3) che il primo grafico occuperà riga 1/colonna 1 e riga 1/colonna 2, il secondo grafico occuperà riga 2/colonna 1 e il terzo grafico occuperà riga 2/colonna 2, con questa logica

1     1
2	3

e questo risultato

fare click sull'immagine per ingrandirla










Nel quarto blocco di codice con layout(matrix(c(1,3,2,3), 2, 2, byrow = TRUE)) si predispone di nuovo la suddivisione della finestra in 2 righe per 2 colonne (, 2, 2,) procedendo per riga (byrow = TRUE), ma questa volta specificando con il primo argomento c(1,3,2,3) che il primo grafico occuperà riga 1/colonna 1, il secondo grafico occuperà riga 2/colonna 1 e il terzo grafico occuperà riga 1/colonna 2 e riga 2/colonna 2, con questa logica

1	3 3
2

e con questo risultato

fare click sull'immagine per ingrandirla

Infine, se lo volete, impiegando l'utilità riportata nel post Salvare i grafici di R in un file, potete anche salvare le immagini realizzate con R sotto forma di file .bmp, .jpeg, .png, .pdf, .ps per poterle stampare, archiviare, inserire in una pubblicazione, in un post o in un sito web.

----------

[1] Digitare help(nomedellafunzione) nella Console di R per la documentazione di questa e delle altre funzioni qui impiegate.

Analisi dei gruppi (clustering gerarchico e non gerarchico)

L'obiettivo dell'analisi dei gruppi (cluster analysis o clustering) è concettualmente semplice: verificare la possibile esistenza, in un insieme di oggetti, di sottoinsiemi/gruppi/cluster (di oggetti) particolarmente simili tra loro. L'analisi dei gruppi si applica a dati multivariati ed è un metodo statistico di tassonomia numerica che riveste un ruolo importante nella analisi esplorativa dei dati.

Nonostante alla base dell'analisi dei gruppi vi sia un'idea semplice e logica, esistono numerosi modi per realizzarla. Esistono metodi gerarchici, che danno luogo a una suddivisione ad albero (dendrogramma) in base alla distanza tra i singoli oggetti dell'insieme, e metodi non gerarchici, nei quali l'appartenenza di un oggetto (dell'insieme) ad uno specifico sottoinsieme/gruppo/cluster viene stabilita sulla base della sua distanza dal centro o dalla media dei dati o da un punto rappresentativo del cluster [1]. Nei metodi che procedono bottom-up all'inizio del processo di classificazione ad ogni oggetto viene fatto corrispondere un cluster. In questo stadio gli oggetti sono considerati tutti dissimili tra di loro. Al passaggio successivo i due oggetti più simili sono raggruppati nello stesso cluster. Il numero dei cluster risulta quindi pari al numero di oggetti diminuito di uno. Il procedimento viene ripetuto ciclicamente, fino ad ottenere (all'ultimo passaggio) un unico cluster. Nei metodi che procedono top-down invece inizialmente tutti gli oggetti sono considerati come appartenenti ad un unico cluster, che viene via via suddiviso in ulteriori cluster. E ancora esistono metodi esclusivi, che prevedono che un oggetto possa appartenere esclusivamente a un cluster, e metodi non esclusivi (fuzzy) che prevedono che un oggetto possa appartenere, in modo quantitativamente diverso, a più cluster. Da notare infine che non esiste una regola per stabilire il numero di cluster in cui sono aggregati gli oggetti, che va deciso caso per caso, cosa che complica ulteriormente lo scenario.

Qui ci occupiamo dell'implementazione con R dei due metodi più tradizionali, il metodo di Ward, un metodo gerarchico che misura la somiglianza di due oggetti/punti sulla base della loro distanza euclidea, e il metodo di MacQueen, un metodo non gerarchico che misura l'appartenenza di un oggetto/punto ad uno specifico cluster con l'algoritmo delle k-means. Entrambi, inoltre, sono metodi esclusivi. Come dati impieghiamo i valori di BMI (indice di massa corporea) rilevati a livello europeo alcuni anni fa e pubblicati dall'Istat [2].

Per proseguire è necessario:

→ effettuare il download del file bmi.csv

→ salvare il file nella cartella C:\Rdati\

In alternativa potete anche copiare i dati riportati qui sotto e salvarli in C:\Rdati\ in un file di testo denominato bmi.csv (attenzione, assicuratevi che il vostro editor di testo abbia salvato il file con l'estensione .csv).

Nazione;sottopeso;normale;sovrappeso;obeso

Austria;2.4;49.6;33.3;14.7

Belgio;2.7;48.0;35.3;14.0

Bulgaria;2.2;43.8;39.2;14.8

Cipro;3.9;47.8;33.8;14.5

Croazia;1.9;40.7;38.7;18.7

Danimarca;2.2;50.0;32.9;14.9

Estonia;2.2;43.9;33.5;20.4

Finlandia;1.2;44.1;36.4;18.3

Francia;3.2;49.6;31.9;15.3

Germania;1.8;46.1;35.2;16.9

Grecia;1.9;41.3;39.4;17.3

Irlanda;1.9;42.3;37.0;18.7

Lettonia;1.7;41.8;35.2;21.3

Lituania;1.9;42.5;38.3;17.3

Lussemburgo;2.8;49.3;32.4;15.6

Malta;2.0;37.0;35.0;26.0

Olanda;1.6;49.0;36.0;13.3

Polonia;2.4;42.9;37.5;17.2

Portogallo;1.8;44.6;36.9;16.6

Regno Unito;2.1;42.2;35.6;20.1

Repubblica Ceca;1.1;42.1;37.6;19.3

Romania;1.3;42.9;46.4;9.4

Slovacchia;2.1;43.6;38.0;16.3

Slovenia;1.6;41.8;37.4;19.2

Spagna;2.2;45.4;35.7;16.7

Svezia;1.8;48.3;35.9;14.0

Ungheria;2.9;41.9;34.0;21.2

Inoltre è necessario scaricare dal CRAN il pacchetto aggiuntivo cluster [3]. Copiate questo script, incollatelo nella Console di R e premete ↵ Invio:

# ANALISI DEI GRUPPI (CLUSTERING GERARCHICO E NON GERARCHICO)

#

mydata <- read.table("c:/Rdati/bmi.csv", header=TRUE, sep=";", row.names="Nazione") # importa i dati

mydata # mostra i dati importati

z <- scale(mydata) # standardizza le variabili

z # mostra le variabili standardizzate

#

# clustering gerarchico con il metodo di Ward (distanza euclidea)

#

mat <- hclust(dist(z, method = "euclidean"), method="ward.D") # matrice delle distanze euclidee

#

plot(mat, main="Dendrogramma", xlab="BMI nei Paesi europei", sub="", ylab="Distanza nei valori di BMI dei cluster", cex=0.7) # traccia il dendrogramma

rect.hclust(mat, k=4, border=c("red","blue","green","goldenrod")) # evidenzia i 4 gruppi/cluster

#

# clustering non gerarchico con il metodo di MacQueen (k-means )

#

myclust <- kmeans(z, 4, algorithm = c("MacQueen"), nstart=50) # genera i 4 gruppi/cluster

#

library(cluster) # carica il pacchetto

windows() # apre una nuova finestra

clusplot(z, myclust$cluster, color=TRUE, col.clus=c("blue","goldenrod","red","green"), shade=FALSE, labels=3, lines=0, main = "Grafico dei cluster", xlab="Componente principale 1", sub="", ylab="Componente principale 2", cex=0.6, col.txt="black", col.p="black") # grafico dei cluster per le prime due componenti principali

#

Con la prima riga di codice sono importati i dati, poi mostrati con mydata alla seconda riga. Con la funzione scale() della terza riga per ciascun dato viene calcolata la deviata normale standardizzata z. In pratica questa funzione prima calcola per i dati di ciascuna colonna/variabile la media e la deviazione standard, poi calcola per ciascuno dato x la corrispondente deviata normale standardizzata z come

z = (x – media) / deviazione standard

I valori di z sono poi qui salvati nel nuovo oggetto qui denominato, per comodità, z. Notare che quando il contenuto dell'oggetto viene mostrato alla quarta riga di codice, in coda ai dati compare

attr(,"scaled:center")

 sottopeso    normale sovrappeso      obeso 

  2.103704  44.537037  36.240741  17.111111 

attr(,"scaled:scale")

 sottopeso    normale sovrappeso      obeso 

 0.6111033  3.3717318  2.8989732  3.2322097

Quindi la funzione, oltre a riportare le deviate normali standardizzate z dei valori calcolate, riporta per completezza le medie ("scaled:center") e le deviazioni standard ("scaled:scale") delle quattro variabili/colonne, che sono state impiegate per calcolare i valori di z [4].

Dopo questi calcoli preliminari si passa al clustering gerarchico con il metodo di Ward con la funzione hclust() nella quale l'argomento method può assumere uno dei seguenti valori: "ward.D", "ward.D2", "single", "complete", "average" (= UPGMA), "mcquitty" (= WPGMA), "median" (= WPGMC), "centroid" (= UPGMC) [5]. All'interno della funzione hclust() con la funzione dist() viene calcolata la matrice delle distanze. Qui viene impiegato come metodo di calcolo il metodo di euclideo (method = "euclidean"), che prevede di misurare la distanza tra due punti con il teorema di Pitagora, ma l'elenco dei metodi che si possono impiegare include: "euclidean", "maximum", "manhattan", "canberra", "binary", "minkowski" [6].

Una volta salvata la matrice delle distanze euclidee tra tutte le possibili coppie di dati, nell'oggetto mat, con la funzione plot() viene tracciato il dendrogramma e con la funzione rect.hclust() e l'argomento k=4 (che può essere cambiato al bisogno) viene definito in 4 il numero dei gruppi da evidenziare nel dendrogramma.

fare click per ingrandire










Per il clustering non gerarchico con il metodo di MacQueen viene come prima cosa generata la matrice delle distanze con l'algoritmo k-means impiegando la funzione kmeans(). Quindi dopo avere caricato il pacchetto cluster e aperto una nuova finestra grafica con windows(), viene impiegata per tracciare il grafico la funzione clusplot() nella quale si fanno notare i seguenti argomenti:

→ color=TRUE consente l'impiego dei colori nei dati riportati sul grafico;

→ col.clus nel quale la sequenza dei colori è stata adattata manualmente per avere gli stessi colori del dendrogramma;

→ shade che se posto =TRUE riporta all'interni dei cluster un tratteggio con una densità uguale al numero di punti inclusi nel cluster diviso per l'area dell'ellisse;

→ labels=3 che è uno dei possibili valori di questo argomento (vederli con help(clusplot));

→ sub="" che elimina il sottotitolo previsto di default dalla funzione;

→ col.txt che definisce il colore del testo che compare all'interno del grafico;

→ col.p che definisce il colore impiegato per rappresentare i punti all'interno del grafico.

Da notare che se digitate myclust$cluster viene mostrato il vettore con l'elenco delle nazioni, ciascuna con il cluster cui appartiene

> myclust$cluster

        Austria          Belgio        Bulgaria           Cipro         Croazia 
              4               4               2               4               2 
      Danimarca         Estonia       Finlandia         Francia        Germania 
              4               3               2               4               2 
         Grecia         Irlanda        Lettonia        Lituania     Lussemburgo 
              2               2               3               2               4 
          Malta          Olanda         Polonia      Portogallo     Regno Unito 
              3               4               2               2               3 
Repubblica Ceca         Romania      Slovacchia        Slovenia          Spagna 
              2               1               2               2               2 
         Svezia        Ungheria 
              4               3 

Questo infine è il grafico del clustering non gerarchico con il metodo di MacQueen realizzato mediante la funzione clusplot()

fare click per ingrandire










Per un confronto contenente anche il dettaglio numerico i dati originali forniti dall'Istat sono stati importati in un foglio elettronico nel quale sono stati ordinati e suddivisi manualmente impiegando i seguenti criteri:

→ è stato identificato un paese che si discosta in modo significativo da tutti gli altri essendo l'unico con una quota di obesi inferiore al 10% e l'unico con una quota di soggetti sovrappeso superiore al 40%;

→ sono stati raccolti in un unico gruppo i paesi con una percentuale di soggetti di peso normale superiore del 5% circa alla media dei restanti paesi, e precisamente superiore al 47%.

→ sono stati raggruppati insieme tutti i paesi rimanenti, distinguendo tra quelli con una quota di obesi uguale o superiore al 10% e inferiore al 20% e quelli con una quota di obesi uguale o superiore al 20% e inferiore al 30%. 

Dalla tabella, che riconduce i risultati del clustering ai dati originari, risulta evidente che nel cluster rosso è dominante la presenza di peso normale, nel cluster blu è dominante la presenza di obeso, nel cluster giallo, collocato nella parte inferiore dell'obeso, è dominante la presenza di sovrappeso, mentre un caso a sé stante è quello della Romania, con valori estremi sia nel sovrappeso sia nell'obeso.

Questa analisi, semplice ma di buon senso, coincide con quella fornita dai due metodi di clustering, a parte il fatto che il metodo di MacQueen include Svezia e Olanda nel cluster normale cerchiato in rosso     mentre il metodo di Ward li include nel cluster sovrappeso cerchiato in giallo    . Qui vale la pena di notare tre cose:
→ la prima è che due metodi, impiegando differenti criteri per il raggruppamento dei dati in cluster, possono fornire risultati diversi;
→ la seconda cosa da notare è che Svezia e Olanda sembrano essere a cavallo tra due cluster, e questo potrebbe deporre a favore dell'impiego in questo caso del clustering con metodi non esclusivi (fuzzy) assegnando i due paesi, salomonicamente, contemporaneamente ai due cluster, anche se in modo quantitativamente diverso
→ la terza cosa da notare infine è che, come già accennato, non esiste una "regola aurea" per stabilire il numero di cluster in cui aggregare gli oggetti, anche se qui quattro cluster sembra essere una scelta adeguata.

Su questi temi vedremo di tornare. 

Nel frattempo per riutilizzare lo script dovete semplicemente:

→ nella prima riga di codice sostituire c:/Rdati/bmi.csv con posizione e nome del vostro file di dati, quindi togliere , row.names="Nazione" o, se le vostre osservazioni/righe hanno un identificativo, sostituire a Nazione il nome della variabile che lo contiene; 
→ nella funzione plot() correggere opportunamente i valori degli argomenti xlab e ylab. 

----------

[1] Massart DL, Vandeginste BGM, Deming SN, Michotte Y, Kaufman L. Chemometrics: a textbook. Elsevier, New York, 1988, ISBN 0-444-42660-4, pp. 371-384.

[2] Vedere il post Indice di massa corporea (BMI).

[3] Vedere il manuale di riferimento del pacchetto Package 'cluster'. URL consultato il 27/09/2019: http://bit.ly/2neaTna

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

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

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

Efficienza e accuratezza diagnostica [1]

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\

Il file contiene la variabile classe che indica la classe di appartenenza dell'osservazione/riga (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 la variabile quantitativa valore, che riporta il risultato del test, come illustrato in questo esempio che include le prime cinque e le ultime cinque osservazioni/righe del file:

classe;valore
0;19
0;22
0;26
0;28
0;29
.....
1;235
1;237
1;237
1;242
1;242

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 

vp <- sort(unlist(pred@tp), decreasing=TRUE) # veri positivi 

fp <- sort(unlist(pred@fp), decreasing=TRUE) # falsi positivi 

fn <- sort(unlist(pred@fn), decreasing=FALSE) # falsi negativi 

vn <- sort(unlist(pred@tn), decreasing=FALSE) # veri negativi 

#

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
          ..........
49      73 771  55  72 798
50      74 766  48  77 805
51      75 763  44  80 809
52      76 760  42  83 811
53      77 755  39  88 814

          .........

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 75 avremo 763 VP, 44 FP, 80 FN e 809 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
                                    .........
49      73 0.914590747 0.935521688 0.855619480 0.9334140 0.9172414   0.9251179
50      74 0.908659549 0.943728019 0.857527476 0.9410319 0.9126984   0.9262972
51      75 0.905100830 0.948417351 0.858413331 0.9454771 0.9100112   0.9268868
52      76 0.901542112 0.950762016 0.857151996 0.9476309 0.9071588   0.9262972
53      77 0.895610913 0.954279015 0.854662700 0.9508816 0.9024390   0.9251179

                                    .........

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"

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

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

e i valori delle grandezze (sensibilità, specificità, valori predittivi) corrispondenti al valore soglia 75 così individuato.

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, come magistralmente illustrato da Gerhardt e da Keller oltre trent'anni fa [2].

Per una versione lievemente più compatta del codice, che tabula in modo più compatto i dati e che impiega al posto degli istogrammi i kernel density plot, vedere nell'Indice il post Efficienza e accuratezza diagnostica [2].

Se strutturate i vostri dati seguendo le semplici indicazioni riportate sopra per il file bayes.csv per riutilizzare lo script dovete solamente, nella prima riga di codice, sostituire c:/Rdati/bayes.csv con posizione e nome del vostro file di dati.

----------

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