Grafici prima-dopo (slopegraph)

In un esempio Campbell [1] riporta i dati di VEMS (misura del Volume Espiratorio Massimo nel primo Secondo) effettuate in 5 soggetti asmatici prima (al tempo t0) e dopo (al tempo t1) l'assunzione di un broncodilatatore, dati riportati in questa tabella:

t0	t1	Differenza
1.5	1.7	-0.2
1.7	1.9	-0.2
2.1	2.2	-0.1
1.6	1.9	-0.3
2.4	2.4	0

In un caso come questo di dati appaiati può essere interessante integrare l'analisi statistica con una rappresentazione grafica che, se da un lato è semplicemente la versione minimalista dei grafici a linee, è in grado di fornire un significativo valore aggiunto all'analisi numerica.

La definizione di grafici prima-dopo è giustificata dal fatto che nel caso dei dati appaiati il disegno sperimentale tipicamente prevede che sullo stesso caso oggetto di indagine sia effettuata la stessa misura due volte, prima e dopo uno specifico trattamento. In questo modo la differenza tra la prima e la seconda misura è in linea di principio imputabile al solo trattamento deliberatamente previsto dal disegno sperimentale.

Nel pacchetto CGPfunctions che impiegheremo, che è necessario scaricare dal CRAN, seguendo la concettualizzazione del tema ad opera di Edward Tufte [2] il grafico viene denominato slopegraph, ponendo in questo modo l'enfasi sul fatto che la grafica è particolarmente utile a evidenziare la pendenza (slope) della linea che unisce i punti, che nel nostro caso prima-dopo sono solamente due, ma che in realtà possono essere molti.

Per effettuare la rappresentazione grafica i dati devono essere strutturati in modo differente da quello della tabella riportata da Campbell sulla pagina di un libro [1], ma devono essere invece strutturati nella forma più generale qui illustrata

Caso	Tempo	Valore
paziente_A	t0	1.5
paziente_A	t1	1.7
paziente_B	t0	1.7
paziente_B	t1	1.9
paziente_C	t0	2.1
paziente_C	t1	2.2
paziente_D	t0	1.6
paziente_D	t1	1.9
paziente_E	t0	2.4
paziente_E	t1	2.4

che corrisponde alla forma nella quale record e campi sono estratti da un database, dove a ogni dato numerico corrisponde un record, e i record possono essere aggregati tramite un identificatore univoco [3] qui riportato come variabile Caso e rappresentato da uno specifico paziente. Una struttura dati di questo genere segue il concetto di "atomizzazione dell'informazione" e consente non solo di ottenere una generalizzazione della soluzione grafica, ma soprattutto di impiegare direttamente i dati estratti da un database.

Per proseguire ora è necessario:

→ effettuare il download del file FEV1slopegraph.csv 

→ salvare il file nella cartella C:\Rdati\

Per il download del file vedere link e istruzioni riportati alla pagina Dati. In alternativa potete copiare le undici righe riportate qui sotto aggiungendo un ↵ Invio al termine dell'ultima riga e salvarle in C:\Rdati\ in un file di testo denominato FEV1slopegraph.csv (di default i file di testo sono salvati con l'estensione .txt ma noi qui salviamo il file con l'estensione .csv).

Caso;Tempo;Valore

paziente_A;t0;1.5

paziente_A;t1;1.7

paziente_B;t0;1.7

paziente_B;t1;1.9

paziente_C;t0;2.1

paziente_C;t1;2.2

paziente_D;t0;1.6

paziente_D;t1;1.9

paziente_E;t0;2.4

paziente_E;t1;2.4

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

# GRAFICO PRIMA-DOPO (SLOPEGRAPH) dati importati da un file

#

library(CGPfunctions) # carica il pacchetto necessario

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

mydata # mostra i dati

windows() # apre una nuova finestra

#

newggslopegraph(dataframe=mydata, Grouping=Caso, Times=Tempo, Measurement=Valore, Title="VEMS prima (t0) e dopo (t1) l'assunzione di un broncodilatatore", SubTitle ="", Caption=NULL, XTextSize=15, YTextSize=4, DataTextSize=4, LineThickness=2) # traccia il graficoprima-dopo (slopegraph)

#

Potete anche utilizzare in alternativa al precedente quest'altro script che costruisce manualmente il dataframe, copiatelo e incollatelo nella Console di R e premete ↵ Invio:

# GRAFICO PRIMA-DOPO (SLOPEGRAPH) dataframe costruito manualmente

#

library(CGPfunctions) # carica il pacchetto necessario

Caso <- c("paziente_A", "paziente_A", "paziente_B", "paziente_B", "paziente_C", "paziente_C", "paziente_D", "paziente_D", "paziente_E", "paziente_E") # vettore con i casi

Tempo <- c("t0", "t1", "t0", "t1", "t0", "t1", "t0", "t1", "t0", "t1") # vettore con i tempi

Valore <- c(1.5, 1.7, 1.7, 1.9, 2.1, 2.2, 1.6, 1.9, 2.4, 2.4) # vettore con i valori

mydata <- data.frame(Caso, Tempo, Valore) # genera il dataframe

mydata # mostra i dati

windows() # apre una nuova finestra

#

newggslopegraph(dataframe=mydata, Grouping=Caso, Times=Tempo, Measurement=Valore, Title="VEMS prima (t0) e dopo (t1) l'assunzione di un broncodilatatore", SubTitle="", Caption=NULL, XTextSize=15, YTextSize=4, DataTextSize=4, LineThickness=2) # traccia il grafico prima-dopo (slopegraph)

#

Dopo avere importato i dati dal file con il primo script, o avere costruito manualmente il dataframe con il secondo script, e avere mostrato i dati importati

> mydata # mostra i dati 

         Caso Tempo Valore

1  paziente_A    t0    1.5

2  paziente_A    t1    1.7

3  paziente_B    t0    1.7

4  paziente_B    t1    1.9

5  paziente_C    t0    2.1

6  paziente_C    t1    2.2

7  paziente_D    t0    1.6

8  paziente_D    t1    1.9

9  paziente_E    t0    2.4

10 paziente_E    t1    2.4

il grafico è realizzato in entrambi i casi in modo identico mediante la funzione newggslopegraph() nella quale gli argomenti sono:

→ dataframe, la tabella contenente le tre variabili da elaborare;

→ Grouping, la variabile che contiene l'identificativo univoco in base al quale aggregare i dati (il caso);

→ Times, la variabile che mette in sequenza, sull'asse orizzontale, i dati di uno stesso caso;

→ Measurement, la variabile che contiene i singoli dati di un caso, che vengono riportati in scala sull'asse verticale;

→ Title, il titolo del grafico, e un possibile sottotitolo (SubTitle), qui lasciato vuoto;

→ Caption, una possibile didascalia, qui non impiegata (= NULL);

→ XTextSize, la dimensione del testo per la variabile Times;

→ YTextSize, la dimensione del testo per la variabile Grouping;

→ DataTextSize, la dimensione del testo per la variabile Measurement;

→ LineThickness, lo spessore della linea che unisce i due valori di ciascun caso.

I colori delle linee sono generati automaticamente dalla funzione newggslopegraph().

Per ulteriori approfondimenti potete digitare help(newggslopegraph) nella Console di R ovvero, cosa che raccomando sempre, scaricate e consultate il manuale di riferimento del pacchetto [2].

Il grafico è piuttosto interessante in quanto evidenzia una tendenza all'aumento dei valori dopo l'assunzione del farmaco, tendenza che potrà poi essere ulteriormente valutata in termini di significatività statistica impiegando un test per il confronto tra dati appaiati [4].

----------

[1] Campbell MJ, Machin D. Medical Statistics. A Commonsense Approach. John Wiley & Sons, New York, 1993, ISBN 0-471-93764-9, p. 142.

[2] Vedere il manuale di riferimento del pacchetto: Package ‘CGPfunctions’.

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

[3] Vedere anche il post Importazione dei dati da un file .csv.

[4] Vedere il post Test parametrici e non parametrici per dati appaiati.

Grafici a scatola con i baffi (boxplot)

Il modo più semplice per generare i grafici a scatola con i baffi (box-and-whiskers plot, denominazione in genere abbreviata in boxplot) è impiegare la funzione boxplot() inclusa nel pacchetto base di R che utilizza i seguenti criteri di rappresentazione [1]:

→ le scatole (box) includono il 50% delle osservazioni;

→ il bordo inferiore delle scatole corrisponde al 25° percentile o primo quartile (Q₁);

→ la linea interna alle scatole corrisponde al 50° percentile o secondo quartile (Q₂) ovvero alla mediana;

→ il bordo superiore delle scatole corrisponde al 75° percentile o terzo quartile (Q₃);

→ i baffi (whiskers) corrispondono al valore minimo (baffo inferiore) e al valore massimo (baffo superiore) osservati dopo avere escluso gli outliers (vedi sotto);

→ la differenza interquartile [2] viene definita come IQR = Q₃ – Q₁ ovvero come differenza tra il valore corrispondente al terzo quartile (Q₃) e il valore corrispondente al primi quartile (Q₁);

→ i valori inferiori a Q₂ - 1.5 · IQR e i valori superiori a Q₂ + 1.5 · IQR sono considerati outliers (dati anomali o dati aberranti), sono esclusi dal computo dei valori minimo e massimo (vedi sopra), e sono riportati come punti singoli separati.

I boxplot forniscono una analisi non parametrica dei dati complementare a quella numerica. Qui li impieghiamo per effettuare l'analisi della concentrazione degli eritrociti (globuli rossi) nel sangue rilevata in 202 atleti australiani riportata nella colonna/variabile rcc della tabella ais inclusa nel pacchetto DAAG. Accertatevi di avere installato il pacchetto o in alternativa procedete come indicato in [3]. La conclusione? 

Con R è possibile realizzare grafici a scatola con i baffi con una sola e semplicissima riga di codice!

Iniziamo con queste quattro righe di codice, copiatele e incollatele nella Console di R e premete ↵ Invio:

# GRAFICI A SCATOLA CON I BAFFI suddivisi in base a un fattore

#

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

windows() # apre e inizializza una nuova finestra grafica 

#

boxplot(rcc~sport, data=ais, horizontal=FALSE, main="Eritrociti per sport praticato", xlab="Sport praticato", ylab="Eritrociti in 10^12/L", notch=FALSE, col="yellow") # eritrociti per ciascuno sport praticato

#

Le prime due righe si limitano a:

→ caricare con la funzione library() il pacchetto DDAG che include il set di dati ais;

→ aprire una finestra grafica con la funzione  windows().

Per tracciare il boxplot viene impiegata una sola riga di codice con la funzione boxplot() e con questi argomenti:

→ rcc~sport che indica che la rappresentazione degli eritrociti sotto forma di boxplot deve essere effettuata aggregando (~) i valori della variabile eritrociti (rcc) nei sottoinsiemi corrispondenti ai valori del fattore/variabile sport (B_Ball, Field, Gym, Netball, Row, Swim, T_400m, T_Sprnt, Tennis, W_Polo);

→ data=ais che specifica il nome della tabella che contiene la variabile rcc e la variabile sport - notare che questo argomento è superfluo se vi riferite alle variabili della tabella ais indicandole con il nome completo ais$nomedellavariabile, ovvero nel nostro caso ais$rcc, ais$sex, ais$sport, e così via;

→ horizontal=FALSE che indica che i boxplot devono essere orientati verticalmente;

→ main="..." che riporta il titolo del grafico;

→ xlab="..." che riporta l'etichetta da applicare all'asse x delle ascisse;

→ ylab="..." che riporta l'etichetta da applicare all'asse y delle ordinate:

→ notch=FALSE che esclude dai boxplot la rappresentazione dell'incisura che riporta i limiti di confidenza al 95% della mediana dei valori osservati;

→ col="..." che definisce il colore di riempimento dei boxplot.

Sull'asse orizzontale apparentemente non c'è spazio sufficiente per riportare le denominazioni di tutti gli sport: ma se "afferrate" il lato sinistro della finestra con il mouse ed estendete la finestra grafica in orizzontale vedete comparire le denominazioni mancanti.

In questo secondo script rielaboriamo i dati aggregando i valori per sesso con l'argomento rcc~sex:

# GRAFICI A SCATOLA CON I BAFFI con i limiti di confidenza della mediana

#

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

windows() # apre e inizializza una nuova finestra grafica 

#

boxplot(rcc~sex, data=ais, horizontal=FALSE, main="Eritrociti per sesso", xlab="Sesso", ylab="Eritrociti in 10^12/L", notch=TRUE, col="green") # eritrociti per sesso

#

Impiegando l'argomento notch=TRUE ora sui lati dei boxplot sono comparse le tacche (notch) o se preferite le incisure che rappresentano i limiti di confidenza al 95% della mediana.

Questo corrisponde ad un test per la significatività della differenza tra le mediane: se le tacche o incisure di due boxplot non si sovrappongono, pur sovrapponendosi in parte le distribuzioni dei dati, come in questo caso, la conclusione è che le mediane delle due distribuzioni differiscono significativamente.

Il fatto interessante è quindi che una rappresentazione grafica può essere impiegata non solo per effettuare una analisi esplorativa dei dati, ma addirittura per effettuare un test statistico (un confronto tra mediane) che ci conferma il fatto che la mediana della concentrazione degli eritrociti nelle atlete è inferiore a quella dagli atleti (è all'incirca di 4.4 · 10¹²/L contro 5.0 · 10¹²/L) e che la differenza tra le mediane, all'incirca di 0.6 · 10¹²/L, è statisticamente significativa come risulta dal confronto fra i due boxplot.

Da notare che quando nella funzione boxplot() si pone l'argomento notch=TRUE potrebbe comparire nella Console di R un messaggio che avverte che in alcuni casi le incisure sono uscite dai bordi della scatola e che suggerisce di valutare l'opportunità di sostituire l'argomento con notch=FALSE.

Questo vi accadrà quando eseguite questo terzo script, che applica l'argomento notch=TRUE ai boxplot differenziati per sport:

# GRAFICI A SCATOLA CON I BAFFI con i limiti di confidenza della mediana

#

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

windows() # apre e inizializza una nuova finestra grafica

#

boxplot(rcc~sport, data=ais,  horizontal=FALSE, cex.axis=0.8, las=2, main="Eritrociti per sport praticato", xlab="Sport praticato", ylab="Eritrociti in 10^12/L", notch=TRUE, col="green") # eritrociti per ciascuno sport praticato

#

In questo caso il grafico

evidenzia la comparsa del problema per gli sport Field, Gym, T_sprnt, Tennis e nella Console di R viene riportato questo messaggio:

> boxplot(rcc~sport, data=ais,  horizontal=FALSE, cex.axis=0.8, las=2, main="Eritrociti per sport praticato", xlab="Sport praticato", ylab="Eritrociti in 10^12/L", notch=TRUE, col="green") # eritrociti per ciascuno sport praticato

Messaggio di avvertimento:

In (function (z, notch = FALSE, width = NULL, varwidth = FALSE,  :

alcune tacche sono andate fuori dai cardini ('box'): potresti impostare notch=FALSE

> #

Il problema è causato dal fatto che per gli sport in questione il numero delle osservazioni è troppo ridotto. In questi casi vi sono solamente due modi per superare il problema:

→ mettere l'argomento notch=FALSE, come consigliato nella nota "potresti impostare notch=FALSE", e quindi rinunciare al "test grafico" di significatività;

→ aumentare adeguatamente, sempre che sia possibile farlo, il numero delle osservazioni.

Da notare che nella funzione boxplot() sono stati aggiunti due argomenti che consentono di ricavare lo spazio necessario per riportare sotto ai boxplot degli eritrociti le denominazioni degli sport senza dover ridimensionare la finestra grafica:

→ l'argomento cex.axis=0.8 che riduce lievemente la dimensione dei caratteri;

→ l'argomento las=2 che ruota verticalmente le etichette.

Ora copiate questo script, incollatelo nella Console di R e premete ↵ Invio:

# GRAFICI A SCATOLA CON I BAFFI con due grafici nella stessa immagine

#

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

windows() # apre e inizializza una nuova finestra grafica

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

#

boxplot(ais$pcBfat~ais$sport, horizontal=TRUE, boxwex=0.4, cex.axis=0.8, las=1, xlab="Grasso corporeo in %", ylab="Sport praticato", notch=FALSE, col="yellow") # valori della percentuale di grasso corporeo aggregati per sport

#

boxplot(ais$wt~ais$sport, horizontal=TRUE, boxwex=0.4, cex.axis=0.8, las=1, xlab="Peso corporeo in kg", ylab="Sport praticato", notch=FALSE, col="yellow") # [2] valori del peso corporeo aggregati per sport

#

In questo caso, per riportare due grafici affiancati:

→ con par(mfrow=c(1,2)) è stata predisposta la suddivisione della finestra in 1 riga e 2 colonne;

→ i boxplot sono riportati in orizzontale (horizontal=TRUE);

→ la larghezza dei box è stata ridotta (boxwex=0.4);

→ è impiegato l'argomento cex.axis=0.8 che riduce lievemente la dimensione dei caratteri;

→ le etichette delle ascisse sono riportate in orizzontale (las=1), 

→ non vengono riportate le incisure (notch=FALSE) che rappresentano i limiti di confidenza al 95% della mediana.

Inoltre qui, in assenza dell'argomento data=ais, nella funzione boxplot() sono stati  riportati per intero – cioè con il prefisso ais$ che indica la tabella che le contiene – i nomi della variabili.

Con quest'ultimo script riportiamo in un unico grafico i boxplot suddivisi in base a due fattori: sport e sesso. Questa volta la novità consiste in una legenda liberamente posizionabile: dopo avere tracciato i grafici separati per sport praticato e per sesso, lo script rimane in attesa. A questo punto bisogna posizionare il mouse dove si vuole che compaia la legenda [4] e fare click con il tasto sinistro del mouse per farla comparire e terminare lo script.

# GRAFICI A SCATOLA CON I BAFFI suddivisi in base a due fattori

#

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

windows() # apre e inizializza una nuova finestra grafica 

#

boxplot(rcc~sex+sport, data=ais, horizontal=FALSE, boxwex=0.4, cex.axis=0.8, las=2, main="Eritrociti per sport praticato e per sesso", xlab="", ylab="Eritrociti in 10^12/L", notch=FALSE, col=c("green", "yellow")) # eritrociti per ciascuno sport praticato

#

legend(locator(1), legend=c("sesso f","sesso m"), fill=c("green", "yellow")) # posiziona la legenda

#

Il codice non consente di spostare la legenda. Se non si è soddisfatti della sua posizione, è necessario rieseguire l’intero blocco di codice e fare nuovamente click con il tasto sinistro del mouse nel punto in cui si vuole posizionare la legenda.

Nella funzione boxplot() sono ora impiegati tre argomenti per risolvere il problema dello spazio necessario per riportare nella finestra grafica i boxplot degli eritrociti per tutti gli sport e per entrambi i sessi:

→ l'argomento boxwex=0.4 che riduce la larghezza dei boxplot;

→ l'argomento cex.axis=0.8 che riduce lievemente la dimensione dei caratteri;

→ l'argomento las=2 che ruota verticalmente le etichette.

La legenda viene realizzata mediante la funzione legend() impiegando come argomenti:

→ la funzione locator() che resta in attesa, legge la posizione del cursore grafico quando viene fatto click con il tasto sinistro del mouse, e posiziona la legenda;

→ legend= che specifica le due righe di testo che compaiono nella legenda;

→ fill= che determina la comparsa accanto alle righe di testo di due quadrati dei colori specificati.

Chiudo con tre suggerimenti che potrebbero essere utili:

→ trovate come combinare adeguatamente in un'unica immagine più grafici a scatola con i baffi nel post Inserire più grafici nella stessa immagine;

→ potete combinare i grafici a scatola con i baffi con il kernel density plot dei dati impiegando i grafici a violino (violin plot) [5];

→ quando i dati sono pochi, e preferite conservare il dettaglio dei singoli valori, potete impiegare in alternativa i grafici a punti (dotplot).

Infine se siete interessati ad una estetica alternativa potete impiegare le funzioni incluse nel pacchetto ggplot2 [6].

----------

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

[2] In italiano oltre che differenza interquartile si impiega anche scarto o ampiezza o range interquartile, in inglese si impiega Inter Quartile Range, da cui IQR.

[3] Vedere il post Il set di dati ais nel quale trovate anche come caricare i dati della tabella senza impiegare il pacchetto DAAG. 

[4] Più precisamente fare click con il tasto sinistro del mouse nel punto in cui si vuole posizionare l'angolo superiore sinistro della legenda.

[5] Richiede l'installazione del pacchetto aggiuntivo ggplot2.

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

[6] Vedere il post Grafici a scatola con i baffi (boxplot) con ggplot.