GGPLOT2     HEALTH     R     EMISSIONS     STATISTICS     GIS   WRF Open Air Air Pollution  

viernes, 24 de febrero de 2012

#ggplot2 #R Comparación de modelos de pronósticos de CO y NOx en Buenos Aires, Bogotá y Santiago, con regresiones y redes neuronales en R, ggplot2 y openair

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)

  • <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


# 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.
#
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))
    }
  }
}

obtenemos

multiplot(p7, p8, p9, p10, p11, p12,cols=3)



multiplot(p7, p8, p9, p10, p11, p12,cols=3)





Referencias:
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