Grafici avanzati con R – Approfondimento

Visto l’interesse suscitato dai passati articoli, abbiamo deciso di andare avanti nell’approfondire nuove tecniche di rappresentazione grafica ed elaborazione dei dati. Siamo ormai arrivati alla versione 1.1 del nostro grafico delle temperature. Fino a ieri, eravamo in grado di rappresentare la serie storica (scaricandola da internet e rendendola fruibile a R), e di sovraimporre la serie smussata.

Abbiamo utilizzato Lowess come tecnica per smussare i dati (sebbene sappiamo che Lowess non sia uno smussatore tout-court). Oggi vedremo come aggiungere altri elementi utili al grafico quali:

• una retta di regressione;
• la lettura e l’evidenziazione dell’ultimo dato climatico disponibile;
• la visualizzazione del trend decadale.

Nel caso non l’abbiate ancora fatto, vi consigliamo di procurarvi il software statistico R, disponibile sia per Windows che per Linux. E’ un software open source, ormai ben consolidato ed affidabile negli algoritmi e nella precisione di calcolo (viene abitualmente utilizzato in campo accademico, come alternativa ai software proprietari dal costo spesso proibitivo).

Prima di iniziare, ringraziamo il sito Processtrends.com per i preziosi consigli e suggerimenti.

Questa nuova versione del grafico (la 1.1) corregge alcuni errori presenti nella versione precedente che, quindi, è da ritenersi superata e da ignorare. Analizziamo il codice, passo a passo.


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

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


In questa fase istruiamo R con due comandi base: da un lato fissiamo l’arrotondamento da noi desiderato, dall’altro lato carichiamo la libreria relativa alle serie storiche, perchè ci tornerà utile più avanti.
Nel caso in cui non abbiate installato la libreria, vi ricordiamo rapidamente il comando da utilizzare:

install.packages("tseries")

Procediamo adesso con una fase molto importante:

#############################
# INPUT E PREPARAZIONE DATI #
#############################

dati<- "http://www.remss.com/data/msu/monthly_time_series/RSS_Monthly_MSU_AMSU_Channel_TLT_Anomalies_Land_and_Ocean_v03_2.txt"

RSS_mensile <- read.table(dati,skip = 3,sep = "",dec=".",row.names = NULL,header = FALSE,as.is = T,colClasses = c(rep("numeric",3),rep("NULL", 8)),comment.char = "#",na.strings = c("*", "-",-99.9, -999.9),col.names = c("Anno", "Mese", "RSS_anom", rep("",8)))

anno_fraz <- RSS_mensile$Anno + (RSS_mensile$Mese-1)/12

RSS_data_frame<-data.frame(RSS_mensile,anno_fraz)
attach(RSS_data_frame)
c<-nrow(RSS_data_frame)
ultimo_mese <- RSS_data_frame$Mese[c]
ultimo_anno <- RSS_data_frame$Anno[c]
ultimo_dato <- RSS_anom[c]

RSS <- ts(RSS_anom,start=c(1979,1),frequency=12)

Come possiamo vedere, questo segmento di codice scarica i dati relativi a RSS, allocando l’intera tabella all’interno della variabile dati. Tuttavia a noi interessano soltanto le prime tre colonne che intitoliamo “Anno”, “Mese” e “RSS_anom”.

Richiamiamo la vostra attenzione su questi due importanti comandi:

1. colClasses = c(rep("numeric",3),rep("NULL", 8)): in questo modo suggeriamo a R di interpretare le prime tre colonne come dato numerico, le restanti invece non ci interessano;
2. comment.char = "#",na.strings = c("*", "-",-99.9, -999.9): con questo invece ordiniamo a R di convertire in stringa vuota tutte quelle celle della tabella che contengano uno dei caratteri elencati (asterisco, “-” ecc.). Questo è un comando utilizzabile sempre, qualora si importi una tabella in R.

I dati forniti sono su base mensile, per consentirne però una più semplice rappresentazione grafica, andiamo a calcolare l’anno come frazione di tempo e inseriamo i dati nel vettore anno_fraz. A questo punto creiamo una nuova matrice di dati che contenga tutte le colonne della serie RSS_mensile più il vettore colonna appena creato anno_fraz. Abbiamo quasi terminato, dal punto di vista della formattazione dei dati, non ci rimane altro che individuare alcune grandezze utili, come ultimo_mese, ultimo_anno, ultimo_dato.

Creiamo quindi la serie storica RSS (per ulteriori informazioni su questo passaggio, vi invitiamo a leggere gli articoli già pubblicati).

######################
# GRAFICO PRINCIPALE #
######################


plot(anno_fraz,RSS_anom,type="s",col="grey",xlab="",ylab = "°C - Anomalia termica",xlim=c(1979, 2010), ylim=c(-0.6,max(RSS)),cex.axis=0.95,cex.lab=0.95)
lines(RSS,type="h",col="lightgrey") # GRAFICO DI FONDO
abline(h=0,col="darkgrey")

Questa sezione non la commenteremo, in quanto ormai l’abbiamo affrontata più volte, ci concentriamo invece sulla regressione e sulla visualizzazione dei relativi dati.


###############
# REGRESSIONE #
###############


lm_fit <- lm(RSS_anom~anno_fraz)
a <- coef(lm_fit)[1]
b <- coef(lm_fit)[2]
yr1 <- min(anno_fraz)
yr2 <- max(anno_fraz)
y1 <- a+b*yr1
y2 <- a+b*yr2
x_val <- c(yr1,yr2)
y_val <- c(y1,y2)
lines(x_val,y_val,type="l",col="red")
b <- signif(b, 3)

Raccogliendo alcuni consigli, compresi quelli del gentile lettore Tommaso, individuiamo la retta di regressione con il comando lm() e volendo possiamo rappresentare la retta con il comando:

abline(lm_fit, col=”red”, lwd=2)

Nel caso volessimo alleggerire l’elaborazione e quindi utilizzare il comando appena descritto, sarebbe sufficiente anteporre il simbolo # alle righe immediatamente successive al comando lm_fit <- lm(RSS_anom~anno_fraz).

Per il momento, tuttavia, vi consigliamo di non escludere alcuna riga dall’elaborazione, in quanto le variabili individuate ci servono per poter visualizzare le informazioni necessarie:


##################
# ULTIMA LETTURA #
##################

points(yr2, ultimo_dato, pch=20, col = "blue")
points(1995, -0.5, pch=20, col = "blue")
note2 <- paste("Ultima osservazione ",ultimo_mese, "/", ultimo_anno, " = ", ultimo_dato," °C")
text(1995, -0.5, note2, pos = 4, col = "blue", cex = 0.7)

#########
# TREND #
#########

note3 <- paste("Trend = ",b," °C / anno ")
text(1995,-0.4, note3,pos=4, col = "red", cex = 0.7)


A questo punto il codice è abbastanza self-explanatory. Nella variabile yr2 abbiamo registrato l’ultima data disponibile nella serie. Andiamo a visualizzarla con il comando points alle coordinate, quindi, di yr2 e ultimo_dato. Una nota sul calcolo del trend: abbiamo utilizzato la variabile interna b, ovvero il coefficiente angolare della retta di regressione e l’abbiamo visualizzato con il comando text.

Ecco il risultato:

Vi lasciamo con un “compito” da svolgere a casa 🙂

Vi proponiamo una variante del grafico appena esaminato. Tale variante può essere utile per visualizzare in linea il dato dell’anomalia termica (o per lo sviluppo di un widget web). A voi la sperimentazione: