mercoledì 30 dicembre 2009

Mappa dei Social Networks per nazione in R

Come ultimo post dell'anno 2009, ho deciso di presentare alcune funzioni di R che finora non avevo mai trattato, e per fare pratica con i dataframe in R. L'argomento che ho scelto come esempio da analizzare è la popolarità dei Social Network, idea che mi è venuta leggendo questo articolo. Non voglio certo soffermarmi in pensieri e riflessioni su questo genere di siti, essendo l'obiettivo di questo blog l'analisi statistica.
Osservando tutti i grafici presenti nel sito che ho linkato, una domanda viene spontanea: ma in definitiva, qual è il social network più popolare per ciascuna nazione? La risposta a questa domanda è l'argomento di questo post.

Innanzitutto è molto interessante lo strumento messo a disposizione da Google: Google Insights for Search. Pur essendo ancora in versione beta, tramite la semplice interfaccia che da sempre caratterizza il marchio Google, possiamo ottenere molte informazioni riguardanti essenzialmente le query più cercate nel celebre motore di ricerca. Questi dati vengono presentati per nazione, ed ecco perchè questo strumento ci sarà molto utile.

Per prima cosa cerchiamo di ottenere i dati che ci servono: nella spazio di ricerca, scriviamo i nomi dei 12 principali social network che vogliamo confrontare: MySpace, Facebook, Hi5, Friendster, LinkedIn, Orkut, Last.fm, LiveJournal, Xanga, Bebo, Imeem e Twitter. Non dobbiamo far altro che clickare su Search, e quindi sul link Download as a CSV, per ottenere la tabella di dati che ci serve.
In che formato sono i dati relativi alle nazioni? I dati sono stardardizzati: fissato come valore 100 la nazione con più query, a tutte le altre viene assegnato un valore in proporzione (stesso metodo di standardizzazione visto in questo post). Ottenuto il file .CSV, possiamo caricarlo in R, utilizzando questo comando:


social <- read.csv("C:/~/Desktop/report-5.csv", header=T, skip=321, nrows=467-322, sep=";")


Studiamo la sintassi:
- read.csv: funzione adatta proprio per i file in formato .csv di Excel
- "C:/...": path dove è posizionato il file; ricordo che in R vanno usati i backslash (/) al posto degli slash (\) per individuare il percorso di cartelle
- header=T: viene specificato nel file il nome delle colonne
- skip=321: R non deve leggere le prime 321 righe, ma partire dalla 322esima (perchè da questa riga in poi sono posizionati i dati che ci servono
- nrows=467-322: il numero di righe che ci serve arriva fino a 467. Specifichiamo quindi il numero di righe da analizzare a partire dalla 322esima, con la semplice differenza
- sep=";": il simbolo usato per separare, in ogni riga, i valori delle colonne

Ecco come appare il file appena caricato (clickate su View Plain e ingrandite la finestra, per motivi di spazio non entra tutto qui):


social[1:6,]
Region facebook orkut twitter imeem bebo friendster hi5 myspace last.fm linkedin livejournal xanga
1 Turkey 100 0 14 1 0 0 4 3 6 3 0 0
2 Tunisia 97 0 20 1 0 0 4 3 0 16 0 0
3 Croatia 77 0 17 0 0 0 0 10 59 15 2 0
4 Italy 75 0 18 0 0 0 1 13 20 33 3 0
5 Colombia 65 0 15 0 0 0 19 13 14 2 1 0
6 Venezuela 64 0 35 0 0 0 10 6 9 2 1 0
...


A questo punto, caricato il dataframe di dati, dobbiamo individuare per ciascuna nazione (ciascuna riga) il social network principale. Lavorare con i dataframe in R è davvero difficoltoso; la mia soluzione a questo problema è la seguente: per ogni riga andiamo a cercare il valore massimo, e scriviamo tale valore in un'altra colonna aggiunta al dataframe. Quindi confrontiamo la colonna dei massimi con le colonne dei singoli Social Network: quando i due valori coincidono, facciamo scrivere in una seconda colonna aggiunta il nome del Social Network.
Sono sicuro che esistano soluzioni più lineari ed eleganti della mia, che tuttavia funziona a dovere. Per fare questo, dobbiamo utilizzare due funzioni finora mai viste: il ciclo for..to..do, e lo statement if..then..else.

Punto primo: creiamo la colonna max con i massimi valori di ogni riga:


for(i in 1:145){
social$max[i] <- max(social[i, 2:13])}


Studiamo la sintassi:
- for(i in x:y){ codice }: questo è il codice generale del ciclo for..to..do in R. Tra parentesi tonde si sceglie un contatore (i) che assumerà i valori compresi tra i valori specificati x e y. Ad ogni ciclo, i aumenta di una unità, e vengono eseguite le operazioni scritte tra le parentesi graffe. Quando i raggiungerà il valore y (quindi dopo n = y - x cicli), l'operazione si conclude.
- max: è la funzione che sceglie il valore massimo tra i valori scritti tra parentesi. Questi valori sono, per ogni ciclo, le colonne dalla e alla 13 della riga i (dalla riga 1 alla riga 145)
- social$max[i]: questa è la colonna aggiuntiva al nostro dataframe. Per ogni ciclo viene calcolato il valore massimo della riga, e scritto nella posizione i-esima della colonna max del dataframe.

A questo punto il dataframe è il seguente:


social[1:6,]
Region facebook orkut twitter imeem bebo friendster hi5 myspace last.fm linkedin livejournal xanga max
1 Turkey 100 0 14 1 0 0 4 3 6 3 0 0 100
2 Tunisia 97 0 20 1 0 0 4 3 0 16 0 0 97
3 Croatia 77 0 17 0 0 0 0 10 59 15 2 0 77
4 Italy 75 0 18 0 0 0 1 13 20 33 3 0 75
5 Colombia 65 0 15 0 0 0 19 13 14 2 1 0 65
6 Venezuela 64 0 35 0 0 0 10 6 9 2 1 0 64
...


Come vedete è comparsa l'ultima colonna max che riporta il valore massimo di ogni riga.

Punto secondo: adesso dobbiamo dare un nome a questo valore massimo, e queste deve essere scelto da uno dei nomi delle altre colonne (ossia uno dei Social Networks).
Ecco il codice:


for(i in 1:13){
for(k in 1:145){
if(social$max[k]==social[k,i]){social$MaxS[k]=names(social)[i]}
}
}


Il codice è reso un pò più complesso dal fatto che sono concatenati due cicli for e una scelta if. Analizziamo prima il codice if, che in generale si presenta così: if(condizione){ codice }.
Il contatore k è relativo alle righe, mentre il contatore i è relativo alle colonne. La condizione if dice la seguente cosa: considerata una riga k (da 1 a 145) e in particolare il valore della colonna max, quando questo valore è uguale a uno dei valori delle i colonne (confronto con la colonna dalla 2 alle 13), allora alla colonna aggiuntiva MaxS viene dato il nome che è uguale a quello della colonna individuata.
E' difficile da spiegare a parole, ma il risultato è il seguente:


social[1:10]
Region facebook orkut twitter imeem bebo friendster hi5 myspace last.fm linkedin livejournal xanga max MaxS
1 Turkey 100 0 14 1 0 0 4 3 6 3 0 0 100 facebook
2 Tunisia 97 0 20 1 0 0 4 3 0 16 0 0 97 facebook
3 Croatia 77 0 17 0 0 0 0 10 59 15 2 0 77 facebook
4 Italy 75 0 18 0 0 0 1 13 20 33 3 0 75 facebook
5 Colombia 65 0 15 0 0 0 19 13 14 2 1 0 65 facebook
6 Venezuela 64 0 35 0 0 0 10 6 9 2 1 0 64 facebook
7 Albania 52 0 11 0 0 0 11 5 0 0 0 0 52 facebook
8 Bosnia and Herzegovina 52 0 10 0 0 0 1 5 22 0 0 0 52 facebook
9 France 48 0 18 0 0 0 1 15 11 23 3 0 48 facebook
10 United Kingdom 48 3 82 1 23 0 2 35 41 45 12 2 82 twitter
...


E' stata aggiunta la colonna MaxS, con il nome della colonna con valore più alto nel dataframe importato in R.

A questo punto possiamo passare alla fase grafica. Abbiamo già visto in questo post come creare una mappa politica con le varie regioni colorate. Il procedimeto risulterebbe del tutto analogo; potremmo eventualmente utilizzare anche le library maps e maptools, che forniscono alcune mappe come il sito gadm.
Vediamo invece un metodo alternativo per ottenere una mappa delle nazioni, uscendo dal programma R, e utilizzando le Google Chart API.

Questo è uno strumento messo a disposizione da Google, che ci permette di ottenere dei grafici, semplicemente trasformando i nostri dati in un link.
Andiamo a vedere quali dati ci servono, dalla guida fornita da Google:
- elenco delle nazioni in formato ISO 3166-1
- scelta dei colori da utilizzare
- elenco di valori relativi alle nazioni

Cominciamo dall'ultimo punto. Da R a Google Chart non possiamo specificare il nome del Social Network per ciascuna nazione, ma dobbiamo indicarlo sotto forma di numero. Per far questo, estraiamo la colonna MaxS, trasformiamola in un fattore, e quindi nuovamente in un valore numerico relativo. Prima di far questo però, mettiamo il dataframe in ordine alfabetico in base alla colonna Region:


social <- social[order(social$Region),]

social[1:6,]
Region facebook orkut twitter imeem bebo friendster hi5 myspace last.fm linkedin livejournal xanga max MaxS
72 Afghanistan 8 1 0 0 0 0 1 9 0 0 0 0 9 myspace
7 Albania 52 0 11 0 0 0 11 5 0 0 0 0 52 facebook
65 Algeria 9 0 6 0 0 0 0 1 0 7 0 0 9 facebook
127 Angola 2 6 0 0 0 0 16 2 0 0 0 0 16 hi5
29 Argentina 21 1 13 0 0 0 1 2 11 17 2 0 21 facebook
129 Armenia 1 0 0 0 0 0 0 1 0 0 0 0 1 myspace
...

maxx <- social$MaxS
maxx <- factor(maxx)

maxx
...
Levels: bebo facebook friendster hi5 imeem last.fm linkedin livejournal myspace orkut twitter xanga

maxx <- as.numeric(maxx)

maxx[1:10]
[1] 9 2 2 4 2 9 9 6 11 11


All'interno del vettore maxx ciascun valore indica il livello, nell'ordine prima dichiarato. Ad esempio la prima nazione (Afghanistan) ha livello 9, ossia myspace; e così via.
Adesso dobbiamo rapportare a 100 i livelli, perchè nelle Chart di Google occorre specificare un gradiente di valori. Per far questo utilizziamo la funzione di standardizzazione già studiata in questo post:


library(vegan)
maxx <- decostand(maxx, "max")


Essendo rapportato a 1, moltiplichiamo per 100 e arrotondiamo all'unità; quindi trasformiamo in un vettore:


maxx <- maxx*100
maxx <- round(maxx)
maxx <- as.vector(maxx)
maxx[1:6]
[1] 75 17 17 33 17 75


Ovviamente tutti questi comandi possono essere scritti in un'unica riga. Per motivi "didattici", tendo sempre a separare le singole operazioni.

A questo punto abbiamo i valori da assegnare alle nazioni: sono valori numerici perchè come tali devono essere letti da Google, ma in realtà fanno riferimento ai livelli. Ossia abbiamo che i livelli standardizzati a 100 sono: 8 17 25 33 42 50 58 67 75 83 92 100, che si riferiscono, ciascuno, ai social networks nell'ordine che ho scritto prima.

Ora dobbiamo cercare i codice ISO 3166-1 per ciascuna nazione: ecco il link con la tabella. Trascriviamo ciascun codice per le nazioni che consideriamo, in ordine alfabetico.
Scegliam infine i colori: possiamo decidere di creare un gradiente di colori (Google calcola automaticamente quanti step inserire tra i colori specificati di inizio e fine, tanti quanti sono i livelli), oppure possiamo scegliere 12 colori diversi a nostro piacimento. Utilizziamo ad esempio la seguente sequenza di colori:
- F57DC3 (█ bebo ),
- 54F615 (█ facebook ),
- D14BC7 (█ friendster ),
- 3772E0 (█ hi5 ),
- CD012B (█ imeem ),
- F15239 (█ last.fm ),
- 06D197 (█ linkedin ),
- 7787F0 (█ livejournal ),
- 631E2A (█ myspace ),
- 8AAB5E (█ orkut ),
- C87E99 (█ twitter ),
- FF119B (█ xanga)

Siamo ora in possesso di tutte le informazioni necessarie per creare il grafico con le Google API Chart. Il codice è il seguente:


http://chart.apis.google.com/chart?
cht=t&
chs=440x220&
chtm=world&
chco=FFFFFF,F57DC3,54F615,D14BC7,3772E0,CD012B,F15239,06D197,7787F0,631E2A,8AAB5E,C87E99,FF119B&
chld=AFALDZAOARAMAUATAZBDBYBEBJBOBABWBRBGBFBIKHCMCATDCLCNCOCRCIHRCUCZDKDOECEGSVEREEETFJFIFRGAGMGEDEGHGRGTGNHTHNHKHUINIDIRIQIEILITJMJPJOKZKEKWKGLALVLBLSLYLTMKMGMWMYMLMRMUMXMDMNMAMZMMNANPNLNZNINENGNOOMPKPSPAPGPYPEPHPLPTPRRORURWSASNSLSGSKSIZAKPESLKSDSZSECHSYTWTZTHTGTTTNTRTMUGUAAEGBUSUYUZVEVNYEZMZW&
chd=t:75,17,17,33,17,75,75,50,92,92,67,58,33,17,17,92,92,50,33,17,17,33,92,33,17,92,17,33,17,17,92,50,58,33,33,17,33,17,50,92,8,50,17,17,17,75,50,92,17,33,33,83,33,100,50,58,25,92,92,8,58,17,92,92,17,67,92,92,17,33,50,92,17,17,50,17,17,17,25,17,17,17,92,67,33,17,33,92,17,92,58,92,33,33,92,50,58,58,17,17,17,83,33,25,50,33,75,33,67,17,58,92,17,67,50,50,92,92,50,17,17,17,50,50,92,92,92,42,17,92,17,17,17,92,67,58,92,75,17,17,17,92,17,92,92&
chf=bg,s,EAF7FE


Analizziamolo, come di consueto (ciascun parametro è separato dal successivo da una e-commerciale &):
- link alle API
- cht: tipo di grafico (in questo caso una mappa)
- chs: dimensioni (per questo tipo di grafico, purtroppo, le dimensioni massime sono quelle indicate)
- chtm: quale mappa visualizzare
- chco: sequenza di colori da utilizzare; il primo colore è il bianco, quello dei territori non scelti;
- chld: elenco senza spazi o separatori, dei codici ISO delle nazioni da colorare
- chd=t:: elenco dei livelli di colore. Ad esempio la prima nazione (AF, Afghanistan) ha livello 75, che è il nono tra i possibili valori; quindi verrà colorato col nono colore specificato; e così via
- chf: colore del mare e background

Il risultato è il seguente:



Mappa ingrandita con i comandi width e height del codice HTML:



La mappa è molto piccola, ma il formato che Google mette a disposizione è solo questo; nonostante ciò ho utilizzato questo strumento per "imparare qualcosa di nuovo". Possiamo verificare la corretteza dei colori, usando la legenda che ho scritto sopra e richiamando il dataframe social da leggere.
Se qualche lettore vuole divertirsi a ottenere questa stessa mappa con il metodo descritto nel mio post precedente, sarò ben felice di caricare qui la sua immagine :)

Spero di non aver fatto errori, ma stavolta il codice da usare era davvero tanto!!

Felice anno nuovo a tutti!!

EDIT: Google Chart interpreta male i colori. Sebbene siano correttamente posizionati, non rispettano quelli della legenda di sopra. Non appena riesco a individuare l'errore, correggerò il post.

EDIT n°2: Ho scoperto di recente un'altra fantastica api di google. Si tratta di Google Fusion Table. Accedendo con l'account Google, è possibile caricare i propri file Excel, e visualizzarli con i più comuni strumenti grafici, come istogrammi, grafici a torta, eccetera. Tra questi, è anche disponibile la visualizzazione della mappa, che potete vedere qui sotto.
A seguito dell'analisi che abbiamo fatto, ho creato una tabella contenente un numero corrispondente a ciascun social network in base all'ordine sempre utilizzato (quindi bebo ha valore 1, facebook ha valore 2, e così via). Spostandovi col mouse, potete vedere per ciascuna nazione qual è il social network più utilizzato.



P.S. Non chiedetemi come mai al posto di comparire "USA" a volte compaia "Georgia"!! :)

mercoledì 23 dicembre 2009

Standardizzazione e diamond chart

In questo post vediamo come creare un diamond plot, o diamond chart, o spider chart, o diagramma di Kiviat, ossia uno dei tanti metodi grafici per confrontare diversi parametri tra vari oggetti in studio. Per questo esempio, ho scelto di confrontare tra loro 3 diversi modelli di auto, ottenendo i dati dal sito gaadi.com. I 3 modelli che ho scelto sono:
- Fiat Grande Punto 1.4 Emotion Pack
- Lamborghini Gallardo Spyder Base
- Mercedes Benz SLK Class SLK 350

I parametri che andremo a confrontare invece sono:
- Kilometri per litro (Mileage)
- Capacità del serbatoio (FuelTank)
- Velocità massima (MaxSpeed)
- Cilindrata (Displacement)
- Potenza in cavalli (PowerPS)

Per ottenere il diamond plot che cerchiamo, installiamo e carichiamo la libreria plotrix. Qui sono disponibili due funzioni tra loro simili, sebbene richiedano due diverse sintassi. Cominciamo a studiare la funzione radial.plot().
Prima di procedere, però dobbiamo standardizzare i dati. La standardizzazione si rende necessaria dal momento che i valori coprono un range troppo elevato per poter essere visualizzato in un unico grafico (da circa 5 a circa 5000, la rappresentazione grafica sarebbe quasi illeggibile). Perciò dobbiamo standardizzare i dati, e utilizzeremo la funzione decostand della libreria vegan. Per prima cosa creiamo una matrice con i dati:


Mileage <- c(11.8, 5.4, 7)
FuelTank <- c(45, 90, 70)
MaxSpeed <- c(180, 315, 252)
Displacement <- c(1368, 4961, 3498)
PowerPS <- c(90, 520, 275.8)

dati1 <- matrix(c(Mileage, FuelTank, MaxSpeed, Displacement, PowerPS), nrow=3)

dati1
[,1] [,2] [,3] [,4] [,5]
[1,] 11.8 45 180 1368 90.0
[2,] 5.4 90 315 4961 520.0
[3,] 7.0 70 252 3498 275.8


A questo punto carichiamo la libreria vegan, e utilizziamo la funzione decostand. Leggendo l'help relativo (help(decostand)), osserviamo che sono a disposizione numerosi metodi di standardizzazione. In questo caso ho scelto di utilizzare il metodo "max": per ogni variabile, il valore più alto viene rapportato ad 1, e gli altri vengono calcolati con una semplice proporzione; dobbiamo specificare MARGIN=2, in modo tale che vengano standardizzati tra loro i valori di ogni colonna (scegliendo MARGIN=1 avrebbe standardizzato tra loro i valori di ogni riga):


library(vegan)
dati1ST <- decostand(dati1, "max", MARGIN=2)

dati1ST
[,1] [,2] [,3] [,4] [,5]
[1,] 1.0000000 0.5000000 0.5714286 0.2757509 0.1730769
[2,] 0.4576271 1.0000000 1.0000000 1.0000000 1.0000000
[3,] 0.5932203 0.7777778 0.8000000 0.7050998 0.5303846
attr(,"decostand")
[1] "max"


Con la matrice così costruita, possiamo utilizzare la funzione radial.plot():


radial.plot(dati1ST,labels=c("Mileage", "FuelTank", "MaxSpeed", "Displacement", "PowerPS") ,rp.type="p", main="Diamond plot auto",line.col=2:4,show.grid=FALSE,lwd=1:3, radial.lim=c(0,1.2))


Prima di osservare il risultato, analizziamo la sintassi:
- dati1ST: specificare il nome della matrice contenente i dati; la matrice deve essere organizzata in modo che le variabili considerate (kilometraggio, capacità, etc.) siamo disposte in colonne;
- labels: elenco dei nomi delle variabili, che verrano sistemati sul grafico
- rp.type="p": specifica di disegnare le linee che congiungono i vari punti
- main: titolo della finestra
- line.col: colori da utilizzare
- show.grid=F: non mostra la griglia
- lwd: spessore delle linee
- radial.lim: range di valori, in questo caso da 0 a 1.2 (essendo il nostro valore massimo pari a 1)

Ecco il risultato (eventualmente si può aggiungere una legenda):



L'altra funzione a disposizione per i grafici a diamante, è la funzione diamondplot(). La sintassi è leggermente diversa: innanzitutto i dati devono essere scritti sottoforma di dataframe, e non di matrice come sopra. Inoltre le colonne sono rappresentate dagli oggetti in studio, e non dalle variabili come per la funzione precedente. Fatte queste precisazioni, ecco il codice che ho utilizzato per ottenere un grafico molto simile al precedente:


Mileage <- c(11.8, 5.4, 7)
FuelTank <- c(45, 90, 70)
MaxSpeed <- c(180, 315, 252)
Displacement <- c(1368, 4961, 3498)
PowerPS <- c(90, 520, 275.8)

dati2 <- matrix(c(Mileage, FuelTank, MaxSpeed, Displacement, PowerPS), ncol=3, byrow=T)
dati2
[,1] [,2] [,3]
[1,] 11.8 5.4 7.0
[2,] 45.0 90.0 70.0
[3,] 180.0 315.0 252.0
[4,] 1368.0 4961.0 3498.0
[5,] 90.0 520.0 275.8

library(vegan)
dati2ST <- decostand(dati2, "max", MARGIN=1)

dati2ST
[,1] [,2] [,3]
[1,] 1.0000000 0.4576271 0.5932203
[2,] 0.5000000 1.0000000 0.7777778
[3,] 0.5714286 1.0000000 0.8000000
[4,] 0.2757509 1.0000000 0.7050998
[5,] 0.1730769 1.0000000 0.5303846
attr(,"decostand")
[1] "max"

dati2DF <- as.data.frame(dati2ST)
colnames(dati2DF) <- c("Fiat", "Lamborghini", "Mercedes")
rownames(dati2DF) <- c("Mileage", "FuelTank", "MaxSpeed", "Displacement", "PowerPS")
dati2DF
Fiat Lamborghini Mercedes
Mileage 1.0000000 0.4576271 0.5932203
FuelTank 0.5000000 1.0000000 0.7777778
MaxSpeed 0.5714286 1.0000000 0.8000000
Displacement 0.2757509 1.0000000 0.7050998
PowerPS 0.1730769 1.0000000 0.5303846

library(plotrix)
diamondplot(dati2DF, bg=gray(0.6), col=cm.colors, name="Diamond plot auto")



lunedì 21 dicembre 2009

Choropleth map in R: coloriamo le mappe geografiche

Dopo una lunga pausa, finalmente torno ad occuparmi del mio blog!
In questo post vedremo come creare una choropleth map, ossia una mappa geografica con le diverse regioni colorate in base ad un certo valore numerico.
Per far questo, ecco ciò di cui abbiamo bisogno:
** libreria sp
** libreria Hmisc
** la mappa bianca: su questo sito è possibile ottenere le mappe in un formato leggibile da R
** i dati da plottare: per questo esempio i dati sono stati presi dal sito Istat

Mettiamoci al lavoro.

Innanzitutto vediamo come caricare in R la mappa bianca, da colorare successivamente. Identifichiamo la nazione di nostro interesse sul sito Gadm, copiamo il link e importiamolo in R con questo codice:


con <- url("http://gadm.org/data/rda/ITA_adm1.RData")
print(load(con))
close(con)


Il primo problema che si pone è quello di capire in che modo verranno colorate le regioni. Con la funzione spplot della library(sp), sarà sufficiente creare un vettore di colori, ciascuno dei quali verrà assegnato alla singola regione. In che ordine vengono presentate le regioni? Richiamiamolo con il seguente codice:


summary(gadm$NAME_1)
Abruzzo Apulia Basilicata
1 1 1
Calabria Campania Emilia-Romagna
1 1 1
Friuli-Venezia Giulia Lazio Liguria
1 1 1
Lombardia Marche Molise
1 1 1
Piemonte Sardegna Sicily
1 1 1
Toscana Trentino-Alto Adige Umbria
1 1 1
Valle d'Aosta Veneto
1 1


Osservate che le regioni sono state messe in ordine alfabetico, ad eccezione della Puglia che viene chiamata Apulia. Quando assegneremo i colori alle regioni, dobbiamo seguire questo ordine.

A questo punto dal sito Istat cerchiamo i valori da analizzare. In questo esempio, andremo a dividere le regioni italiana in base alla percentuale di donne nubili. Cerchiamo i valori sul sito e importiamoli in R (clickate su view plain e allargata la finestra per visualizzare correttamente questa tabella):


data <- read.table("clipboard", header=T)

data
Regioni Maschi_Celibi Totale_Maschi Femmine_Nubili Totale_Femmine Maschi_Femmine
1 Abruzzo 284385 648680 242027 685995 1334675
2 Puglia 883027 1979254 809566 2100448 4079702
3 Basilicata 130276 289275 112162 301326 590601
4 Calabria 452617 978789 398402 1029920 2008709
5 Campania 1324589 2820078 1207237 2992884 5812962
6 Emilia_Romagna 944689 2109482 795045 2228497 4337979
7 Friuli_Venezia_Giulia 261186 596265 212969 634671 1230936
8 Lazio 1241962 2703994 1126148 2922716 5626710
9 Liguria 320014 767057 273659 848007 1615064
10 Lombardia 2148215 4762370 1806394 4980306 9742676
11 Marche 334044 763741 278521 805837 1569578
12 Molise 68308 156036 58492 164759 320795
13 Piemonte 917932 2149373 761839 2283198 4432571
14 Sardegna 399198 819518 349360 851483 1671001
15 Sicilia 1113415 2433605 1009358 2604194 5037799
16 Toscana 757715 1787668 649483 1920150 3707818
17 Trentino_Alto_Adige 251109 500811 215719 517846 1018657
18 Umbria 182907 431313 155483 462909 894222
19 Valle_D_Aosta 29314 62451 23671 64614 127065
20 Veneto 1083455 2392663 905460 2492885 4885548


Dal sito Istat abbiamo importato questa tabella con la colonna di donne nubili e la colonna di donne totali, per regione. A questo punto calcoliamo la percentuale relativa per ogni regione di donne nubili:


perc_nubili <- data$Femmine_Nubili / data$Totale_Femmine * 100


Dividiamo adesso questo vettore in un fattore a 5 livelli in base alla percentuale (così da ottenere ad esempio le classi: "altissima percentuale", "alta percentuale", "media percentuale", "bassa percentuale", "bassissima percentuale", di donne nubili). Utilizziamo la funzione cut2 della library(Hmisc):


library(Hmisc)
lev <- cut2(perc_nubili, g=5)


Il parametro g=5 indica proprio il numero di classi. Automaticamente la funzione calcola i cut-off in modo tale da dividere i 20 valori in 5 gruppi, ciascuno con 4 valori.

Ora dobbiamo creare un vettore contenente l'elenco delle classi, relativi alle regioni. Ossia nell'ordine delle regioni visto prima, dobbiamo ottenere un vettore che assegni il livello a ciascuna regione (ad esempio la prima regione, Abruzzo, appartiene alla quarta classe). Utilizziamo la seguente stringa:


nubili <- as.character(as.numeric(lev))

nubili
[1] "2" "4" "4" "4" "5" "3" "1" "4" "1" "3" "2" "2" "1" "5" "5" "2" "5" "1" "3" "3"


Adesso creiamo un fattore dei livelli così ottenuti, e lo aggiungiamo alla variabile gadm, contenente i dati della mappa:


gadm$nubili <- as.factor(nubili)


Creiamo infine un vettore contenente i 5 colori (tanti quanti sono i livelli creati):


col <- rev(heat.colors(length(levels(lev))))


heat.colors permette di ottenere una palette di colori sulle tonalità rosse; altre palette già impostate in R sono:
- rainbow
- heat.colors
- terrain.colors
- topo.colors
- cm.colors

A questo punto possiamo colorare la nostra mappa. Utilizziamo la funzione spplot della library sp. La sintassi che ho utilizzato è la seguente:
- gadm: è il nome della variabile contenente i dati della mappa
- "nubili": è il vettore contenente i livelli per ciascuna regione
- col.regions: è il vettore contenente i colori da utilizzare
- main: è il titolo della finestra


spplot(gadm, "nubili", col.regions=col, main="Classi regionali di donne nubili")


Il risultato ottenuto è il seguente:


Come si può vedere dalla mappa, le regioni a più bassa percentuale di donne nubili (colore giallo chiaro, classe 1) sono: Umbria, Friuli Venezia Giulia, Piemonte e Liguria; le regioni a più alta percentuale (colore rosso, classe 5) invece sono: Sicilia, Sardegna, Campania, Trentino Alto Adige.

Sul sito Gadm, sono anche disponibili mappe più dettagliate, che dividono l'Italia in base alle province, e anche in base alle città, in modo da ottenere choropleth map ancora più dettagliate.

Per par condicio, ecco il codice e la mappa relativi alla percentuale di maschi celibi:


perc_celibi <- data$Maschi_Celibi / data$Totale_Maschi * 100
data2 <- data.frame(data, perc_celibi)
lev <- cut2(perc_celibi, g=5)

con <- url("http://gadm.org/data/rda/ITA_adm1.RData")
print(load(con))
close(con)

celibi <- as.character(as.numeric(lev))
gadm$celibi <- as.factor(celibi)
col <- rev(heat.colors(length(levels(lev))))
spplot(gadm, "celibi", col.regions=col, main="Classi regionali di uomini celibi")