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.

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.