Alcune cose da sapere sui fattori

Quando importiamo un file di dati in R, tutti i valori che il software identifica come di tipo carattere vengono automaticamente convertiti in tipo fattore. Alcune funzioni di lettura dati, come read.table o read.csv, possiedono un argomento chiamato stringsAsFactors, che, se impostato a FALSE, evita questa conversione.

Una domanda che mi viene posta abbastanza di frequente da chi non è solito usare i modelli statistici è: «che cosa sono i fattori?». La prima volta che ho dovuto a rispondere a tale quesito, devo ammetterlo, mi sono trovato un po’ in difficoltà. Sembra assurdo, ma per chi si occupa di statistica inferenziale (a qualsiasi livello) il concetto è così tanto dato per scontato da rendere complicato trovare le parole giuste per descriverlo a chi non ha mai fatto statistica o non ha mai avuto a che fare con il disegno di un esperimento.

La cosa più difficile da spiagare è quale sia la differenza tra variabili di tipo carattere e variabili di tipo fattore. Infatti, quando i dati vengono importati in R, i caratteri – e solo i caratteri – vengono convertiti in “factor”. Ma che cosa cambia? Non avremmo potuto lasciare tutto com’era e lavorare con le variabili “character” invece che “factor”? Tanto, sempre di stringhe si tratta!

Invece no, non si tratta sempre di stringhe, e adesso proverò a spiegarvi perché.

La differenza tra caratteri e fattori


Il fattore è un tipo di dato che permette di etichettare un gruppo di osservazioni utilizzando delle categorie predefinite. Immaginiamo di aver registrato in una tabella alcuni dati relativi a sette personaggi Disney residenti a Paperopoli e dintorni:

Nome Indirizzo Genere
Paperina Vico II dei ciclamini, 2 F
Paperoga Piazza C. Barks, 33 M
Nonna Papera Viale degli uliveti, s.n. F
Gastone Via della Dea Bendata, 1 M
Paperino Via della marineria, 8

M
Brigitta Via dei ciclamini, 21 F
Zio Paperone Colle del deposito, s.n. M

Questa tabella di dati contiene tre variabili: il nome del personaggio, l’indirizzo di residenza e il sesso. In tutti e tre i casi, le osservazioni sono codificate attraverso dei caratteri, ma solo una di queste variabili potrà essere utilizzata come fattore: il genere. Il nome e l’indirizzo, infatti, rappresentano dei dati univoci e caratteristici di ogni singolo personaggio. Se dovessimo raggruppare i personaggi, di certo non potremo usare l’indirizzo come criterio di raggruppamento, visto che ogni personaggio abita in un luogo diverso, né tantomeno il nome, univoco per definizione. Il genere, al contrario, è una variabile che può assumere solo due valori (F: femmina, M: maschio) e consente di suddividere i personaggi in due gruppi, uno di tre unità (le femmine) e uno di quattro unità (i maschi).

La prima differenza tra il tipo fattore e il tipo carattere consiste quindi nel fatto che il fattore serve per raggruppare le unità sottoposte a osservazione, mentre la semplice variabile di tipo carattere no. Il fattore assume un numero limitato di modalità, dette livelli, che specificano il gruppo di appartenenza di un’unità/individuo.

La seconda differenza tra il tipo fattore e il tipo carattere è che i dati che andranno a rappresentare una variabile fattore non necessariamente devono essere registrati utilizzando dei caratteri. Nell’esempio dei personaggi Disney, per codificare il genere ho usato delle etichette, ovvero F e M, ma queste etichette sono state scelte esclusivamente a mia personalissima discrezione. Nessuno mi avrebbe impedito di scrivere per esteso “femmina” o “maschio” o utilizzare addirittura dei numeri, magari 1 per indicare le femmine e 0 i maschi. Sarebbe invece stato più complicato indicare l’indirizzo attraverso dei numeri e sicuramente i caratteri erano d’obbligo per specificare il nome del personaggio.

La gestione dei fattori in R


Il tipo fattore assume un senso e gioca un ruolo fondamentale quando si deve adattare un modello statistico. Si immagini di voler eseguire un confronto tra gruppi di individui, suddivisi sulla base di una qualche caratteristica. R deve sapere ogni individuo a quale gruppo appartiene, per cui in fase di analisi dei dati sarà necessario specificare una variabile (fattore, ovviamente) che indica al software a quale raggruppamento appartiene ogni unità sottoposta a misurazione.

Da questo indirizzo è possibile scaricare in formato csv il dataset descritto nella precedente tabella. Per importare i dati in R dobbiamo prima di tutto scaricarli in locale e quindi leggerli usando la funzione read.csv2, con l’accorgimento di impostare l’argomento stringsAsFactors come FALSE; saremo noi, in seguito all’importazione, a decidere cosa convertire in fattore e cosa no.

disney <- read.csv2("disney.csv", stringsAsFactors=FALSE)
&#91;/code&#93;

Visualizziamo la struttura del dataset:

&#91;code language="R"&#93;
> str(disney)
'data.frame':	7 obs. of  4 variables:
 $ Nome     : chr  "Paperina" "Paperoga" "Nonna Papera" "Gastone" ...
 $ Indirizzo: chr  "Vico II dei ciclamini, 2" "Piazza C. Barks, 33" "Viale degli uliveti, s.n." "Via della Dea Bendata, 1" ...
 $ Genere   : chr  "F" "M" "F" "M" ...
 $ Sesso    : int  1 0 1 0 0 1 0

Rispetto alla tabella, in questo dataset è presente un’ulteriore variabile, Sesso, che codifica il genere utilizzando dei numeri: 1 per le femmine e 0 per i maschi.

La funzione per ricodificare una variabile come fattore si chiama as.factor. Iniziamo con il ricodificare la variabile di tipo carattere Genere:

disney$Genere <- as.factor(disney$Genere)
&#91;/code&#93;

Vediamo quindi l'esito della ricodifica:

&#91;code language="R"&#93;
> str(disney$Genere)
 Factor w/ 2 levels "F","M": 1 2 1 2 2 1 2

> disney$Genere
[1] F M F M M F M
Levels: F M

Il comando str, che descrive la struttura della variabile, ci informa che Genere adesso è un fattore a due livelli: “F” e “M”. Anche quando visualizziamo tutti i valori contenuti nel vettore, di seguito alla sequenza di valori (alla voce Levels), R stampa l’elenco delle etichette che costituiscono i livelli dei fattori.

Soffermiamoci un attimo sull’esito dell’applicazione del comando str sulla variabile fattore. Di seguito all’elenco dei livelli, “F” e”M”, compare una sequenza di numeri: 1 2 1 2 2 1 2. R codifica i livelli di un fattore numerandoli; il numero 1 viene assegnato al primo livello, il numero 2 al secondo livello, il numero 3 al terzo livello, e così via.

La scelta dei numeri è strettamente legata all’ordine alfabetico delle etichette utilizzate per rappresentare i livelli. In ordine alfabetico, l’etichetta “F” viene prima dell’etichetta “M”, per cui a “F” sarà assegnato il numero 1 e a “M” il numero 2. Per questo motivo, trasformare una variabile di tipo fattore in una variabile di tipo numerico è un’operazione concessa in R:

> as.numeric(disney$Genere)
[1] 1 2 1 2 2 1 2

Ovviamente, nell’operazione di conversione da tipo fattore a tipo numerico, le etichette “F” verranno rimpiazzate dal numero 1 (perché F è il primo livello) e le etichette “M” verranno rimpiazzate dal numero 2 (perché M è il secondo livello).

Proviamo ora a convertire in fattore la variabile numerica Sesso:

> disney$Sesso <- as.factor(disney$Sesso)
> disney$Sesso
[1] 1 0 1 0 0 1 0
Levels: 0 1

Nella variabile fattore Sesso, 0 (che codifica i maschi) è il primo livello, mentre 1 (che codifica le femmine) è il secondo livello. Questo accade semplicemente perché, nella sequenza numerica, 0 viene prima di 1. Ora, quindi, è la categoria “maschio” che costituisce il primo livello, mentre la categoria “femmina” il secondo.

Si noti che, nonostante i dati siano stati etichettati con dei numeri, sui fattori non è consentito applicare alcuna funzione matematica:

> mean(disney$Sesso)
[1] NA
Warning message:
In mean.default(disney$Sesso) :
  argument is not numeric or logical: returning NA

Riconvertire il tipo fattore in tipo numerico


La rappresentazione numerica che R utilizza per codificare i fattori rischia di creare molta confusione. La categoria che noi abbiamo etichettato con 0 viene codificata da R con il numero 1, mentre la categoria che noi abbiamo etichettato con 1 viene codificata da R con il numero 2. Se volessimo riconvertire il fattore Sesso in una variabile numerica, rischieremmo di combinare un vero disastro. Questo è infatti l’esito della conversione:

> as.numeric(disney$Sesso)
[1] 2 1 2 1 1 2 1

Nel riconvertire in formato numerico la variabile, R non utilizza le etichette che noi abbiamo attribuito ai livelli, ma la codifica numerica che lui ha dato ai livelli. Come fare allora?

La soluzione è semplice: per riconvertire una variabile fattore in tipo numerico, prima la si trasforma in carattere e solo succcessivamente le stringhe potranno essere convertite in numeri. Ovvero:

> disney$Sesso <- as.character(disney$Sesso)
> disney$Sesso
[1] "1" "0" "1" "0" "0" "1" "0"
> disney$Sesso <- as.numeric(disney$Sesso)
> disney$Sesso
[1] 1 0 1 0 0 1 0

Oppure, in un unico comando:


disney$Sesso <- as.numeric(as.character(disney$Sesso)) [/code] Per ora concludiamo qui: direi che di roba da digerire ce n'è a sufficienza!

Rimodellare i dati

Quando dobbiamo organizzare una tabella per raccogliere dati provenienti da disegni sperimentali, la prima domanda che dovremmo porci è come organizzarli. All’università mi hanno sempre ripetuto che a ogni riga deve corrispondere un “soggetto”, mentre a ogni colonna una variabile. Questa regola, ampiamente condivisa e sulla base della quale sono programmati i software statistici, in alcuni casi può diventare ambigua. Nello specifico, questo accade quando un soggetto (unità statistica) viene misurato più volte.

Per chiarire il concetto, passiamo subito a un caso pratico. Si tratta di un esempio inventato e abbastanza semplicistico, che nella sua semplicità ci aiuterà a definire meglio il problema.

In un laboratorio di psicofisiologia del sonno si stanno testando gli effetti di una nuova psicoterapia di tipo cognitivo-comportamentale e di un trattamento di tipo farmacologico per la cura dell’insonnia. Entrambi i trattamenti hanno la durata di due mesi.
Sono state selezionate 10 persone insonni con problemi nella fase di addormentamento: 5 vengono trattate con psicoterapia (PT) e 5 con terapia farmacologica (FT).
Il giorno prima dell’inizio del trattamento, a 30 e a 60 giorni viene registrato il numero di ore dormite da ogni individuo.
Per verificare la loro ipotesi, i ricercatori vorrebbero eseguire un’analisi statistica per verificare se le persone rispondono diversamente ai due trattamenti, valutando se questi portano a delle differenze nel numero di ore di sonno.

Dovendo organizzare i dati di questo disegno sperimentale, o si inseriscono i soggetti per riga oppure le variabili per colonna, ma non è possibile soddisfare i due criteri simultaneamente. Vediamo perché.

Nella tabella qui sotto, sono mostrati i dati raccolti nel cosiddetto formato wide.

sogg tratt tempo.0 tempo.30 tempo.60
1 PT 5.60 3.30 8.90
2 PT 4.80 4.70 7.90
3 PT 3.90 5.30 8.50
4 PT 5.80 4.50 8.20
5 PT 3.50 3.70 6.90
6 FT 3.90 4.70 7.00
7 FT 4.10 4.90 5.60
8 FT 5.00 3.80 5.30
9 FT 4.00 4.70 5.30
10 FT 4.00 4.80 6.30

Come possiamo osservare, è vero che a ogni riga corrisponde un soggetto, ma non è vero che a ogni colonna corrisponde una variabile. Infatti, il fattore “tempo” è stato suddiviso in tre colonne, una per ogni livello: 0 giorni (inizio), 30 giorni e 60 giorni. “tempo.0”, “tempo.30” e “tempo.60” non sono delle variabili, ma sono tutti livelli di un’unica variabile. Software come SPSS o STATISTICA richiedono i dati disposti proprio in questo modo.

Un altro modo in cui i dati possono essere organizzati è il formato long, che possiamo vedere qui di seguito.

sogg tempo tratt ore
1 0 PT 5.60
1 30 PT 3.30
1 60 PT 8.90
2 0 PT 4.80
2 30 PT 4.70
2 60 PT 7.90
3 0 PT 3.90
3 30 PT 5.30
3 60 PT 8.50
4 0 PT 5.80
4 30 PT 4.50
4 60 PT 8.20
5 0 PT 3.50
5 30 PT 3.70
5 60 PT 6.90
6 0 FT 3.90
6 30 FT 4.70
6 60 FT 7.00
7 0 FT 4.10
7 30 FT 4.90
7 60 FT 5.60
8 0 FT 5.00
8 30 FT 3.80
8 60 FT 5.30
9 0 FT 4.00
9 30 FT 4.70
9 60 FT 5.30
10 0 FT 4.00
10 30 FT 4.80
10 60 FT 6.30

In questo formato, a ogni colonna corrisponde una variabile, ma a ogni riga corrisponde un’osservazione e non un soggetto. Infatti, dato che i soggetti sono misurati tre volte, per ognuno di essi sono presenti tre righe.

Ognuno di questi formati ha dei pro e dei contro. Il formato wide, essendo più compatto, consente di visualizzare con più facilità la tabella e permette di utilizzare ognuno dei tre tempi come se fosse una variabile diversa (in alcuni casi può tornare utile). Tuttavia, trovo che il formato long sia migliore da un punto di vista concettuale, perché consente di individuare immediatamente quali sono le variabili di uno studio, distinguendo chiaramente la dipendente “ore di sonno” dalla indipendente “tempo” (questo aspetto non è da sottovalutare nell’insegnamento della statistica).

Eccetto casi particolari, le funzioni R richiedono che i dati siano disposti in formato long. Però, capita spesso che i dati siano stati registrati in formato wide. Talvolta ci si può trovare nella situazione di dover passare i propri dati in formato long a un utente di SPSS che, comprensibilmente, li gradirebbe in formato wide. È quindi importante capire come poter passare da un formato all’altro utilizzando R.

Da wide a long


Costruiamo in R la tabella di dati in formato wide:

sonno.wide <- data.frame(
    sogg = factor(c(1,2,3,4,5,6,7,8,9,10)),
    tratt = c("PT","PT","PT","PT","PT","FT","FT","FT","FT","FT"), 
    tempo.0 = c(5.6,4.8,3.9,5.8,3.5,3.9,4.1,5,4,4),
    tempo.30 = c(3.3,4.7,5.3,4.5,3.7,4.7,4.9,3.8,4.7,4.8),
    tempo.60 = c(8.9,7.9,8.5,8.2,6.9,7,5.6,5.3,5.3,6.3)
)
&#91;/code&#93;

Per rimodellare i dati, trasformandoli in formato long, possiamo utilizzare il comando <b>reshape</b>, in questo modo:

[code]
sonno.long <- reshape(sonno.wide, varying=3:5, v.names="ore",
                      timevar="tempo", direction="long")
&#91;/code&#93;

La funzione reshape ha tanti argomenti che possono essere manipolati, io qui ho sfruttato giusto i principali. L'argomento <b>varying</b> richiede di specificare quali sono le colonne che dovranno essere impilate, cioè quelle colonne che contengono i dati (le ore di sonno): tempo.0, tempo.30 e tempo.60. Io qui ho le ho indicate specificando il loro indice di colonna, ma avrei anche potuto passare, invece degli indici, un vettore contenente i nomi delle variabili. Come visto in precedenza, queste tre colonne sono in realtà un'unica variabile, che nel formato wide era stata suddivisa. Ora, nel formato long, a partire da queste tre colonne se ne dovrà ottenere una nuova che contiene i dati. Questa nuova colonna assumerà il nome specificato nell'argomento <b>v.names</b>. Inoltre, al dataset dovrà essere aggiunta una colonna che indica in quale tempo ogni dato è stato registrato; il nome che questa colonna dovrà assumere è stato indicato nell'argomento <b>timevar</b>. Infine, bisogna esplicitare che il dataset dovrà essere convertito in formato long e questo viene dichiarato attraverso l'argomento <b>direction</b>.

Possiamo vedere di seguito un'anteprima del risultato:

[code]
> sonno.long[1:10,]
     sogg tratt tempo ore id
1.1     1    PT     1 5.6  1
2.1     2    PT     1 4.8  2
3.1     3    PT     1 3.9  3
4.1     4    PT     1 5.8  4
5.1     5    PT     1 3.5  5
6.1     6    FT     1 3.9  6
7.1     7    FT     1 4.1  7
8.1     8    FT     1 5.0  8
9.1     9    FT     1 4.0  9
10.1   10    FT     1 4.0 10

Da long a wide


Costruiamo in R la tabella di dati in formato long:

sonno.long <- data.frame(
    sogg = factor(c(1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,6,6,6,7,7,7,8,8,8,9,9,9,10,10,10)),
    tempo = factor(c(0,30,60,0,30,60,0,30,60,0,30,60,0,30,60,0,30,60,0,30,60,0,30,60,0,30,60,0,30,60)),
    tratt = c("PT","PT","PT","PT","PT","PT","PT","PT","PT","PT","PT","PT","PT","PT","PT",
        "FT","FT","FT","FT","FT","FT","FT","FT","FT","FT","FT","FT","FT","FT","FT"),
    ore = c(5.6,3.3,8.9,4.8,4.7,7.9,3.9,5.3,8.5,5.8,4.5,8.2,3.5,3.7,6.9,
        3.9,4.7,7,4.1,4.9,5.6,5,3.8,5.3,4,4.7,5.3,4,4.8,6.3)
)
&#91;/code&#93;

Usiamo ora la funzione reshape:

&#91;code&#93;
sonno.wide <- reshape(sonno.long, v.names="ore", timevar="tempo",
                      idvar="sogg", direction="wide")
&#91;/code&#93;

Stavolta, con <b>v.names</b> ho specificato al comando reshape qual è la variabile da suddividere nelle colonne, e, attraverso <b>timevar</b>, ho indicato sulla base di quale variabile bisogna creare queste colonne; inoltre, ho dovuto indicare a <b>idvar</b> qual è la variabile che codifica per i soggetti. Infine, non dobbiamo dimenticarci di indicare a <b>direction</b> che dobbiamo trasformare i dati in formato wide.

Visualizziamo il risultato:

[code]
> sonno.wide
   sogg tratt ore.0 ore.30 ore.60
1     1    PT   5.6    3.3    8.9
4     2    PT   4.8    4.7    7.9
7     3    PT   3.9    5.3    8.5
10    4    PT   5.8    4.5    8.2
13    5    PT   3.5    3.7    6.9
16    6    FT   3.9    4.7    7.0
19    7    FT   4.1    4.9    5.6
22    8    FT   5.0    3.8    5.3
25    9    FT   4.0    4.7    5.3
28   10    FT   4.0    4.8    6.3

Suddivisione dei dati


Quando i dati sono in formato long, talvolta può essere comodo spezzare il dataset in più oggetti a seconda del livello di una variabile. Per esempio, potremmo avere l’esigenza di suddividere i dati delle persone trattate con psicoterapia da quelli delle persone trattate con terapia farmacologica, magari per creare dei grafici separati per ogni gruppo. Questa operazione può essere realizzata molto facilmente attraverso il comando split:

sonno.split <- split(sonno.long, f=sonno.long$tratt)
&#91;/code&#93;

Il comando richiede, oltre al dataset, la variabile da utilizzare come criterio per la suddivisione (che nell'esempio è la colonna <i>tratt</i>). Come risultato, avremo un oggetto di tipo lista con due slot, uno per ogni livello della variabile <i>tratt</i>.

Per estrarre i dati di ogni slot, possiamo usare il dollaro, per esempio:

[code]
> sonno.split$FT
   sogg tempo tratt ore
16    6     0    FT 3.9
17    6    30    FT 4.7
18    6    60    FT 7.0
19    7     0    FT 4.1
20    7    30    FT 4.9
21    7    60    FT 5.6
22    8     0    FT 5.0
23    8    30    FT 3.8
24    8    60    FT 5.3
25    9     0    FT 4.0
26    9    30    FT 4.7
27    9    60    FT 5.3
28   10     0    FT 4.0
29   10    30    FT 4.8
30   10    60    FT 6.3

Per chi volesse divertirsi, ci sono due pacchetti, reshape e reshape2, che permettono di andare ben oltre quanto presentato in questo post. Qui l’articolo di presentazione dell’autore dei pacchetti.

Save the data!

In questo blog ci occupiamo di dati. Manipolazione, elaborazione, visualizzazione… cerchiamo di affrontare l’argomento sotto diversi punti di vista. Nella pratica, però, non tutto è così semplice come ci viene presentato in un post.

Talvolta, davanti al problema di analizzare dei dati, capita che i dubbi ci assalgano e non sappiamo bene come muoverci. Magari, dopo un primo approccio ai numeri in cui ci sembrava di avere le idee ben chiare, ci perdiamo e non sappiamo più che strada seguire.

Per questo motivo, insieme all’associazione Sardinia Open Data, abbiamo organizzato l’evento Save the data!

Di cosa si tratta?

Save the data! sarà una serata informale dedicata all’analisi collaborativa dei propri dati. Ognuno di noi porterà un computer portatile e i propri dati, con lo scopo di collaborare con gli altri per ricevere consigli e a propria volta mettersi a disposizione di chi ha bisogno di una mano.

Come scrive l’associazione Sardinia Open Data sul suo blog:

« Davanti a un problema, la possibilità di poter ascoltare altri punti di vista e di poter contare su competenze diverse dalle proprie può fornire quel valore aggiunto in grado di metterci sulla giusta strada. Con i dati si possono raccontare storie, ma raramente nei dati è scritta un’unica storia e ancor più raramente esiste un solo modo di raccontarla. Elaborare dati significa fare delle scelte, e le scelte, quando sono frutto di collaborazione, portano a soluzioni più forti. »

Sul blog dell’associazione troverai informazioni più dettagliate su come iscriverti all’evento, che si terrà il 7 marzo a Cagliari, in via Tempio 22, a partire dalle 15:00.

Risolvere pRoblemi

Risolvere problemi è il metodo migliore per imparare qualcosa. Ricordate la scuola?

Con la programmazione in R è un po’ la stessa cosa. Risolvere problemi favorisce la dimestichezza con il programma e facilita l’apprendimento nell’uso dei vari pacchetti. Il suggerimento è di iniziare da cose semplici, come quelle che potete trovare in un manuale o nei siti di discussione.

Di recente l’utente Umesh Acharya ha postato un quesito su CrossValidated, un sito che si occupa di problemi di statistica. Il quesito in realtà è di natura informatica, poiché riguarda la codifica di un plot. Fino a poco tempo fa, il quesito era reperibile qui, ma è stato ormai rimosso, essendo quel genere di domande più adatte a StackOverflow, il sito gemello di CrossValidated, ma più genericamente dedicato a problemi informatici.

Domanda Umesh: “Come posso mettere tre plot in un’unica finestra?”. La domanda non è banale. Umesh ha creato una funzione che consente di “stampare” nel plot la formula di una regressione usando i valori stimati per i parametri, e per farlo ha usato ggplot2, il programma di grafica creato da Hadley Wickham, mago della programmazione in R (su InsulaR se n’era già parlato qui).

A differenza del sistema grafico base, in ggplot2 non è possibile ripartire le finestre grafiche con il comando par (che avete visto usare in precedenti post). Bisogna agire sulle griglie. E bisogna saperlo.

Dunque, abbiamo un dataset con 6 variabili, che vogliamo rappresentare a coppie di due, coordinate in una regressione. Come si rappresentano i tre grafici in un’unica finestra? Iniziamo a costruire il nostro dataset, che chiameremo ‘tom’ (toy object model):

a1 <- seq(from=1, to=10, by=1)
a2 <- seq (from=2, to=20, by=2)
b1 <- seq(from=3, to=30, by=3)
b2 <- seq(from=4, to=40, by=4)
c1 <- seq(from=5, to=50, by=5)
c2 <- seq(from=2, to=29, by=3)

# combiniamo i vettori in una matrice con il comando cbind
dataframe <- cbind(a1,a2,b1,b2,c1,c2)

# trasformiamo la matrice in un dataframe con etichette per l’intestazione delle variabili
tom <- as.data.frame(dataframe)

E adesso la bella funzione creata dal nostro Umesh per inserire il risultato della regressione nel grafico:

## function to create equation expression
lm_eqn = function(x, y, df){
    m <- lm(y ~ x, df);
    eq <- substitute(italic(y) == b %.% italic(x) + a,
    list(a = format(coef(m)[1], digits = 2),
    b = format(coef(m)[2], digits = 2)))
    as.character(as.expression(eq));
}

Chiamiamo il pacchetto:

library (ggplot2)

Creiamo i tre grafici in ggplot2:

################ a1 and a2 couple of variables
P1<-ggplot(tom, aes(x=a1, y=a2)) +
geom_point(shape=1) +      # Use hollow circles

geom_smooth(method=lm,    # Add linear regression line

se=FALSE)     # Don't add shaded confidence region

P2<- P1 + labs (x= "a1 value", y = "a2 value")## add x and y labels

P3 <- P2 + theme(axis.title.x = element_text(face="bold", size=20)) +
labs(x="a1 value")                            #label and change font size

P4 <- P3 + scale_x_continuous("a1 value",
limits=c(-10,60),
breaks=seq(-10, 60, 10))  ##adjust axis limits and breaks

##add regression equation using annotate
P5a <- P4 + annotate("text", x = 30, y = 10, label = lm_eqn(tom$a1, tom$a2, tom), color="black", size = 5, parse=TRUE)

P5a

Variabili a1 e a2

################ b1 and b2 couple of variables
P1b<-ggplot(tom, aes(x=b1, y=b2)) +
geom_point(shape=1) +      # Use hollow circles
geom_smooth(method=lm,    # Add linear regression line
se=FALSE)     # Don't add shaded confidence region

P2b<- P1b + labs (x= "b1 value", y = "b2 value")## add x and y labels

P3b <- P2b + theme(axis.title.x = element_text(face="bold", size=20)) +
labs(x="b1 value")                            #label and change font size
P4b <- P3b + scale_x_continuous("b1 value",
limits=c(-10,60),
breaks=seq(-10, 60, 10))  ##adjust axis limits and breaks

##add regression equation using annotate
P5b <- P4b + annotate("text", x = 40, y = 20, label = lm_eqn(tom$b1, tom$b2, tom), color="black", size = 5, parse=TRUE)

P5b

Variabili b1 e b2

################ c1 and c2 couple of variables
P1c<-ggplot(tom, aes(x=c1, y=c2)) +
geom_point(shape=1) +      # Use hollow circles
geom_smooth(method=lm,    # Add linear regression line
se=FALSE)     # Don't add shaded confidence region

P2c<- P1c + labs (x= "c1 value", y = "c2 value")## add x and y labels

P3c <- P2c + theme(axis.title.x = element_text(face="bold", size=20)) +
labs(x="c1 value")                            #label and change font size

P4c <- P3c + scale_x_continuous("c1 value",
limits=c(-10,60),
breaks=seq(-10, 60, 10))  ##adjust axis limits and breaks

##add regression equation using annotate
P5c <- P4c + annotate("text", x = 40, y = 15, label = lm_eqn(tom$c1, tom$c2, tom), color="black", size = 5, parse=TRUE)

P5c

Variabili c1 e c2

E adesso il gran finale. Come rappresentare i tre grafici in un’unica finestra?

Se il nostro Umesh, peraltro bravo programmatore, fosse stato meno pigro, avrebbe scoperto subito che CrossValidated non era il sito che faceva per lui, e avrebbe scoperto in fretta almeno due soluzioni al suo problema.

Le soluzioni che qui propongo sono tratte dal sito Cookbook for R, curato da Wiston Chang, una autentica miniera di gemme.

La prima soluzione utilizza il sistema di gestione delle griglie, e fa uso di due librerie aggiuntive a ggplot2 (librerie che dovrete scaricare se volete usarle).

library(grid)
library(gridExtra)

# semplice ed elegante (il titolo lo scegliete voi)
grid.arrange(P5a, P5b, P5c, ncol = 3, main = "Your main title")
# after saving, dev.off()
# se volete cambiare le dimensioni del titolo, seguite le istruzioni

### changing dimensions of the title
grid.arrange(P5a, P5b, P5c, ncol = 3,
main=textGrob("Your main title", gp=gpar(fontsize=20,font=3)))

Combinazione grafici

La seconda soluzione fa ricorso ad una funzione creata da Chang, e messa cortesemente a disposizione degli utenti dal creatore (il risultato è analogo a quello riportato sopra).


# Multiple plot function
# ggplot objects can be passed in …, or to plotlist (as a list of ggplot objects)
# – cols: Number of columns in layout
# – layout: A matrix specifying the layout. If present, ‘cols’ is ignored.
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
# from: http://www.cookbook-r.com/

multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) { require(grid) # Make a list from the ... arguments and plotlist plots <- c(list(...), plotlist) numPlots = length(plots) # If layout is NULL, then use 'cols' to determine layout if (is.null(layout)) { # Make the panel # ncol: Number of columns of plots # nrow: Number of rows needed, calculated from # of cols layout <- matrix(seq(1, cols * ceiling(numPlots/cols)), ncol = cols, nrow = ceiling(numPlots/cols)) } if (numPlots==1) { print(plots[[1]]) } else { # Set up the page grid.newpage() pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout)))) # Make each plot, in the correct location for (i in 1:numPlots) { # Get the i,j matrix positions of the regions that contain this subplot matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE)) print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row, layout.pos.col = matchidx$col)) } } } # Una volta che avete caricata la funzione nell’area di lavoro, l’uso è semplicissimo # (Il numero di colonne (in questo caso 3, ma può variare a piacimento): multiplot(P5a, P5b, P5c, cols=3) [/code] That's all, folks! Antonello Preti Stay Tuned for our next episode!

Creare una mappa colorata dell’Italia con mapIT

Attenzione: questo post è stato aggiornato rispetto alla pubblicazione originale per adeguare il codice agli ultimi aggiornamenti apportati alle funzioni descritte. Si ringrazia Marco Bertoletti per la segnalazione e il supporto nello sviluppo del pacchetto mapIT. [4 giugno 2015]
In R ci sono diverse pacchetti che permettono di disegnare mappe. Ormai ho perso il conto, ma quelle che sicuramente meritano almeno una citazione sono sp e ggmap. Nonostante la grossa disponibilità di funzioni R dedicate a questo tema, devo ammettere che, quando in passato mi sono trovato davanti al problema di rappresentare una semplicissima mappa dell’Italia con le regioni colorate a seconda dell’intensità di una variabile, sono entrato in crisi.

La maggior parte delle funzioni che avevo trovato erano molto semplici da usare per rappresentare aree degli Stati Uniti, ma dovendo raffigurare l’Italia le cose si complicavano parecchio. Mi aspettavo che, con R, realizzare una choropleth map (così si chiama il grafico in cui un’area geografica è rappresentata colorando in maniera diversificata le varie porzioni di territorio) fosse cosa assai banale. Non fu affatto così… trovai le procedure descritte nei tutorial fin troppo complicate rispetto alla banalità del grafico che volevo realizzare (qui si trova comunque una bella guida).

Mi resi conto di non essere l’unico utente R ad avere questo problema. Diverso tempo fa, su Statistica@Ning, Lorenzo di Blasio aveva proposto una soluzione in un bel tutorial che descriveva come costruire una mappa delle regioni italiane con ggplot2. Sintetizzando il codice proposto da Lorenzo, io avevo assemblato una funzione per creare la mappa in maniera rapida e semplificata; infine, Nicola Sturaro del gruppo MilanoR aveva preso in mano il codice e l’aveva decisamente migliorato e completato, inserendolo dentro un pacchetto: mapIT.

Attualmente, il pacchetto mapIT risiede in un repository su GitHub; per installarla si utilizza devtools, di cui abbiamo già parlato:

library(devtools)
install_github("quantide/mapIT")

La prima volta che utilizzai mapIT fu per un lavoro realizzato in collaborazione con Claudia Foti, studentessa all’Università di Messina, che aveva raccolto dei dati relativi alla valutazione di vini da parte di alcune guide specializzate. Avevamo la necessità di visualizzare, per ogni regione, il numero di cantine i cui vini erano stati recensiti.

Di seguito sono riportati i dati, organizzati in un data frame a due colonne, la prima che indica la regione e la seconda il numero di cantine.

vino <- data.frame(
    Regione = c("Abruzzo", "Basilicata", "Calabria", "Campania", "Emilia-Romagna",
                "Friuli-Venezia Giulia", "Lazio", "Liguria", "Lombardia", "Marche",
                "Molise", "Piemonte", "Puglia", "Sardegna", "Sicilia", "Toscana",
                "Trentino-Alto Adige", "Umbria", "Valle d\'Aosta", "Veneto"),
    Num.Cantine = c(22, 8, 9, 35, 24, 74, 19, 8, 41, 29, 5, 191, 22, 14, 40, 173,
                    57, 29, 6, 92)
)
&#91;/code&#93;

I nomi delle regioni possono essere scritti sia in minuscolo che in maiuscolo; spazi e altri caratteri non alfabetici vengono ignorati. Si può quindi scrivere indifferentemente: &lsquo;Trentino-Alto Adige&rsquo;, &lsquo;Trentino Alto Adige&rsquo; o &lsquo;TrentinoAltoAdige&rsquo;. Per le regioni con denominazione bilingue, viene riconosciuta la sola dicitura in italiano.

Per costruire la mappa delle regioni italiane, il pacchetto mapIT mette a disposizione l'omonima funzione <b>mapIT()</b>. Il primo argomento da passare alla funzione è la variabile numerica (Num.Cantine) e il secondo la variabile che specifica a quale regione italiana ogni dato si riferisce (Regione). Eventualmente, può essere passato un terzo argomento, costituito dal dataset dal quale le variabili andranno estrapolate.

Inoltre, sono disponibili dei parametri aggiuntivi per modificare l'aspetto grafico, che vanno passati attraverso una lista di nome <b>graphPar</b>. Io qui ne utilizzo uno, <b>guide.label</b>, che specifica l'etichetta che dovrà essere utilizzata come titolo della legenda a lato.

[code language="r"]
library(mapIT)
mapIT(Num.Cantine, Regione, data=vino, graphPar = list(guide.label="Numero\nCantine"))

La sequenza di escape ‘\n’, utilizzata nella stringa passata a guide.label, serve per mandare accapo il testo. Questo è il grafico risultante:

Mappa regioni italiane

Facile, no? È bastato caricare il pacchetto e lanciare una brevissima istruzione!

Il grafico può essere personalizzato agendo sull’argomento graphPar. Questa lista può contenere una lunga sequenza di dettagli che definiscono l’aspetto finale del grafico (per i dettagli, si veda l’help della funzione). Una delle prime cose che sicuramente vorremo fare sarà alterare i colori. Per alterare i colori bisogna specificare, nella lista da passare a graphPar, il colore da attribuire al valore minimo (low) e il colore da attribuire al valore massimo (high):

gp <- list(guide.label="Numero\nCantine", low="#fff0f0", high="red3")
&#91;/code&#93;

Per comodità, ho salvato questi valori nell'oggetto gp; si noti che i colori possono essere specificati sia utilizzando il <a href="https://it.wikipedia.org/wiki/Lista_dei_colori" target="_blank">codice esadecimale</a> ("#fff0f0") che utilizzando le parole chiave che definiscono i colori ("red3").

Ora, rilanciamo il comando passando la lista gp all'argomento graphPar:

[code language="r"]
mapIT(Num.Cantine, Regione, data=vino, graphPar=gp)

Mappa regioni italiane in rosso

Potete giocare con i colori per trovare la configurazione che preferite maggiormente. Per individuare il codice esadecimale dei colori si può usare un applicativo web come RGB color picker.

Per rendere il grafico in bianco e nero bisogna sempre giocare con i valori low e high di graphPar. In questo caso, per rendere il grafico un pochino più accattivante, potremmo sfruttare i temi di ggplot2; il primo esempio proposto di seguito sfrutta il tema theme_bw e produce il grafico di sinistra, mentre il secondo esempio sfrutta il tema theme_grey e produce il grafico di destra.


library(ggplot2)

# Tema: black and white
gp <- list(guide.label="Numero\nCantine", low="white", high="gray20", theme=theme_bw()) mapIT(Num.Cantine, Regione, data=vino, graphPar=gp) # Tema: grey gp <- list(guide.label="Numero\nCantine", low="white", high="gray20", theme=theme_grey()) mapIT(Num.Cantine, Regione, data=vino, graphPar=gp) [/code] Mappa regioni italiane in bianco e nero

Ci sono ancora diverse funzionalità da implementare e non escluderei che in futuro alcune cose possano cambiare. Chi avesse idee su come migliorare mapIT o dovesse riscontrare qualche malfunzionamento, può aprire un issue su GitHub in modo da darci una mano nello sviluppo.