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;

→ 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\

Trovate il file, con le istruzioni per scaricare questo e gli altri file di dati impiegati nei post, alla pagina Dati.

 

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 mediante la  funzione prediction() un oggetto che contiene una pre-elaborazione dei dati

#

cut <- sort(unlist(pred@cutoffs), decreasing=FALSE) # calcola i valori soglia 

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

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

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

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

#

tabpred <- data.frame(cut, vp, fp, fn, vn) # genera la tabella con i 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 la sensibilità per ciascun valore soglia

spec <- vn/(vn+fp) # calcola la specificità per ciascun valore soglia

effi <- sens*spec # calcola la efficienza diagnostica per ciascun valore soglia

vptp <- vp/(vp+fp) # calcola il valore predittivo del test positivo per ciascun valore soglia

vptn <- vn/(vn+fn) # calcola il valore predittivo del test negativo per ciascun valore soglia

accu <- (vp+vn)/(vp+fp+fn+vn) # calcola la accuratezza diagnostica per ciascun valore soglia

#

tabprob <- data.frame(cut, sens, spec, effi, vptp, vptn, accu) # genera la tabella con i 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.7) # 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.7) # 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

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

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

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

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

https://cran.r-project.org/web/packages/ROCR/ROCR.pdf

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

Grafici a linee (line chart)

Vediamo come rappresentare i grafici a linee impiegando una funzione e un set di dati inclusi nella installazione base di R. Si tratta della funzione plot() [1], utile per una molteplicità di rappresentazioni

Se nella Console di R digitate data() potete vedere l'elenco dei set dati che sono caricati automaticamente all'avvio del programma e che include il set di dati WorldPhones qui impiegato.

Copiate e incollate nella Console di R questo script e premete ↵ Invio:

# GRAFICI A LINEE set di dati di R

#

mydata <- data.frame(WorldPhones) # trasforma il set di dati WorldPhones in dataframe

mydata # mostra i dati

#

windows() # apre una nuova finestra

#

# traccia il grafico a linee per la prima variabile

plot(mydata$N.Amer, type="o", pch=19, lty=1, col="red", ylim=c(0,80000), axes=FALSE, ann=FALSE)

# sovrappone i grafici a linee per le variabili successive

lines(mydata$Europe, type="o", pch=20, lty=2, col="blue")

lines(mydata$Asia, type="o", pch=21, lty=3, col="brown")

lines(mydata$S.Amer, type="o", pch=22, lty=4, col="darkolivegreen")

lines(mydata$Oceania, type="o", pch=23, lty=5, col="gold1")

lines(mydata$Africa, type="o", pch=24, lty=6, col="green4")

lines(mydata$Mid.Amer, type="o", pch=25, lty=7, col="cornflowerblue")

#

# traccia gli assi con le rispettive scale e riporta le etichette sull'asse orizzontale

axis(1, at=1:7, labels=c("1951","1956","1957","1958","1959","1960","1961"))

axis(2, las=1, at=20000*0:80000)

#

# riporta titolo ed etichetta dell'asse delle x

title(main="Numero di apparecchi telefonici nel mondo", col.main="black", font.main=1)

title(xlab="Anno della rilevazione", col.lab="black")

#

# riporta una legenda con i nomi delle variabili e i simboli delle linee impiegati

legend(5, 70000, c("Nord America", "Europa", "Asia", "Sud America", "Oceania", "Africa", "Centro America"), cex=0.8, col=c("red","blue", "brown", "darkolivegreen", "gold1", "green4", "cornflowerblue"), pch=19:25, lty=1:7, bty="y")

#

Nella prima riga di codice dal set di dati WorldPhones con la funzione data.frame() viene creata una tabella (dataframe) che viene assegnata (<-) all'oggetto mydata, quindi con mydata sono mostrati i dati importati:

> mydata # mostra i dati

     N.Amer Europe Asia S.Amer Oceania Africa Mid.Amer

1951  45939  21574 2876   1815    1646     89      555

1956  60423  29990 4708   2568    2366   1411      733

1957  64721  32510 5230   2695    2526   1546      773

1958  68484  35218 6662   2845    2691   1663      836

1959  71799  37598 6856   3000    2868   1769      911

1960  76036  40341 8220   3145    3054   1905     1008

1961  79831  43173 9053   3338    3224   2005     1076

che sono rappresentati dalle rilevazioni, effettuate tra l'anno 1951 e l'anno 1961, del numero di apparecchi telefonici, suddivisi per le principali aree geografiche del mondo.

Dopo avere aperto la finestra grafica con windows() nel successivo blocco di codice viene generato mediante la funzione plot() il grafico a linee che riporta in ascisse (x) i descrittori delle righe (1951, 1952, … 1961) e impiega come argomenti:

→ mydata$N.Amer la prima variabile del set di dati, posta in ordinate (y);

→ type="o" che consente di sovrascrivere i grafici successivi;

→ pch=19 quale dei simboli dei punti di R deve essere impiegato per i punti da rappresentare;

→ lty=1 quale degli stili delle linee di R deve essere impiegato per la linea da rappresentare;

→ col="red" quale colore impiegare per i punti e le linee da rappresentare;

→ ylim=c(0,80000) limiti inferiore e superiore dell'asse delle y;

→ axes=FALSE che indica di non rappresentare gli assi del grafico che verranno poi configurati manualmente nelle righe successive;

→ ann=FALSE che indica di non rappresentare titolo ed etichette dell'asse delle x e delle y che verranno poi configurati manualmente nelle righe successive.

Al grafico così creato, con le sei righe di codice successive sono sovrapposti i grafici a linee delle altre sei variabili che vogliamo rappresentare: Europe, Asia, S.Amer, Oceania, Africa, Mid.Amer.

Viene poi tracciato l'asse delle x axis(1....) ponendo nelle posizioni che vanno dalla 1 alla 7 (at=1:7) le etichette "1951", "1956", "1957", "1958", "1959", "1960", "1961" riportate nell'argomento labels. Per l'asse delle ordinate axis(2....) le tacche sono tracciate in orizzontale (las = 1) ogni 20000 (at=20000) per la scala che va da 0 a 80000 (0:80000).

Titolo ed etichetta dell'asse delle x riportati nel blocco di codice successivo non richiedono particolari commenti.

Infine la funzione legend() consente di aggiungere una legenda, che prevede i seguenti argomenti:

→ 5, 70000 sono le coordinate x e y alle quali viene posizionato l'angolo superiore sinistro della legenda;

→ "Nord America", "Europa", "Asia", "Sud America", "Oceania", "Africa", "Centro America" sono i nomi da riportare;

→ cex.axis = 0.8 specifica la dimensione dei caratteri da impiegare;

→ "red","blue", "brown", "darkolivegreen", "gold1", "green4", "cornflowerblue" sono i colori di R che ovviamente riprendono nell'ordine quelli delle linee riportate nel grafico;

→ pch=19:25 sono i simboli dei punti di R impiegati nel grafico;

→ lty=1:7 sono gli stili delle linee di R impiegati nel grafico;

→ bty="y" che è l'opzione di default della funzione plot() e consente di tracciare il riquadro contenente la legenda, viene qui riportato per ricordare che con bty="n" è possibile eliminarlo.

Da notare che l'esempio è stato sviluppato esclusivamente per illustrare le funzioni e gli argomenti per rappresentare un grafico a linee. Ma si fa notare nel set di dati fornito, tra il 1951 e il 1956 intercorrono 5 anni, mentre gli intervalli successivi sono di un anno: questo distorce la rappresentazione, aumentando la pendenza del primo segmento di retta rispetto al reale. Ovviamente si raccomanda di non introdurre mai distorsioni di questo o di altro genere nella rappresentazione dei propri dati.

I grafici a linee vengono spesso impiegati per rappresentare variabili con pochi dati: potrebbe pertanto essere utile disporre di un esempio nel quale i dati sono inseriti manualmente. Copiate e incollate nella Console di R questo script e premete ↵ Invio:

# GRAFICI A LINEE dataframe costruito manualmente

#

N.Amer <- c(45939, 60423, 64721, 68484, 71799, 76036, 79831) # vettore con la prima variabile

Europe <- c(21574, 29990, 32510, 35218, 37598, 40341, 43173) # vettore con la seconda variabile

Asia <- c(2876, 4708, 5230, 6662, 6856, 8220, 9053) # vettore con la terza variabile

S.Amer <- c(1815, 2568, 2695, 2845, 3000, 3145, 3338) # vettore con la quarta variabile

Oceania <- c(1646, 2366, 2526, 2691, 2868, 3054, 3224) # vettore con la quinta variabile

Africa <- c(89, 1411, 1546, 1663, 1769, 1905, 2005) # vettore con la sesta variabile

Mid.Amer <- c(555, 733, 773, 836, 911, 1008, 1076) # vettore con la settima variabile

mydata <- data.frame(N.Amer, Europe, Asia, S.Amer, Oceania, Africa, Mid.Amer) # combina i vettori nel dataframe mydata

row.names(mydata) <- c(1951, 1956, 1957, 1958, 1959, 1960, 1961) # aggiunge i nomi delle righe

mydata # mostra i dati

#

windows() # apre una nuova finestra

#

# traccia il grafico a linee per la prima variabile

plot(N.Amer, type="o", pch=19, lty=1, col="red", ylim=c(0,80000), axes=FALSE, ann=FALSE)

# sovrappone i grafici a linee per le variabili successive

lines(Europe, type="o", pch=20, lty=2, col="blue")

lines(Asia, type="o", pch=21, lty=3, col="brown")

lines(S.Amer, type="o", pch=22, lty=4, col="darkolivegreen")

lines(Oceania, type="o", pch=23, lty=5, col="gold1")

lines(Africa, type="o", pch=24, lty=6, col="green4")

lines(Mid.Amer, type="o", pch=25, lty=7, col="cornflowerblue")

#

# traccia gli assi con le rispettive scale e riporta le etichette sull'asse orizzontale

axis(1, at=1:7, labels=c("1951", "1956", "1957", "1958", "1959", "1960", "1961"))

axis(2, las=1, at=20000*0:80000)

#

# riporta titolo ed etichetta dell'asse delle x

title(main="Numero di apparecchi telefonici nel mondo", col.main="black", font.main=1)

title(xlab="Anno della rilevazione", col.lab="black")

#

# riporta una legenda con i nomi delle variabili e i simboli delle linee impiegati

legend(5, 70000, c("Nord America", "Europa", "Asia", "Sud America", "Oceania", "Africa", "Centro America"), cex=0.8, col=c("red","blue", "brown", "darkolivegreen", "gold1", "green4", "cornflowerblue"), pch=19:25, lty=1:7, bty="y")

#

Dopo avere inserito mediante la funzione c() i valori delle singole variabili, questi sono combinati mediante la funzione data.frame() nella tabella (dataframe) mydata, alla quale sono infine aggiunti con la funzione row.names() i nomi delle righe che verranno riportati sull'asse orizzontale del grafico. Il resto del codice è identico a quello dello script precedente.

Da notare come ultima cosa che in quest'ultimo script avendo eseguito la funzione row.names() si possono impiegare direttamente i nomi delle variabili (ad esempio N.Amer) senza eseguire la funzione attach(), mentre nello script precedente, nel quale la funzione attach() non è stata eseguita, si è reso necessario specificare per le variabili il nome completo (ad esempio mydata$N.Amer).

Entrambi gli script possono essere facilmente riutilizzati, il primo sostituendo al set di dati WorldPhones i propri dati, organizzati in modo analogo, il secondo adattando opportunamente numero, nomi e valori delle variabili che confluiscono nell'oggetto mydata.

Nota bene: in alternativa alle funzioni di base qui illustrate è possibile realizzare grafici a linee anche impiegando le funzioni del pacchetto ggplot2, per questo rimando al post Grafici a linee con ggplot.

----------

[1] Digitate help(plot) nella Console di R per la documentazione della funzione plot() e digitate help(nomedellafunzione) nella Console di R per la documentazione delle altre funzioni impiegate nello script.