Para el congreso Argentina Ambiental 2012 he estado terminando un resumen extendido, y me tiene muy contento el hecho de aplicar estadística dentro del ámbito de los sistemas ambientales.
Para este congreso realicé un pronóstico de CO y NOx para las ciudades de Buenos Aires, Bogotá y Santiago, comparando los resultados de una regresión y redes neuronales. A continuación presentaré dos gráficos que comparan los resultados de forma muy simple y didáctica.
Entonces, con las técnicas que no voy a presentar acá genere las regresiones y redes neuronales, y obtuve un CO y NOx pronosticado con cada técnica para cada ciudad. Pensé en graficar con ggplot2 aplicando un factor específico, y se me ocurrió en excel clasificar las el valor absoluto de la diferencia de los pronósticado con lo original, según percentiles del original, así:
abs(diferencia)<P(20)=1; es decir, el valor absoluto de la diferencia debe cae dentro del percentil 20 del valor original para que de el valor de 1, y entonces generé la escala
abs(diferencia)
Y la fórmula en Excel 2010 es la siguiente:
SI(ABS(H2-S2)<PERCENTIL.INC($H$2:$H$3379;0.2);1;SI(Y(ABS(H2-S2)>PERCENTIL.INC($H$2:$H$3379;0.2);(ABS(H2-S2)<PERCENTIL.INC($H$2:$H$3379;0.4)));0.8;SI(Y(ABS(H2-S2)>PERCENTIL.INC($H$2:$H$3379;0.4);(ABS(H2-S2)<PERCENTIL.INC($H$2:$H$3379;0.6)));0.6;SI(Y(ABS(H2-S2)>PERCENTIL.INC($H$2:$H$3379;0.6);(ABS(H2-S2)<PERCENTIL.INC($H$2:$H$3379;0.8)));0.4;SI(ABS(H2-S2)>PERCENTIL.INC($H$2:$H$3379;0.8);0.2;0)))))
y ahí tenemos nuestro factor
Luego:
setwd("C:/r")
library(openair)
#Uso la librería openair ya que es ideal para trabajar con datos en formato típico
#attach es para manejar directamente las variables, pero acá no es necesario
aire1<-import("rpud-co.csv",sep=";",date.name="fecha",time.name="fecha" )
attach(aire1)
names(aire1)
aire2<-import("rpud-nox.csv",sep=";",date.name="fecha",time.name="fecha" )
attach(aire2)
names(aire2)
aire3<-import("rbai-co.csv",sep=";",date.name="fecha",time.name="fecha" )
attach(aire3)
names(aire3)
aire4<-import("rbai-nox.csv",sep=";",date.name="fecha",time.name="fecha" )
attach(aire4)
names(aire4)
aire5<-import("rbog-co.csv",sep=";",date.name="fecha",time.name="fecha" )
attach(aire5)
names(aire5)
aire6<-import("rbog-nox.csv",sep=";",date.name="fecha",time.name="fecha" )
attach(aire6)
names(aire6)
library(ggplot2)
Y ahora podemos graficar
Referencias:
Para este congreso realicé un pronóstico de CO y NOx para las ciudades de Buenos Aires, Bogotá y Santiago, comparando los resultados de una regresión y redes neuronales. A continuación presentaré dos gráficos que comparan los resultados de forma muy simple y didáctica.
Entonces, con las técnicas que no voy a presentar acá genere las regresiones y redes neuronales, y obtuve un CO y NOx pronosticado con cada técnica para cada ciudad. Pensé en graficar con ggplot2 aplicando un factor específico, y se me ocurrió en excel clasificar las el valor absoluto de la diferencia de los pronósticado con lo original, según percentiles del original, así:
abs(diferencia)<P(20)=1; es decir, el valor absoluto de la diferencia debe cae dentro del percentil 20 del valor original para que de el valor de 1, y entonces generé la escala
abs(diferencia)
- <P(20)=1
- P(20)-P(40)=0.8
- P(40)-P(60)=0.6
- P(60)-P(80)=0.4
- >P(80)=0.2
Y la fórmula en Excel 2010 es la siguiente:
SI(ABS(H2-S2)<PERCENTIL.INC($H$2:$H$3379;0.2);1;SI(Y(ABS(H2-S2)>PERCENTIL.INC($H$2:$H$3379;0.2);(ABS(H2-S2)<PERCENTIL.INC($H$2:$H$3379;0.4)));0.8;SI(Y(ABS(H2-S2)>PERCENTIL.INC($H$2:$H$3379;0.4);(ABS(H2-S2)<PERCENTIL.INC($H$2:$H$3379;0.6)));0.6;SI(Y(ABS(H2-S2)>PERCENTIL.INC($H$2:$H$3379;0.6);(ABS(H2-S2)<PERCENTIL.INC($H$2:$H$3379;0.8)));0.4;SI(ABS(H2-S2)>PERCENTIL.INC($H$2:$H$3379;0.8);0.2;0)))))
y ahí tenemos nuestro factor
Luego:
setwd("C:/r")
library(openair)
#Uso la librería openair ya que es ideal para trabajar con datos en formato típico
#attach es para manejar directamente las variables, pero acá no es necesario
aire1<-import("rpud-co.csv",sep=";",date.name="fecha",time.name="fecha" )
attach(aire1)
names(aire1)
aire2<-import("rpud-nox.csv",sep=";",date.name="fecha",time.name="fecha" )
attach(aire2)
names(aire2)
aire3<-import("rbai-co.csv",sep=";",date.name="fecha",time.name="fecha" )
attach(aire3)
names(aire3)
aire4<-import("rbai-nox.csv",sep=";",date.name="fecha",time.name="fecha" )
attach(aire4)
names(aire4)
aire5<-import("rbog-co.csv",sep=";",date.name="fecha",time.name="fecha" )
attach(aire5)
names(aire5)
aire6<-import("rbog-nox.csv",sep=";",date.name="fecha",time.name="fecha" )
attach(aire6)
names(aire6)
library(ggplot2)
Y ahora podemos graficar
> p7 <-ggplot(aire6)+aes(x=bog_nox_regpro, y=CO)+geom_point(binwidth = 1, aes(size = factor(bog_dif_co_reg),colour = factor(bog_dif_co_reg)))+scale_colour_brewer(type="seq", palette=1)+annotate("text", x = 150,y = 15, label = "R2 = 0.45", colour = "black", size=8)+ opts(title="NOx Bogotá Regresión")+opts(legend.position="none")+ xlab("NOx Pronosticado (ppb)")+ylab("NOx Original (ppb)") > p10 <- ggplot(aire6)+aes(x=bog_nox_neupro, y=CO)+geom_point(binwidth = 1, aes(size = factor(bog_dif_co_neu),colour = factor(bog_dif_co_neu)))+scale_colour_brewer(type="seq", palette=1)+annotate("text", x = 150,y = 15, label = "R2 = 0.45", colour = "black", size=8)+ opts(title="NOx Bogotá Neuronal")+opts(legend.position="none")+ xlab("NOx Pronosticado (ppb)")+ylab("NOx Original (ppb)") > p8 <-ggplot(aire4)+aes(x=bai_nox_reg_pro, y=ba_nox_original)+geom_point(binwidth = 1, aes(size = factor(bai_dif_nox_reg),colour = factor(bai_dif_nox_reg)))+scale_colour_brewer(type="seq", palette=1)+annotate("text", x = 140,y = 300, label = "R2 = 0.40", colour = "black", size=8)+ opts(title="NOx Buenos Aires Regresión")+opts(legend.position="none")+ xlab("NOx Pronosticado (ppb)")+ylab("NOx Original (ppb)") > p11 <-ggplot(aire4)+aes(x=bai_nox_neu_pro, y=ba_nox_original)+geom_point(binwidth = 1, aes(size = factor(bai_dif_nox_neu),colour = factor(bai_dif_nox_neu)))+scale_colour_brewer(type="seq", palette=1)+annotate("text", x = 100,y = 300, label = "R2 = 0.33", colour = "black", size=8)+ opts(title="NOx Buenos Aires Neuronal")+opts(legend.position="none")+ xlab("NOx Pronosticado (ppb)")+ylab("NOx Original (ppb)") > p9 <-ggplot(aire2)+aes(x=p_nox_reg_pro, y=p24nox)+geom_point(binwidth = 1, aes(size = factor(p_dif_nox_reg),colour = factor(p_dif_nox_reg)))+scale_colour_brewer(type="seq", palette=1)+annotate("text", x = 170,y = 400, label = "R2 = 0.73", colour = "black", size=8)+ opts(title="NOx Santiago Regresión")+opts(legend.position="none")+ xlab("NOx Pronosticado (ppb)")+ylab("NOx Original (ppb)") > p12 <- ggplot(aire2)+aes(x=p_nox_neu_pro, y=p24nox)+geom_point(binwidth = 1, aes(size = factor(p_dif_nox_neu),colour = factor(p_dif_nox_neu)))+scale_colour_brewer(type="seq", palette=1)+annotate("text", x = 170,y = 400, label = "R2 = 0.80", colour = "black", size=8)+ opts(title="NOx Santiago Neuronal")+opts(legend.position="none")+ xlab("NOx Pronosticado (ppb)")+ylab("NOx Original (ppb)"
> p7 <-ggplot(aire6)+aes(x=bog_nox_regpro, y=CO)+geom_point(binwidth = 1, aes(size = factor(bog_dif_co_reg),colour = factor(bog_dif_co_reg)))+scale_colour_brewer(type="seq", palette=1)+annotate("text", x = 150,y = 15, label = "R2 = 0.46", colour = "black", size=8)+ opts(title="NOx Bogotá Regresión")+opts(legend.position="none")+ xlab("NOx Pronosticado (ppb)")+ylab("NOx Original (ppb)")
> p10 <- ggplot(aire6)+aes(x=bog_nox_neupro, y=CO)+geom_point(binwidth = 1, aes(size = factor(bog_dif_co_neu),colour = factor(bog_dif_co_neu)))+scale_colour_brewer(type="seq", palette=1)+annotate("text", x = 150,y = 15, label = "R2 = 0.45", colour = "black", size=8)+ opts(title="NOx Bogotá Neuronal")+opts(legend.position="none")+ xlab("NOx Pronosticado (ppb)")+ylab("NOx Original (ppb)")
> p8 <-ggplot(aire4)+aes(x=bai_nox_reg_pro, y=ba_nox_original)+geom_point(binwidth = 1, aes(size = factor(bai_dif_nox_reg),colour = factor(bai_dif_nox_reg)))+scale_colour_brewer(type="seq", palette=1)+annotate("text", x = 140,y = 300, label = "R2 = 0.40", colour = "black", size=8)+ opts(title="NOx Buenos Aires Regresión")+opts(legend.position="none")+ xlab("NOx Pronosticado (ppb)")+ylab("NOx Original (ppb)")
> p11 <-ggplot(aire4)+aes(x=bai_nox_neu_pro, y=ba_nox_original)+geom_point(binwidth = 1, aes(size = factor(bai_dif_nox_neu),colour = factor(bai_dif_nox_neu)))+scale_colour_brewer(type="seq", palette=1)+annotate("text", x = 100,y = 300, label = "R2 = 0.33", colour = "black", size=8)+ opts(title="NOx Buenos Aires Neuronal")+opts(legend.position="none")+ xlab("NOx Pronosticado (ppb)")+ylab("NOx Original (ppb)")
> p9 <-ggplot(aire2)+aes(x=p_nox_reg_pro, y=p24nox)+geom_point(binwidth = 1, aes(size = factor(p_dif_nox_reg),colour = factor(p_dif_nox_reg)))+scale_colour_brewer(type="seq", palette=1)+annotate("text", x = 170,y = 400, label = "R2 = 0.73", colour = "black", size=8)+ opts(title="NOx Santiago Regresión")+opts(legend.position="none")+ xlab("NOx Pronosticado (ppb)")+ylab("NOx Original (ppb)")
> p12 <- ggplot(aire2)+aes(x=p_nox_neu_pro, y=p24nox)+geom_point(binwidth = 1, aes(size = factor(p_dif_nox_neu),colour = factor(p_dif_nox_neu)))+scale_colour_brewer(type="seq", palette=1)+annotate("text", x = 170,y = 400, label = "R2 = 0.80", colour = "black", size=8)+ opts(title="NOx Santiago Neuronal")+opts(legend.position="none")+ xlab("NOx Pronosticado (ppb)")+ylab("NOx Original (ppb)")
Y con esta fórmula sensacional de la receta de cocina para R, el Cookbook for R
> multiplot <- function(..., plotlist=NULL, cols) {
+ require(grid)
+
+ # Make a list from the ... arguments and plotlist
+ plots <- c(list(...), plotlist)
+
+ numPlots = length(plots)
+
+ # Make the panel
+ plotCols = cols # Number of columns of plots
+ plotRows = ceiling(numPlots/plotCols) # Number of rows needed, calculated from # of cols
+
+ # Set up the page
+ grid.newpage()
+ pushViewport(viewport(layout = grid.layout(plotRows, plotCols)))
+ vplayout <- function(x, y)
+ viewport(layout.pos.row = x, layout.pos.col = y)
+
+ # Make each plot, in the correct location
+ for (i in 1:numPlots) {
+ curRow = ceiling(i/plotCols)
+ curCol = (i-1) %% plotCols + 1
+ print(plots[[i]], vp = vplayout(curRow, curCol ))
+ }
+
+ }
obtenemos
multiplot(p7, p8, p9, p10, p11, p12,cols=3)
multiplot(p7, p8, p9, p10, p11, p12,cols=3)
R Development Core Team (2011). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL http://www.R-project.org/.
H. Wickham. ggplot2: elegant graphics for data analysis. Springer New York, 2009.
Teetor Paul. Cookbook for R. Pub. O'Really Media. 2011. 978-0-596-80915-7
