Grafici avanzati con R

Oggi continuiamo il nostro viaggio all’interno dell’incredibile mondo di R. Fin’ora abbiamo utilizzato alcune delle potenti funzioni matematiche e statistiche, relegando la parte grafica solo all’ultima fase. Ebbene, con l’articolo di oggi, invece, andremo decisamente più a fondo nelle funzioni grafiche (notevoli, va detto) di R.

Per l’esperimento di oggi utilizzeremo la serie di dati UAH aggiornati ad Aprile 2009. Dal momento che rispetto alle tecniche utilizzate in precedenza, faremo ricorso a metodologie più complesse, nei limiti del possibile cercherò di dettagliare ogni passaggio. R è un software open source, disponibile sia per Linux che per Microsoft Windows.

Oltre alle spiegazioni in questo articolo, troverete anche delle note all’interno del codice, nella seguente forma:

# Questo è un commento

Nel momento in cui andrete a inserire il codice sottostante, l’interprete di R ignorerà le righe di commento. Con l’esperimento di oggi otterremo il successivo grafico, nel quale oltre ad una rappresentazione evoluta della serie di dati, troviamo rappresentato il trend e una funzione di smussamento dei dati (LOWESS o LOESS, ovvero Local WEighted Scatterplot Smoothing).

Qui trovate i dati, da qui invece potete scaricare l’intero script.

Iniziamo con alcune operazioni preliminari:

####################
# INIZIALIZZAZIONE #
####################

library(tseries)
options(digits=6) # APPROSSIMAZIONE ALLA 6° CIFRA DECIMALE

Per il momento nulla di particolare, abbiamo caricato la libreria dedicata alle serie storiche (l’abbiamo già utilizzata in passato qui) e abbiamo impostato la precisione di calcolo fino alla sesta cifra decimale.

#############################
# INPUT E PREPARAZIONE DATI #
#############################
dati <- read.table ("http://www.climatemonitor.it/Download/R/UAH.csv",header=FALSE,sep="", dec=".",col.names=c("Anno","Mese","Anomalia"))
UAH_mensile <- data.frame(dati)
attach (UAH_mensile)
c <- nrow (UAH_mensile)
ultimo_mese <- UAH_mensile$Mese[c]
ultimo_anno <- UAH_mensile$Anno[c]
ultimo_valore <- Anomalia[c]
UAH <- ts (Anomalia,start=c(1978,12),frequency=12)

La fase di input e preparazione dei dati è decisamente più complessa rispetto a quanto abbiamo visto fino ad oggi. Con il comando read.table andiamo a leggere il nostro file, indicando a R che non abbiamo un header, che i decimali sono individuati dal punto e che le nostre colonne (in tutto tre) desideriamo chiamarle “Anno”, “Mese” e “Anomalia”. Fin qui la nostra serie di dati. Per lavorare in modo più agevole d’ora in poi creeremo un cosiddetto data frame, ovvero un contenitore entro il quale inseriremo tutti i nostri dati, quindi non solo la nostra serie delle anomalie.

Per rendere più snelli i calcoli anche in futuro, individuiamo tre misurazioni importanti, tramite il comando nrow andiamo a contare il numero di osservazioni presenti nella nostra serie storica, in particolare rimarrà memorizzata nel vettore l’ultima osservazione disponibile. In questo modo dal vettore c appositamente creato, potremo estrarre i dettagli dell’ultima osservazione disponibile (mese, anno e anomalia). Infine creiamo la nostra serie storica UAH a partire dalla colonna “Anomalia”.

######################
# GRAFICO PRINCIPALE #
######################
plot (UAH,type="s",col="darkgrey",xlab="",ylab = "°C - Anomalia termica",xlim=c(1979, 2010), ylim=c(-0.6, 0.80))
lines (UAH,type="h",col="lightgrey") # GRAFICO DI FONDO
abline (h=0,col="darkgrey")

Con questa serie di comandi, invece, creiamo il grafico principale che consiste nella serie di dati (anomalie termiche) riprodotta due volte. Notate che con il comando plot creiamo solo l’andamento della serie, con il comando lines invece andiamo a creare un grafico a barre che funge da riempitivo al primo (in pratica lo utilizziamo come sfondo colorato). Infine aggiungiamo la linea base.

##########
# TITOLI #
##########
mtext ("Climate Monitor - www.climatemonitor.it",side=1,line=2,adj=1,cex=0.7)
mtext ("Dati: http://vortex.nsstc.uah.edu/public/msu/t2lt/tltglhmam_5.2",side=1,line=3,adj=1,cex=0.7)
title (main="Anomalie termiche - Dati UAH",cex=0.9)

Inseriamo titoli e commenti. Qui è tutto piuttosto semplice. Mi soffermo brevemente sul comando mtext, come vedete ha qualche parametro. Rimandandovi al manuale ufficiale di R per un dettaglio maggiore, richiamatelo così:

help(mtext)

vediamo che side stabilisce il lato nel quale comparirà il testo (1=in basso, 2=a sinistra, 3=in alto), con line selezioniamo il numero di riga (ne abbiamo 3 in tutto a disposizione), con adj modifichiamo la formattazione (1=a destra, 0=a sinistra) e infine cex per determinare la grandezza del testo.

###############
# REGRESSIONE #
###############
fit <- lm (UAH~Anno)
a <- coef(fit)[1]
b <- coef(fit)[2]
Anno1 <- min (Anno)
Anno2 <- max (Anno)
y1 <- a+b*Anno1
y2 <- a+b*Anno2
x_coord <- c (Anno1,Anno2)
y_coord <- c (y1,y2)
lines (x_coord,y_coord,type="l",col="red")

Arriviamo al cuore dei calcoli. Come avrete già notato di tratta di una semplice regressione lineare che ci consente di individuare i valori tipici a,b e coefficiente angolare. Infine con alcuni passaggi matematici otteniamo le coordinate della retta di regressione che poi possiamo rappresentare tramite il comando lines. Scendendo più in dettaglio, creiamo la serie fit che altro non è se non la regressione lineare calcolata tramite lm della serie UAH rispetto al tempo.

##########
# LOWESS #
##########
lines (lowess(UAH,f=.1),col="blue",lwd=2)

Infine creiamo una serie smussata con il comando lowess.
Il grafico è sicuramente perfettibile, in un prossimo intervento vedremo di integrare ulteriori informazioni all’interno dello stesso grafico (senza appesantirne la rappresentazione). Potremmo valutare di inserire metodi alternativi di smussamento. Ripeto, il grafico è in continua evoluzione, l’obiettivo è quello di arrivare ad una forma definitiva che possa diventare lo standard per le nostre rappresentazioni. Proprio per questo motivo, rimango in attesa delle vostre preziose segnalazioni e dei vostri suggerimenti.

Una raccomandazione: non limitatevi ad effettuare il copia e incolla del codice, ma scrivete riga per riga cercando di pensare al singolo comando ed eventualmente a possibili migliorie. Infine, se decideste di copiare e incollare il codice (piuttosto che digitarlo) fate attenzione all’ultima riga copiata, spesso alcuni sistemi operativi non inviano il ritorno a capo, con il rischio che l’ultima riga di comando non venga letta.

Bibliografia

– Cleveland, W.S. (1979) “Robust Locally Weighted Regression and Smoothing Scatterplots,” Journal of the American Statistical Association, 74, 829–836.
– Cleveland, W.S., Devlin, S.J. (1988) “Locally-Weighted Regression: An Approach to Regression Analysis by Local Fitting,” Journal of the American Statistical Association, 83, 596–610.
– Kelly O’Day, ProcessTrends.com