Archivio dell'autore: Davide Massidda

Davide Massidda

Informazioni su Davide Massidda

Dottore di ricerca in psicologia sperimentale, mi occupo di consulenza statistica nel campo delle scienze psicologiche, sanitarie e sociali.

Gestione avanzata dei fattori

Nell’ultimo post abbiamo parlato delle variabili di tipo fattore, per cercare di capire in cosa queste si differenzino dalle più semplici variabili di tipo carattere o numerico. Ora vedremo qualche funzione per manipolare i fattori, per ricodificarli e adattarli alle nostre esigenze.

Procuriamoci subito dei dati su cui lavorare. Sul portale open data del Comune di Cagliari è disponibile un dataset che contiene tutti gli incidenti stradali del 2014 avvenuti nel territorio comunale in cui sia intervenuta una pattuglia del corpo di Polizia Locale. Qui c’è la pagina di descrizione del dataset, mentre questo è il link diretto al file csv.

incidenti <- read.csv2("http://www.comune.cagliari.it/resources/cms/documents/Cagliari_PoliziaLocale_SinistriStradali_2014.csv")
&#91;/code&#93;

Il dataset contiene un bel po' di variabili (descritte sommariamente <a href="http://www.comune.cagliari.it/portale/it/opendata_visualizza_contenuto.page?modelId=45&contentId=SCH80418" target="_blank">qui</a>); noi ci concentreremo sul fattore <b>tipologia</b>, che classifica gli incidenti per tipo. Nel seguito di questo post eseguiremo una serie di operazioni su questo fattore; per evitare di alterare la variabile originale, salviamo la colonna &ldquo;tipologia&rdquo; del dataset &ldquo;incidenti&rdquo; in un nuovo oggetto, esterno al data frame:

[code language="R"]
tipoIncidenti <- incidenti$tipologia
&#91;/code&#93;

<h3 style="margin-top: 1.5em; margin-bottom: 0em;">I livelli del fattore</h3>
<hr style="margin-top: 0.4em; margin-bottom: 1.5em; width: 80%; background-color: black; height: 2px;">

<b>Prima domanda:</b> quanti livelli presenta il fattore tipologia? La risposta ce la fornisce il comando <b>nlevels</b>:

[code language="R"]
> nlevels(tipoIncidenti)
[1] 5

Seconda domanda: quali sono i livelli del fattore tipologia? La risposta ce la fornisce il comando levels, che estrae il vettore di etichette:

> levels(tipoIncidenti)
[1] "Danneggiamento"                         
[2] "Incidente stradale con feriti"          
[3] "Incidente stradale mortale"             
[4] "Incidente stradale per caduta di pedone"
[5] "Incidente stradale senza feriti"

Possiamo anche estrarre ogni singolo livello utilizzando le parentesi quadre. Per esempio, per estrarre il quarto livello:

> levels(tipoIncidenti)[4]
[1] "Incidente stradale per caduta di pedone"

Il fattore raggruppa gli incidenti sulla base di cinque categorie: la prima comprende generici danneggiamenti derivanti da buche sul manto stradale, caduta di rami dagli alberi, urto contro cassonetti, ecc., mentre le altre categorie comprendono incidenti stradali veri e propri.

Terza domanda: quante osservazioni sono presenti per ogni livello del fattore? La risposta ce la fornisce il comando table:

> table(tipoIncidenti)
tipoIncidenti
                         Danneggiamento           Incidente stradale con feriti 
                                    170                                     571 
             Incidente stradale mortale Incidente stradale per caduta di pedone 
                                      5                                      66 
        Incidente stradale senza feriti 
                                    634

Nel caso in cui le etichette che descrivono i nomi dei fattori siano lunghe (come in questo caso), possiamo associare il comando cbind al comando table, in modo da avere un output più ordinato:

> cbind(table(tipoIncidenti))
                                        [,1]
Danneggiamento                           170
Incidente stradale con feriti            571
Incidente stradale mortale                 5
Incidente stradale per caduta di pedone   66
Incidente stradale senza feriti          634

Ricodificare i livelli


Le etichette che costituiscono i nomi dei livelli sono piuttosto lunghe. Soprattutto quando si vogliono visualizzare i dati in un grafico, è bene che le etichette siano quanto più concise possibile. Possiamo rinominare i livelli sfruttando sempre il comando levels, utilizzandolo per attribuire al fattore un nuovo vettore di nomi:

newlev <- <- c("Danni", "con feriti", "con morti", "caduta pedone", "senza feriti")
levels(tipoIncidenti) <- newlev
&#91;/code&#93;

Attenzione: all'interno del vettore <i>newlev</i>, i nomi devono essere disposti nello stesso ordine in cui sono disposti i livelli nel fattore, pena l'attribuire ai livelli le etichette sbagliate.

Verifichiamo il risultato:

[code language="R"]
> levels(tipoIncidenti)
[1] "Danni"         "con feriti"    "con morti"     "caduta pedone"
[5] "senza feriti"

> table(tipoIncidenti)
tipoIncidenti
        Danni    con feriti     con morti caduta pedone  senza feriti 
          170           571             5            66           634

Possiamo ricodificare i livelli anche singolarmente. Per esempio, per trasformare l’etichetta “Danni” in “danni”:

> levels(tipoIncidenti)[1] <- "danni"
> levels(tipoIncidenti)
[1] "danni"         "con feriti"    "con morti"     "caduta pedone"
[5] "senza feriti"

Accorpare i livelli


Immaginiamo di voler raggruppare i sinistri in due categorie: quelli con feriti o morti (“tipo 1”) e quelli senza alcuna conseguenza fisica per le persone coinvolte (“tipo 2”). Per eseguire questa operazione possiamo sfruttare ancora una volta la funzione levels, attribuendo la medesima etichetta ai livelli che vogliamo accorpare:

levels(tipoIncidenti) <- c("tipo 2", "tipo 1", "tipo 1", "tipo 2", "tipo 2")
&#91;/code&#93;

Verifichiamo il risultato:

&#91;code language="R"&#93;
> levels(tipoIncidenti)
[1] "tipo 2" "tipo 1"

> table(tipoIncidenti)
tipoIncidenti
tipo 2 tipo 1 
   870    576

Come si vede, ci sono 576 incidenti di tipo 1 (somma di “con feriti” e “con morti”) e 870 incidenti di tipo 2 (somma di “danni”, “caduta pedone” e “senza feriti”).

Riordinare i livelli


La precedente aggregazione dei livelli ha prodotto un risultato per alcuni versi sconveniente: l’etichetta “tipo 2” costituisce il primo livello, mentre l’etichetta “tipo 1” il secondo livello. Ecco cosa succede se proviamo a visualizzare in un grafico a barre il numero di incidenti per ogni tipologia:

> barplot(table(tipoIncidenti))

barplot incidenti

I livelli vengono visualizzati nell’ordine con cui sono stati codificati, per cui il tipo 2 compare prima del tipo 1. Per ovviare a questo problema possiamo usare la funzione factor, che consente di rifattorizzare una variabile:

> tipoIncidenti <- factor(tipoIncidenti, levels=c("tipo 1","tipo 2"))
&#91;/code&#93;

Adesso i livelli sono stati riordinati, per cui sarà &ldquo;tipo 1&rdquo; a comparire per primo:

&#91;code language="R"&#93;
> table(tipoIncidenti)
tipoIncidenti
tipo 1 tipo 2 
   576    870
> barplot(table(tipoIncidenti))

barplot incidenti

Si noti che la stessa funzione factor, al pari della funzione levels, può essere usata per modificare i nomi dei livelli:

> tipoIncidenti <- factor(tipoIncidenti, labels=c("Tipo 1","Tipo 2"))
&#91;/code&#93;

<h3 style="margin-top: 1.5em; margin-bottom: 0em;">Fattori con livelli ordinati</h3>
<hr style="margin-top: 0.4em; margin-bottom: 1.5em; width: 80%; background-color: black; height: 2px;">

Per come l'abbiamo concepito fino a questo momento, un fattore non è altro che una variabile qualitativa su <b>scala nominale</b>. In realtà, il tipo fattore è pensato anche per registrare dati misurati su <b>scala ordinale</b>. Ripartiamo dall'inizio e riprendiamo l'originale vettore &ldquo;tipologia&rdquo; del dataset &ldquo;incidenti&rdquo;, salviamolo nell'oggetto tipoIncidenti abbreviando i nomi dei livelli:

[code language="R"]
tipoIncidenti <- factor(incidenti$tipologia, labels=c("danni", "con feriti", "con morti", "caduta pedone", "senza feriti"))
&#91;/code&#93;

Ammettiamo che sia possibile dare un ordinamento alle categorie del fattore sulla base della gravità: al livello più basso ci sarà la categoria "Danneggiamento", poi "Incidente stradale per caduta di pedone", poi "Incidente stradale senza feriti", quindi "Incidente stradale con feriti" e infine "Incidente stradale mortale".

Per prima cosa, creiamo un vettore che contiene i nomi dei livelli disposti nell'ordine desiderato:
&#91;code language="R"&#93;
lev <- c("danni", "caduta pedone", "senza feriti", "con feriti", "con morti")
&#91;/code&#93;

Si noti che avremmo potuto ottenere lo stesso risultato anche utilizzando il comando levels:

&#91;code language="R"&#93;
> lev <- levels(tipoIncidenti)&#91;c(1,4,5,2,3)&#93;
> lev
[1] "danni"         "caduta pedone" "senza feriti"  "con feriti"   
[5] "con morti"

A questo punto, rifattorizziamo il vettore utilizzando il comando factor, impostando a TRUE l’argomento ordered:


tipoIncidenti <- factor(tipoIncidenti, levels=lev, ordered=TRUE) [/code] Il comando str, adesso, indicherà che tipoIncidenti è una variabile di tipo fattore, ma i cui livelli sono ordinati, dove danni < caduta pedone < senza feriti < con feriti < con morti. [code language="R"] > str(tipoIncidenti)
Ord.factor w/ 5 levels “danni”<"caduta pedone"<..: 3 3 3 3 3 1 2 4 3 4 .. [/code] Si noti che lo stesso risultato avremmo potuto ottenerlo con la funzione ordered:


> tipoIncidenti <- ordered(tipoIncidenti, levels=lev) [/code] L'ordinamento dei livelli di un fattore è una sottigliezza che assume una sua utilità esclusivamente nell'adattamento di un modello statistico. È comunque importante che una variabile misurata su scala ordinale sia codificata nel modo corretto, soprattutto quando ci si passa i dati tra colleghi e certi dettagli rischiano di essere omessi e talvolta male interpretati.

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.

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.

Visualizzare relazioni curvilinee

Un modello statistico ha lo scopo di descrivere e prevedere uno o più fenomeni, semplificandoli, per renderli comprensibili all’uomo. Sotto questo punto di vista, i modelli statistici lineari costituiscono una soluzione ottimale per descrivere la realtà: sono semplici, facilmente interpretabili e adattabili a molte situazioni.

Nella sua formulazione più semplice, un modello lineare descrive la relazione tra due variabili attraverso una retta. Ci sono dei casi, però, in cui la relazione tra due variabili non può essere descritta attraverso una retta, ma è necessario interpolare i dati servendosi di una curva.

Descrivere una relazione di tipo curvilineo non implica necessariamente utilizzare un modello non lineare. La “linearità” di un modello, infatti, è determinata dalla forma matematica che questo assume e non dalla forma della relazione tra le variabili. Mi rendo conto che la cosa può apparire come un paradosso, ma è così: una relazione non lineare potrebbe benissimo essere descritta attraverso un modello lineare, purché questo modello sia lineare nei parametri.

In questo post vedremo un esempio di relazione non lineare e vedremo come visualizzare un modello di questo tipo utilizzando le funzioni basilari di R.

Un esempio: il compito di linea numerica


Una relazione non lineare è quella che si presenta quando un bambino piccolo, intorno ai 6/7 anni, è sottoposto al cosiddetto compito di linea numerica. Il compito consiste nel chiedere al bambino di posizionare una serie di numeri compresi entro un certo intervallo lungo un segmento. La versione più comune del compito richiede di posizionare lungo il segmento dei numeri che possono andare da 0 a 100.

Immaginiamo di presentare a un bambino trenta numeri compresi tra 0 e 100, che costituiscono gli stimoli, e di chiedergli di indicare in quale punto della linea dovrebbe essere posizionato ogni numero. Quindi, convertiamo ogni posizione indicata dal bambino nel corrispettivo numero che realmente si trova in quel punto della linea.

Registriamo nella variabile stimulus i numeri presentati e nella variable response le stime del bambino.

# Numeri presentati:
stimulus <- c(5,7,10,12,14,17,19,22,24,27,29,31,34,36,39,
            41,44,46,49,51,53,56,58,61,63,66,68,70,73,75)

# Risposte del bambino:
response <- c(38.22,42.67,40.63,45.81,44.09,44.15,45.7,52.6,
            49.51,48.13,49.68,51.83,49.4,53.48,50.06,48.76,
            53.22,51.54,53.04,56.08,53.24,50.9,52.12,53.11,
            55.95,54.81,53.96,54.62,50.92,55.61)
&#91;/code&#93;

Come ci ha insegnato Antonello Preti (<a href="">qui</a>), la primissima cosa che dovremmo fare è visualizzare la relazione tra le due variabili.

[code language="r"]
plot(stimulus, response)

Scatterplot tra stimolo e risposta nel compito di linea numerica

La relazione non sembra propriamente lineare, infatti i punti assumono un andamento leggermente curvilineo e tale curva sembra accentuata soprattutto per stimoli bassi piuttosto che elevati.

In ogni caso, se il bambino non ha risposto in maniera totalmente casuale, le sue risposte dovrebbero in qualche modo dipendere dagli stimoli presentati. Proviamo ad adattare un modello lineare:

model1 <- lm(response ~ stimulus)
summary(model1)
&#91;/code&#93;

Il modello presenta un indice R² pari a 0.6965. Non male. Vediamo allora come la retta di regressione interpola i punti:

&#91;code language="r"&#93;
plot(stimulus, response)
abline(model1, col="red")
&#91;/code&#93;

<img src="http://www.insular.it/wp-content/gallery/post-images/nlt-plot1.png" alt="Modello lineare risposta = a+b*stimolo" width="60%">

La retta non sembra interpolare benissimo i punti. Per carità, il modello non è malaccio, ma se diamo un'occhiata alla letteratura scientifica scopriremo che la relazione tra le due variabili non è lineare ma sembra essere <b>logaritmica</b>.

Quindi, a questo punto sorge un problema: come possiamo applicare un modello lineare se la relazione è non lineare (nel caso logaritmica)? Semplice: si trasforma <i>stimulus</i>, in modo da rendere lineare la sua relazione con <i>response</i>. L'assunzione è che la relazione tra il logaritmo naturale di <i>stimulus</i> e <i>response</i> sia lineare.

Adattare questo modello con R è molto semplice:

[code language="r"]
model2 <- lm(response ~ log(stimulus))
summary(model2)
&#91;/code&#93;

L'indice R² è decisamente superiore rispetto al precedente: 0.8275. La relazione potebbe davvero essere di tipo logaritmico. Come possiamo visualizzare questo modello? La funzione <b>abline</b> disegna rette, ma qui noi dobbiamo rappresentare una curva.

Per rappresentare un modello del genere si utilizza la funzione <b>lines</b>, in questo modo:

[code language="r"]
plot(stimulus, response)
lines(stimulus, predict(model2), col="blue")

Modello lineare risposta = a+b*log(stimolo)

La funzione predict estrapola i valori previsti dal modello. Alla funzione lines bisogna passare la variabile già disposta sull’asse orizzontale (stimulus) e i valori previsti dal modello, estrapolati proprio attraverso la funzione predict. Il risultato è la curva visualizzata in azzurro, che sembra interpolare i punti decisamente meglio rispetto alla retta del modello precedente.