Aprendiendo de los datos-Análisis Exploratorio-Part 1

El título de la entrada se refiere hacer con los datos algunas gráficas, de las cuales uno puede ir analizando posibles relaciones o visualizando conexiones entre ciertas variables. Esto se llama “Análisis Exploratorio”.

Muchas de las técnicas de Machine Learning tienen su contra parte gráfica; es decir, existe algún recurso gráfico para visualizar el resultado del uso de los algoritmos.

En la entrada trato de mostrar ejemplos de las herramientas gráficas estándar, tanto con R project, como con Python. De cierto modo trato de ir de ejemplos gráficos sencillos a herramientas gráficas más sofisticadas, hasta quizás acercarme a mostrar algunos ejemplos de visualización de datos.

Esto último es quizás sumamente atractivo como campo de trabajo y herramienta para trasmitir descubrimientos o hacer amigable ciertas relaciones entre variables analizadas, ya que permite diseñar algunas gráficas interactivas. Lo malo es que en general las mejores herramientas de visualización; son herramientas web y requieren explicaciones totalmente distintas. Pero mi recomendación es visitar la página de D3 y hacer algunos de los ejemplos para darse una idea general. Ahora hay muchas herramientas de visualización de datos, pero al final ciertas gráficas son las básicas y suelen ser las más útiles.

Análisis Exploratorio

Como todo visionario, muchos años antes de que el poder de las computadoras estuvieran en nuestra sala, teléfono y en el día a día, Jhon W. Tukey auguró desde antes de 1977 en su libro “Exploratory Data Analysis”, la relevancia de dedicarse a la investigación de diseño de herramientas gráficas para visualizar relaciones estadísticas. De cierto modo, dejando como primera etapa hacer una análisis exploratorio entre las variables analizar y posteriormente hacer una análisis confirmatorio.

Así que la relevancia que tiene hacer un previo análisis gráfico es fundamental, pero no siempre es determinante. Con esto quiero decir que si bien una gráfica facilita la visualización de relaciones entre variables, esto no lo confirma. Por lo cual es bueno posteriormente hacer pasar nuestra hipótesis por alguna técnica estadística que respalde nuestro descubrimiento.

Existe un libro sobre la construcción de gráficas en español del Dr. Juan Carlos Correa de la Universidad de Medellín [2], el cual está disponible en la red de manera gratuita. De este libro tomo los principios de William Playfair sobre el análisis gráfico:

Principios.

  • Los métodos gráficos son un modo de simplificar lo tedioso y lo complejo.
  • Los hombres ocupados necesitan algún tipo de ayuda visual.
  • Un gráfico es más  accesible que una tabla.
  • El método gráfico es concordante con los ojos.
  • El método gráfico ayuda al cerebro, ya que permite entender y memorizar mejor.

Teniendo en cuenta que los anteriores principios fueron pensados en 1800, nos debería de hacer pensar en la relevancia que tienen el trasmitir conocimiento por métodos gráficos. Por eso no es sorpresa que siempre contar con algunas gráficas en nuestros reportes o aplicaciones facilita la visualización y entendimiento de los fenómenos o hechos estudiados.

Análisis Exploratorio con R

En R Project, desde su creación se contó con herramientas gráficas diversas y al pasar de los años se han agregado herramientas o librerías sumamente útiles, como lattice, ggplot y ggvis y ahora interfaces gráficas como Shiny.

En esta entrada hago 3 ejemplos, el primero usando las gráficas nativas o base de R, el segundo haciendo uso la librería lattice y el último haciendo uso de la librería ggplot. Y muestro como este último ejemplo se sirve para visualizar de manera interactiva las gráficas haciendo uso de la librería ggvis. Para gráficas de una sola variable hago uso de las nativas y ggplot, para ver relaciones entre varias variables hago uso de lattice y ggplot.

Para cada uno de los ejemplos indico como extraer los datos, el primer ejemplo fue motivado por un ejemplo compartido por el Dr. Yanchang Zhao y el segundo ejemplo fue motivado por Institute for Digital Research and Education de la universidad de California y la referencia [3,4].

 Ejemplo con librerías nativas en R.

Los datos ha usar son Iris, los cuales se encuentran por default en R, los gráficos que explico en breve son para explorar la información de una sola variable.

Boxplot

#Ejemplo
#Cargo los datos para todos los ejemplos
data(iris)
#Tamaño de los datos
dim(iris)
#Tipo de variables
str(iris)
#Conocer las variables
names(iris)

#Boxplot
range(iris$Sepal.Length)
boxplot(iris$Sepal.Length, main="Boxplot de la Longitud de Sepal-Datos Iris", ylab="Tamaño", col="2")

La gráfica anterior los beneficios que tienen es que es fácil de calcular en el sentido computacional, por otro lado de manera rápida uno puede detectar si la distribución de los datos es simétrica o no, puede notar varios detalles sobre los cuantiles y la mediana. Para más detalles pueden verse de manera fácil en wikipedia, si subestiman este sitio como referencia, están subestimando mucha información de calidad. En este ejemplo se hace notar que no es simétrica la distribución y no se detectan valores atípicos.

Boxplot

Jitter-Gráficas de Puntos

Esta gráfica es fácil e informativa, cuando se tienen una buena cantidad de datos se puede ver los puntos de aglomeración, lo cual permite identifica como se comporta la variable analizada. En el ejemplo no muestra claramente algún comportamiento apreciable.

#Ejemplo Jitter
#Se usan los mismo datos de boxplot

stripchart(iris$Sepal.Length, method='jitter',vertical=TRUE,pch=1, col="2", main="Gráfica de Puntos", ylab="Rango de Datos")

#Método 2
stripchart(iris$Sepal.Length, method='overplot',vertical=TRUE,pch=1, col="3", main="Gráfica de puntos", ylab="Rango de Datos")

Gr´ficas_puntosGraficas_puntos2

Histograma y Densidad

Las dos gráficas que más noto que se usan para describir propiedades de alguna variable son los histogramas y las gráficas de densidad. La primera tienen varios detalles que es recomendable revisar, no los comento a detalle, pero hacer variar la cantidad de clases afecta a la construcción del gráfico, por otro lado hacer cambios en considerar la estimación de la frecuencia o de la frecuencia relativa también cambio no de forma , sino de escala. Así que esos detalles se pueden revisar en información de ayudan con el comando help() en la consola de R.

La idea central de estos dos gráficos es presentar una estimación de la densidad de la variable analizada, por eso la forma del histograma y la gráfica de la densidad tienen forma similar.

#Histograma
#Pongo varios ejemplos de como ir agregando parámetros a la gráfica

#Gráfica 1
hist(iris$Sepal.Length)

#Gráfica 2
hist(iris$Sepal.Length,col="2",ylab="Frecuencia", xlab="Longitud del Sépalo", main="Histográma de la Longitud de Sépalo-Datos Iris")

#Gráfica 3
hist(iris$Sepal.Length,col="2",ylab="Frecuencia", xlab="Longitud del Sépalo", main="Histográma de la Longitud de Sépalo-Datos Iris", labels = TRUE)

#Gráfica 4
hist(iris$Sepal.Length,col="2",ylab="Frecuencia", xlab="Longitud del Sépalo", main="Histográma de la Longitud de Sépalo-Datos Iris", labels = TRUE,border="7", nclass=15)

#Gráfica 5-Considerando las frecuencias relativas
hist(iris$Sepal.Length,col="2",ylab="Frecuencia relativa", xlab="Longitud del Sépalo", main="Histográma de la Longitud de Sépalo-Datos Iris", labels = TRUE,border="7", nclass=15, probability = TRUE)
hist(iris$Sepal.Length,col="2",ylab="Frecuencia relativa", xlab="Longitud del Sépalo", main="Histográma de la Longitud de Sépalo-Datos Iris", labels = TRUE,border="7", probability = TRUE)

Hist_1Hist_3
Hist_2

La gráfica de densidad hace una estimación de la densidad, para ello existen  varias funciones ya que la densidad no necesariamente es gausiana o de forma de campana. Lo que determina la aproximación es el tipo de “kernel” que se emplea para estimarla, en general el kernel que se considera como default es gaussiano, pero puede ser modificado. Para ver detalles de las opciones se puede consultar en la información de ayuda en R, poniendo el comando help(density).

#Gráfica de Densidad
#Si se pide estimar la densidad se obtienen lo siguiente
 
density(iris$Sepal.Length)

#Call:
# density.default(x = iris$Sepal.Length)

#Data: iris$Sepal.Length (150 obs.); Bandwidth 'bw' #= 0.2736

# x y 
# Min. :3.479 Min. :0.0001495 
# 1st Qu.:4.790 1st Qu.:0.0341599 
# Median :6.100 Median :0.1534105 
# Mean :6.100 Mean :0.1905934 
# 3rd Qu.:7.410 3rd Qu.:0.3792237 
# Max. :8.721 Max. :0.3968365

#Para la gráfica se hace uso de la funación plot()

plot(density(iris$Sepal.Length), xlab="Rango de Valores", ylab="Densidad", main="Gráfica de la Densidad", col="3",type="b", add=TRUE)

Density

Pie y Barras

Para estas dos gráficas se requiere tomar los datos para formar tablas, las cuales facilitan la gráfica tanto de barras como de pie. Principalmente se usan con variables categóricas, para estos datos la gráfica de barras no muestra gran funcionalidad por que la cantidad es la misma para las tres categorías.

#Se usan los mismos datos

table(iris$Species)
#setosa versicolor virginica 
# 50 50 50 
 
pie(table(iris$Species), radius = 0.9, col=rainbow(24),main="Gráfica Circular", clockwise=TRUE)
barplot(table(iris$Species), main="Gráfica de Barras",col="2", border="7", ylab="Cantidad")

Pie_1

Barras

La ventaja de usar gráficos para representar la información o el comportamiento de las variables analizadas resulta ser más atractivo cuando varios tipos de gráficas se combinan. En el caso de gráficos para una sola variable dos ejemplos rápidos se obtienen al combinar boxplot y jitter, histogramas y densidad.

#Combinación de gráficas

#Gráfica combinada 1
boxplot(iris$Sepal.Length, main="Boxplot de la Longitud de Sepal-Datos Iris", ylab="Tamaño", col="2")
stripchart(iris$Sepal.Length, method='jitter',vertical=TRUE,pch=1, main="Gráfica de Puntos", ylab="Rango de Datos", add=TRUE)

#Gráfica combinada 2
hist(iris$Sepal.Length,col="2",ylab="Frecuencia relativa", xlab="Longitud del Sépalo", main="Histográma de la Longitud de Sépalo-Datos Iris", labels = TRUE,border="7", probability = TRUE)
lines(density(iris$Sepal.Length), xlab="Rango de Valores", ylab="Densidad", main="Gráfica de la Densidad", col="3",type="b", add=TRUE)

Comb_1

Comb_2

En lo siguiente voy a replicar las gráficas haciendo uso de la librería ggplot2, a primera vista si no se conoce el modelo bajo el cual trabaja ggplot2 puede parecer bastante engorroso pasar de las gráficas base de R a las que proporciona esta librería.

Mi impresión es la siguiente, resulta muy fácil cuando uno aprende a usar R hacer el histograma y agregarle detalles a la gráfica con solo usar la función hist(), en cambio en ggplot2 se requiere hacer una combinación de comandos, ejemplo ggplot()+geom_hist()+…

Esto a primera vista hace parecer que ggplot2 complica mucho el hacer una simple gráfica, pero no es así. El valor o la relevancia de ggplot2 es cuando uno hacer gráficas más elaboradas, ya que cierta parte del procesamiento de los datos o de la detección de patrones de varias variables resulta inmediatamente visible, lo cual es hasta ahora para mi muy complicado hacerlo solo con las gráficas básicas de R. Así que invito a leer el libro de Hadley Wickham ( referencia [1]) o visitar el sitio de ggplot y por supuesto irse familiarizando con ” la gramática de gráficas”.

Boxplot

Como mencioné replico las gráficas en ggplot2, pero para eso hago en general dos versiones , haciendo uso de la función qplot y de ggplot. No explico más sobre las gráficas, pero comparto el código para replicar cada una.

#Boxplot
qplot(data=iris, y=iris$Sepal.Length,x=iris$Species, geom="boxplot", main="Boxplot por tipo de Iris",xlab="Tipo de Iris", ylab="Longitud de Sepal", colour="red")
ggplot(data=iris, aes(y=iris$Sepal.Length))+geom_boxplot(aes(iris$Species), col="red")+ylab("Longitud del Sepal")+xlab("Tipo de Iris")+ggtitle("Boxplot por tipo de Iris")

Boxplot_ggplot

Gráfica de Punto

#Gráficas de Puntos

qplot(data=iris, y=iris$Sepal.Length,x=iris$Species, geom="jitter", main="Jitter de la Longitud de la Sepal por tipo de Íris", ylab="Longitud", xlab="Tipo de Íris")
ggplot(data=iris, aes(y=iris$Sepal.Length,x=iris$Species ))+geom_jitter()+ylab("Longitud")+xlab("Tipo de Íris")+ggtitle("Jitter de longitud de Sepal")

Jitter_ggplot

Histograma

#Hitograma
qplot(data=iris, x=iris$Sepal.Length, geom="Histogram", col="red", main="Histograma", ylab="Cantidad", xlab="Rango de Valores")
qplot(data=iris, x=iris$Sepal.Length, main="Histograma", xlab="Rango de Valores", ylab="Frecuencia")+geom_histogram(col="red")

#Para generar el histograma con la medición de la densidad
ggplot(data=iris, aes(x=iris$Sepal.Length))+geom_histogram(aes(y=..density..), col="red")+xlab("Rango de Valores")+ylab("Frecuencia Relativa")+ggtitle("Histograma")

Hist_ggplot

Densidad

#Density
qplot(data=iris, x=iris$Sepal.Length, geom="Density",main="Densidad", ylab="Frecuencia Relativa", xlab="Rango de Valores")
#Para rellenar el área
ggplot(data=iris, aes(x=iris$Sepal.Length))+geom_density(colour="red", fill="orange")+ylab("Densidad")+xlab("Rango de Valores de la Longitud")+ggtitle("Densidad de la Longitud de la Sepalo")

Densidad_ggplot

Gráficas de Barras y de pie

#Barras

qplot(data=iris, factor(iris$Species), geom="bar", main="Ejemplo de gráfica de Barras", ylab="Cantidad",xlab="Tipo de Íris", fill=factor(iris$Species))

#Pie
ggplot(data=iris, aes(x=factor(iris$Species),fill=factor(iris$Species)))+geom_bar(width=1)+coord_polar()+xlab("Clasificado por Tipo de Especia")+ylab("")

Bar_ggplotPie_ggplot

Combinar gráficas en ggplot2, es muy sencillo solo se agrega a la gráfica que se está haciendo la otra gráfica que se desea. Ejemplo agregar la densidad al histograma resulta sencillo, ya que solo se necesita agregar  “+geom_density()”. Hago dos ejemplos para mostrar como funciona ggplot2.

#Boxplot y jitter

ggplot(data=iris, aes(y=iris$Sepal.Length))+geom_boxplot(aes(iris$Species), col="red")+geom_jitter(aes(iris$Species))+ylab("Longitud del Sepal")+xlab("Tipo de Iris")+ggtitle("Boxplot por tipo de Iris")

#Histogramas y Densidad

ggplot(data=iris, aes(x=iris$Sepal.Length))+geom_histogram(aes(y=..density..), col="red")+geom_density(colour="blue", size=1.5)+xlab("Rango de Valores")+ylab("Frecuencia Relativa")+ggtitle("Histograma")

Boxplot-jitter_ggplot

Hist_density_ggplot

Lattice…otra librería.

Cuando se requiere hacer una exploración de datos con varias variables, no resulta muy fácil pensar en como hacer eso y peor aún si se requiere hacer una comparación entre esas variables por varios años o varios periodos.

En general la mayoría de problemas involucran varias variables de distinta naturaleza, es decir; algunas pueden ser variables categóricas, indicadoras otras variables cuantitativas o cualitativas; en fin;  hacer un análisis de este tipo de datos requiere un cierto tiempo y las gráficas nativas a mi parecer hacer resultan no ser la mejor herramienta para trabajar.

En este caso las dos librerías, ggplot2 y lattice creo que ayudan bastante, detrás de estas dos librerías están dos teorías distintas sobre la exploración o visualización, es por eso que cada una requiere su tiempo de entrenamiento y sobre todo de investigación y aprendizaje. Dejo en las referencias los textos base para aprender a fondo como funcionan estas librerías.

Para las gráficas hago uso de un conjunto de datos que se descargarán automáticamente en R con una línea de código.

#Carga de datos y exploración básica

#Se cargan los datos desde el servidor
hsb2 <- read.table('http://www.ats.ucla.edu/stat/r/modules/hsb2.csv', header=T, sep=",")
#Se revisan los aspectos generales de los datos
head(hsb2)
dim(hsb2)
str(hsb2)

#Agrego una nueva variable
hsb2$gsex=factor(hsb2$female,labels=c("Male","Female"))
summary(factor(hsb2$gsex))
#Male Female 
# 91 109

 Boxplot y gráfica de puntos.

#Boxplot
bwplot(~read,hsb2, main="Boxplot de lectura", xlab="Lectura")
bwplot(~read|gsex,hsb2,main="Boxplot de lectura", xlab="Lectura")

#dotplot
dotplot(~read,hsb2,main="Gráfica de puntos de lectura", xlab="Lectura")
dotplot(~read|gsex,hsb2,main="Gráfica de puntos de lectura", xlab="Lectura")

Boxplot_latticeBoxplot_lattice2Dot_lattice2

Histogramas y Densidad

#Histograma
histogram(~math,hsb2,main="Histogramas por Género", xlab="Lectura",ylab="Porcentage Total")
histogram(~math|gsex,hsb2, ylab="Porcentage Total", xlab="Math", main="Histogramas por Género")

#Densidad
densityplot(~math, hsb2,main="Densidad de la variable math", xlab="Rango de valores", ylab="Densidad")
densityplot(~math|gsex, hsb2,main="Densidad de la variable math", xlab="Rango de valores", ylab="Densidad")

hist_latticeHist_lattice2Density_latticeDensity_lattice2

Cuantiles y Scatter Plot

Estas dos gráficas, no las había comentado. La primera se puede calcular para cualquier variable y es un modo rápido de identificar el comportamiento de la distribución de la variable, principalmente se puede validar el comportamiento de las colas de la distribución. La segunda es apropiada para detectar la posible relación entre dos variables y para explorar muchas variables se puede construir un panel de scatterplot para detectar visualmente entre qué variables es posibles que existe alguna relación.

#qqplot
qqmath(~math,hsb2, main="Gráfica de cuantiles o qqplot", xlab="qnorm",ylab="Variable Math")
qqmath(~math|gsex,hsb2, main="Gráfica de cuantiles o qqplot", xlab="qnorm",ylab="Variable Math")

#scatter plot
xyplot(write~read,hsb2, main="Scatter Plot de las variables write vs read",xlab="Variable read",ylab="Variable write")
xyplot(write~read|gsex,hsb2, main="Scatter Plot de las variables write vs read",xlab="Variable read",ylab="Variable write")

qqplot_latticeqqplot_lattice2Scatterplot_latticeScatterplot_lattice2

Gráficas en ggplot2

Lo único que hago en lo siguiente es replicar el tipo de gráficas separadas por una variable categórica.

Boxplot y dotplot

#boxplot
ggplot(data=hsb2,aes(y=read),colour=factor(hsb2$gsex))+geom_boxplot(aes(x=hsb2$gsex,fill=factor(hsb2$gsex)))+xlab("Genero")+ylab("Variable read")+
 ggtitle("Boxplot de la variable read vs genero")+theme(plot.title = element_text(lineheight=.8, face="bold"))

#dotplot
ggplot(data=hsb2,aes(y=read))+geom_jitter(aes(x=hsb2$gsex,colour=factor(hsb2$gsex)))+xlab("Genero")+ylab("Variable read")+
 ggtitle("Jitter de la variable read vs genero")+theme(plot.title = element_text(lineheight=.8, face="bold"))

Boxplot_ggplot_2Dotplot_ggplot

Histogramas y Densidad

#Histograma
ggplot(data=hsb2,aes(x=read))+geom_histogram(colour="black",aes(fill=factor(gsex)))+facet_grid(.~gsex)+xlab("Genero")+ylab("Variable read")+
 ggtitle("Histogramas de la variable read vs genero")+theme(plot.title = element_text(lineheight=.8, face="bold"))

#Density

ggplot(data=hsb2,aes(x=read))+geom_density(colour="black",aes(fill=factor(gsex)))+facet_grid(.~gsex)+xlab("Genero")+ylab("Variable read")+
 ggtitle("Densidad de la variable read vs genero")+theme(plot.title = element_text(lineheight=.8, face="bold"))

Hist_ggplot2Density_ggplot2

Cuantiles y Scatterplot

#qqnorm

ggplot(data=hsb2,aes(sample=read))+geom_point(stat="qq",aes(colour=factor(gsex)))+facet_grid(.~gsex)+xlab("Genero")+ylab("Variable read")+
 ggtitle("Cuantiles de la variable read vs genero")+theme(plot.title = element_text(lineheight=.8, face="bold"))

#Scatter plot

ggplot(data=hsb2,aes(x=read, y=write))+geom_point(aes(colour=factor(gsex)))+facet_grid(.~gsex)+xlab("Genero")+ylab("Variable read")+
 ggtitle("Scatterplot de la variable read vs write")+theme(plot.title = element_text(lineheight=.8, face="bold"))

Qqplot_ggplotScatterplot_ggplotLos ejemplos anteriores muestran la construcción de gráficos donde podemos ver la relación entre variables o explorar con mayor detalle. El siguiente paso entre el análisis exploratorio y lo que ahora es muy sencillo por las capacidades de computo es hacer análisis exploratorio con cierto nivel de interactividad.

Los casos más llamativos suelen ser las simulaciones, estas se pueden hacer tomando como parámetro una variable temporal. Ejemplo, si se desea ver el comportamiento de cierta población sobre la frontera de dos países se puede usar un mapa y ver con respecto al tiempo los flujos o movimientos de las personas.

Estas simulaciones a mi me parecen gratas para tener idea del fenómeno que se explora más que de la naturaleza de las variables.

En R se pueden hacer gráficos con cierta capacidad de interactividad está ggvis o shiny, así que los ejemplos siguientes los hago con ggvis. Para aprender a usar Shiny la recomendación es ver la documentación en la página oficial.

Página oficial Shiny

Gráficas de Puntos

#Puntos por Género
slider<-input_slider(1,100)

hsb2 %>% 
  ggvis(~as.factor(female),~read, 
        fill :="red"
       )%>% 
       layer_points(size :=slider)%>%
       layer_points(stroke := "black", fill := NA, size :=slider )%>%
       add_axis("x",title="Genero, 0 = Masculino y 1 = Femenino")%>%
       add_axis("y",title="Lectura")

La respuesta que se obtienen del código anterior es:

puntos_ggvis

#Densidad
hsb2 %>% 
  ggvis(~read, 
        fill :="red"
       )%>% 
       layer_densities(adjust = input_slider(0.1, 2))%>%
       add_axis("x",title="Densidad de la variable Lectura")%>%
       add_axis("y",title="")

La imagen que se obtiene son las siguientes:

Densidad_1.png

Variando la estimación se tiene:

Densidad_2.png

Gráfica de Barras

#Densidad
hsb2 %>% 
  ggvis(~read, 
        fill :="red"
       )%>% 
       layer_bars()%>%
       add_axis("x",title="Grafica de Barras")%>%
       add_axis("y",title="")

EjempBarras_ggvis.png

Más ejemplos de como funciona ggvis se pueden consultar en la página oficial:

Página Oficial de ggvis

Referencias:

1.- Ggpplot2 -Hadley  Wickham.

2.-Lattice – Deeppayan Sarkar.

 

 

 

 

Anuncios

Split-Apply-Combine- Parte 3

A terminar el ejemplo!!!

En la entrada Split-Apply-Combine parte 2, mostré unos ejemplos sencillos de como usar librerías de R y de Python para manipular los datos y obtener otro tipo de relaciones tales que de manera directa explorando la información sin transformarla no pueden ser detectadas dichas relaciones.

La relevancia de transformar la información y encontrar nuevas relaciones es parte fundamental para un buen análisis, esto requiere un poco de ensayo y error, de jugar con los datos y el tipo de datos, de buscar conexiones para tratar de descubrir alguna relación. Pero lo descubierto no implica que las “relaciones”  determinan algo relevantes, lo que en teoría uno debe de hacer es validar estos supuesto descubrimientos por medio de alguna técnicas estadística  para probar que la hipótesis o descubrimiento son por lo menos estadísticamente ciertos.

Lo que he compartido en las dos entradas previas(Parte 1 y Parte 2) es en general se puede considerar solo como exploración y diagnóstico, la idea de Split-Apply-Combine( SAC) es dividir , aplicar y combinar los resultados y por medio de esto obtener información que relevante.

El ejemplo en la referencia [1] da una idea de la metodología de manera fácil, por lo cual replico uno de los ejemplo para incentivar a revisar la referencia y observar como se aplica la metodología. El código que uso no es el mismo que comparten en dicho artículo, pero se obtienen los mismo resultados.

Las etapas son las siguientes:

  1. Se cargan los datos y se analiza la muestra de datos para un solo jugador( Split).
  2. Se construye un modelo el cual se aplica a cada uno de los jugadores(Apply).
  3. Se eligen dos elementos obtenidos del modelo para analizar el comportamiento de los datos en general y se muestran algunas gráficas para expresar los resultados (Combine).

Los datos corresponden a la base de datos baseball que se puede obtener desde la librería plyr.

#Split
#Se cargan las librerías
library(dplyr)
library(reshape)
library(plyr)
library(ggplot2)

#Se cargan los datos y se revisan que no cuente con NA
data("baseball")
head("baseball")
sum(is.na("baseball"))

#Se elige solo los datos correspondiente a baberuth

baberuth=filter(baseball,id=="ruthba01")
head(baberuth)
dim(baberuth)
#str(baberuth) se puede usar este comando para explorar que tipo de variables se tienen en los datos.

#Se agregan un par de variables para analizar lso años jugando y el número de carreras por bateos y el número de ocasiones que pasó a batear.

L=mutate(baberuth,cyear=year-min(year)+1,indice=rbi/ab)

#Se analiza la relación entre las dos nuevas variables

ggplot(L,aes(x=cyear,y=indice))+geom_line(col="red")+ggtitle('Relación entre datos')+
 xlab('Índice')+ylab('Ciclo de años')+geom_smooth(method='lm')

SAC_baseball

El código anterior lo que hace es elegir una muestra de datos y analizar la relación entre un par de variables que se construyen de las variables existentes, para ver la información de la descripción de la base se puede consultar por medio del comando ?—nombre de la base—.

La gráfica muestra una relación  y se hace un ajuste lineal como modelo para estudiar dicha relación. Esto no necesariamente es buen buen modelo, pero como ejemplo resulta fácil de entender y de aplicar.

#Apply
#Se construye el modelo 

model<-function(df){lm(rbi/ab~cyear,data=df)}

#Seleccionan los datos que permitan aplicar el modelo
L2=filter(baseball,ab>25)
#Se agregan las variables para toda la base

L3=mutate(L2,cyear=year-min(year)+1,indice=rbi/ab)
#Se aplica el modelo 
bmodels<-dlply(L3,.(id),model)
#Se observan el data.frame obtenido después de aplicar el modelo
head(bmodels)

Lo que se hace en el código anterior es lo siguiente, ya que se definen las nuevas variables y se analiza su relación se procede agregar esas variables para toda la base y se define una función que es el modelo aplicar a todos los datos, posteriormente se elige de la base aquellos elementos que permiten que el modelo se aplique. Si se aplica el modelo a toda la base sin revisar los valores de la variable ab, se obtendrá un error, esto debido al cálculo del modelo sobre los datos.

Se aplica el modelo por medio de la función dlplyr, el “dl” del comando indica que manda un data.frame para obtener una lista.Por último observo como son los datos por medio de head().

La parte de combinar los resultados implica usar algún estadístico para analizar los modelos, en la primera gráfica se observa que no es un buen modelo, entonces haciendo los ajustes lineales se puede extraer la r cuadrada que indica que tan bueno es el modelo de manera descriptiva. Ya que si la r cuadrada es muy alta indicaría que la relación es muy lineal y que el modelo es hasta cierto punto bueno, pero si es baja indica que no lo es.

Por lo tanto lo que hago, siguiendo al artículo o referencia [1], es extraer tanto la  r cuadrada de cada modelo, como la intercepto de la recta. Sabiendo un poco de geometría, se debe de tener en cuanta que la pendiente de la recta está relacionada con la correlación de las dos variables analizadas y por lo tanto con el valor de la r cuadrada.

#Combine

# Se define la función para extraer la R cuadrada de los modelos

rsq<-function(x)summary(x)$r.squared
#Se aplica el modelo a la lista y se regresa un data.frame
bcoef<-ldply(bmodels,function(x)c(coef(x),rcuadrada=rsq(x)))
head(bcoef)
#Se cambian nombres a las columnas
names(bcoef)[2:3]<-c("intercept","slope")
head(bcoef)

#Se revisa el histograma de los datos
ggplot(bcoef,aes(rcuadrada))+geom_histogram(col="red")+ggtitle('Histograma de R cuadradas')+
 xlab('Valores de R cuadrada')+ylab('Conteo')

#Se asigna el valor de las variables a toda la base
baseballcoef<-merge(baseball,bcoef,by="id")
head(baseballcoef)
dim(baseballcoef)

#Grafica 1
d1=qplot(data=baseballcoef,x=slope,y=intercept)+xlim(-0.05,0.05)+ylim(-2,2)
d1+geom_point(alpha=0.05,aes(size=rcuadrada))+ggtitle('Cruce de pendiente vs interceptos')+
 xlab('Pendientes')+ylab('Interceptos')+geom_hline(yintercept=0,col="red")+geom_vline(xintercept=0,col="red")

#Gráfica 2
d2=qplot(data=baseballcoef,x=slope,y=intercept)+xlim(-0.010,0.010)+ylim(-.25,.25)
d2+geom_point(alpha=0.05,aes(size=rcuadrada))+ggtitle('Cruce de pendiente vs interceptos')+
 xlab('Pendientes')+ylab('Interceptos')+geom_hline(yintercept=0,col="red")+geom_vline(xintercept=0,col="red")

Lo que hace todo el código anterior es extraer un estadístico que permita saber si es bueno el modelo, analizar gráficamente los datos. Por último analizar la posible relación gráfica entre los estadístico extraídos. Las gráficas son las siguientes:

Hist_baseball

SAC3SAC3_2

Entonces en resumen lo que se hizo, aun que fue replicar los cálculos del artículo trato de agregar algún modo distinto de obtenerlo con el código. La técnica al final es clara con este ejemplo, cada paso puede requerir más detalle y exploración sobre la información, pero en resumen la primera de estas últimas gráficas muestra que el modelo es malo, lo cual reta a buscar un mejor modelo para analizar los datos y otros estadísticos.

Ejemplo 2

Retomando los ejemplos de los datos de vuelos desde tres aeropuertos cercanos a la ciudad de NY, tenemos que que la base de datos es considerablemente más grande que el ejemplo anterior. Lo que trato de evitar con el ejemplo es hacer confuso como aplicar la técncias de SAC, así que lo que hago en un primer momento es algo similar al ejemplo anterior.

Tomo los datos, reviso el indicador “tailnum” para elegir los aviones que tienen el mismo número en la cola del avión y después analizo un par de variables que no son tan relevantes pero que permiten hacer un modelo sobre los tiempos de retraso al despegar y al arribar. Lo único que cambia es que la regresión o el modelo lineal es múltiple y es un buen modelo, lo cual no es de sorprender ya que si un avión se retrasa al despegar se espera que también se retrase al arribar.

Las etapas son las siguientes:

  1. Se cargan los datos y se analiza la muestra de datos para un solo registro “tailnum”, de hecho elijo el que aparece más ocasiones en los datos (Split).
  2. Se construye un modelo el cual se aplica a cada uno de los registros (Apply).
  3. Se eligen dos elementos obtenidos del modelo para analizar el comportamiento de los datos en general y se muestran algunas gráficas para expresar los resultados (Combine).
#Split

library(nycflights13)
library(plyr)
library(dplyr)
library(reshape2)
library(ggplot2)

data("flights")
dim(flights)
#[1] 336776 16

#missing value
sum(is.na(flights))
#[1] 60593
flights2=na.omit(flights)

head(flights2,10)
str(flights2)
dim(flights2)
#[1] 327346 16

#Seleccionamos el que tiene mayor número de apariciones

sort(summary(factor(flights2$tailnum)),decreasing = TRUE)

N725MQ=filter(flights2,tailnum=="N725MQ")

head(N725MQ)
dim(N725MQ)

#Agrego dos variables
N725MQ_2=mutate(N725MQ,dep_d=dep_delay-min(dep_delay)+1,arr_d=arr_delay-min(arr_delay)+1)

#Gráfica de cruce entre variables por mes
ggplot(data=N725MQ_2,aes(y=dep_d, x=arr_d))+geom_point(aes(colour=factor(month)))+geom_smooth(method='lm')

La gráfica que se obtiene es la siguiente:

fligths3_cruce

La gráfica muestra que existe una relación linea entre las dos variables que se definieron, pero esto además agrego una revisión de como se comportan los datos con respecto a los meses de vuelo y si bien no es claro una relación muestra que el comportamiento no es del todo aleatorio, en ciertos meses se presenta un retraso mayor que en otros. Un modo de ver esto es correr la regresión lineal solo para los datos de este registro de tailnum y revisar si es significativa, en este caso R cuadrada marca un 0.8 lo cual es bastante bueno.

Entonces lo que ahora hago es aplicar este modelo a los datos de toda la tabla.

#Apply
#Definición del modelo
modelo<-function(df){lm(dep_d~arr_d+month,data=df)}

#Procesamiento de toda la base
flight3=mutate(flights2,dep_d=dep_delay-min(dep_delay)+1,arr_d=arr_delay-min(arr_delay)+1)

#Aplicación del modelo
fmodelo<-dlply(flight3,.(tailnum),modelo)
head(fmodelo)
str(fmodelo)

#Extracción del estadístico
rsq<-function(x){summary(x)$r.squared}

fcoef<-ldply(fmodelo,function(x)c(coef(x),rcuadrada=rsq(x)))

#Extracción de datos del modelo
fcoef=ldply(fmodelo,function(x)c(coef(x),rcuadrada=rsq(x)))
head(fcoef)
names(fcoef)[2]=c("Intercepto")

#Histograma del modelo
ggplot(fcoef,aes(rcuadrada))+geom_histogram(col="red")+ggtitle('Histograma de R cuadradas')+
 xlab('Valores de R cuadrada')+ylab('Conteo')

La gráfica que se obtiene es la siguiente:

flights_hist

El histograma muestra que la regresión es un buen modelo para analizar todos los datos, y que se comporta bien la r cuadrada, es decir que es próxima a 1, lo único extraño es el conjunto de valores que tienen un valor pequeño  o cercano a cero, esto indica que debe de revisarse el para algunos registros. Esto de puede hacer de modo fácil usando la función filter().

 Lo último es combinar la información con la tabla original y analizar las relaciones entre los valores estimados del modelo. Estos datos permiten hacer otra exploración de la información y de cierto modo los datos fueron enriquecidos.

#Combine
#Unión de las dos tablas
flightscoef<-merge(flights2,fcoef,by="tailnum")
head(flightscoef)
dim(flightscoef)
str(flightscoef)

#Gráficas de resultados por aeropuerto de origen
ggplot(flightscoef,aes(rcuadrada))+geom_histogram(color="black",aes(fill=origin))+ggtitle('Histograma de R cuadradas')+
 xlab('Valores de R cuadrada')+ylab('Conteo')

#Gráficas de coeficientes

d1=qplot(data=flightscoef,x=month.y,y=Intercepto)
d1+geom_point(alpha=0.05,aes(size=rcuadrada))+ggtitle('Cruce de pendiente vs interceptos')+
 xlab('Pendientes de Meses')+ylab('Interceptos')+geom_hline(yintercept=0,col="red")+geom_vline(xintercept=0,col="red")

d2=qplot(data=flightscoef,x=arr_d,y=Intercepto)
d2+geom_point(alpha=0.05,aes(size=rcuadrada))+ggtitle('Cruce de pendiente vs interceptos')+
 xlab('Pendientes de Meses')+ylab('Interceptos')+geom_hline(yintercept=0,col="red")+geom_vline(xintercept=0,col="red")

Las gráficas que se obtienen son las siguientes:

flights_hist_origin

flights1_082015Flights2_082015

Lo que se observa es que efectivamente el comportamiento del modelo por aeropuerto de origen es similar, lo cual se ve en el histograma con diferentes colores por origen. Las otras dos gráficas muestran la primera la relación entre los coeficientes de la regresión, entre los cuales se aprecia cuales tienen cierto comportamiento, para profundizar en el análisis de estos coeficientes se requiere un poco más de detalle, pero en general la idea de aplicar SAC es también reflejada en este ejemplo.

Ultimo ejemplo.

En la entrada Split-Apply-Combine Parte 2 mostré que los datos podían ser analizados como series de tiempo después de una serie de transformación sobre la información, lo que trato de hacer con el ejemplo es aplicar un modelo de series de tiempo a los datos por medio de la metodología SAC.

Lo primero que hago es procesar la información y explorar solo los datos como series de tiempo de un solo registro de “tailnum”, esto se muestra gráficamente que tienen cierto comportamiento, el registro elegido es el que cuenta con mayor cantidad de datos o líneas en la base, el cual es “N725MQ”.

El código para explorar la información y los datos es el siguiente:

#Ejemplo de Series de tiempo
#Cargo las librerías requeridas
library(nycflights13)
library(plyr)
library(dplyr)
library(reshape2)
library(ggplot2)

#Revisión de los datos

data("flights")
dim(flights)
#[1] 336776 16

#missing value
sum(is.na(flights))
#[1] 60593
flights2=na.omit(flights)

head(flights2,10)
str(flights2)
dim(flights2)
#[1] 327346 16

#Seleccionamos el que tiene mayor número de apariciones
sort(summary(factor(flights2$tailnum,labels=levels)),decreasing = TRUE)

#Elijo el registro con mayor cantidad de líneas es "N725MQ"

N725MQ=filter(flights2,tailnum=="N725MQ")
head(N725MQ)
dim(N725MQ)

N725MQ_1=mutate(N725MQ,fechas=as.Date(paste(month,day,year,sep="/"),"%m/%d/%Y"))

Dat_1<-N725MQ_1%>%
 group_by(fechas,tailnum)%>%
 summarise(dis_m=mean(distance))

head(Dat_1)
qplot(data=Dat_1,x=fechas,y=dis_m,geom=c("line","point"))+geom_point(colour="red")+scale_x_date()+ggtitle('Seríe de N725MQ')

La gráfica que se obtiene es la siguiente:

SAC_timeseries

Estos datos se pueden tratar como una serie de tiempo, burdamente elijo un modelo para esta seria, el método que elijo es  ets (exponential smoothing), la intención es aplicar el modelo y tomar el estadístico RMSE, el cual la raíz de la media de los errores cuadrados y a grandes rasgos representa la desviación estándar entre los valores predichos y los estimados.

El código para aplicar el modelo a esta serie y obtener el RMSE, es el siguiente:

#Código
#Modelo
TS=ts(Dat_1$dis_m)
S=ets(TS)
S1=summary(S)
#El valor de RMSE es el siguente
S1[2]
201.2376
plot(TS, main="Serie de N725MQ para el año 2013",col="blue", ylab="Valores de la serie", xlab="Tiempo")
lines(S$fitted,col="2")

SerieN725MQ_1

Ahora el trabajo es elegir un cierto grupo de valores de tailnum que permitan aplicar el modelo y explorar como se comporta el estadístico y ver el comportamiento de los parámetros alpha y beta del modelo para cada un de los registros elegidos. Entonces se requiere que se transformen los datos, que se seleccionen cierto grupo de datos, se aplica un modelo y se extrae un estadístico para revisar el comportamiento del modelo.

#Transformando los datos
Datos1<-Vuelos1%>%
 group_by(fechas,tailnum)%>%
 summarise(dis_m=mean(distance))

Dato2=melt(Datos1,id=c("fechas","tailnum"))

Datos3=dcast(Datos1,fechas~tailnum)

Datos3[is.na(Datos3)] <- 0

#Eligiendo los tailnum con mayor actividad
Lista<-flights2%>%
 group_by(tailnum)%>%
 summarise(cantidad=n())

Aviones=filter(Lista,cantidad>365)

Lista_av=Aviones$tailnum

Lista_av2=c('fechas',Lista_av)

#Selección
Datos4=subset(Datos3,select=Lista_av2)

#Gráficas del comportamiento de algunos registros
ggplot(Datos4,aes(x=fechas,y=N324JB))+geom_line()
ggplot(Datos4,aes(x=fechas,y=N713MQ))+geom_line()
ggplot(Datos4,aes(x=fechas,y=N722MQ))+geom_line()

Como ejemplo gráfico algunos de los registros para ver su comportamiento en el año.

N324JBN713MQN722MQSe aprecia en los gráficos que el comportamiento es diferente y que es de esperar que el modelo no sea el mejor para todos los registros.

#Todos los datos
Datos5=melt(Datos4,id="fechas")
head(Datos5)
ggplot(Datos5,aes(x=fechas,y=value,color=variable,group=variable))+geom_line()+
 ggtitle('Todos los datos')+
 ylab('Media de las distancia de vuelos')+xlab('Fechas')

El comportamiento de todos los registros es el siguiente:

Todos_tailnum

Se hace notar que el comportamiento de separa en dos bloques, por lo cual es posible que para los vuelos que muestran mayor actividad el modelo sea más adecuado.

#Más transformaciones

#Ultima transformación
Datos6=dcast(Datos5,variable~fechas) 
head(Datos6)
colnames(Datos6)[1]<-'tailnum'
str(Datos6)

#Modelo
modelo<-function(df){ets(as.numeric(df))}

#Aplicación del modelo
Amodels<-dlply(Datos6,.(tailnum),modelo)

#Estadistico
rmse<-function(x){sqrt(x$mse)}

#Extracción del estadistico
Acoef<-ldply(Amodels,function(x)rmse(x))



ggplot(Acoef,aes(x=tailnum,y=V1))+geom_boxplot(col="red")+
 ggtitle("RMSE")+theme(plot.title=element_text(lineheight = 1,face='bold'))+ylab("RMSE")+xlab("TailNum")

Estas última parte del código hace el cierre de la metodología SAC, ya teniendo los datos preparados para aplicar el modelo se define la función que se va aplicar, se extrae el estadístico y se revisan los resultados con algún gráfico para tener identificado como se comportó el modelo.

Considerando que RMSE indica la desviación estándar de los errores del modelo, lo que se desea es que no sea muy grande ya que indica que los errores no están tan centrados. Esto es burdo pero es un modo sencillo de interpretar el gráfico que analizar como se comporta el modelo, el gráfico es el siguiente:

RMSE_ejemploRevisando en los datos, se observa que por ejemplo el vuelo N335AA tiene un comportamiento no muy favorable para el modelo,  haciendo la gráfica del comportamiento en el año se observa lo siguiente:

N335AAPara el cual se ve complicado que el modelo funcione bien, ya que siendo de suavizado exponencial no se favorable su aplicación.

 Espero que la entrada ilustre de manera general como se hace uso de la metodología SAC, si bien fueron tres entradas para explicar su aplicación la intensión es tener varios ejemplos tanto de gráficas exploratorias , como de gráficas no triviales y la aplicación de algún modelo sencillo que permita visualizar la metodología de manera clara.

Referencias:

1.-http://www.jstatsoft.org/v40/i01/paper

 

 

Split-Apply-Combine-Parte 2

Parte 2.-Procesando los datos y analizando otras posibles relaciones.

En la primera parte solo puse algunos ejemplos de gráficas que se pueden generar de manera directa de los datos sin hacer ningún procesamientos sobre la información.

Existen varias librerías para procesar datos en R project, las que aparecen en mayor cantidad de referencias son reshape2, dplyr y data.table. Para realizar un procesa SAC en los datos, no basta con solo revisar las variables de manera directa, es recomendable hacer un poco de procesamiento de los datos para realizar una exploración más sutil sobre la información.

En la entrada “Algo de procesamiento de datos” compartí algunos ejemplos sencillos del uso de las tres librerías de R y sobre el procesamiento de datos en Python con Pandas. Por lo cual no explico todos los detalles sobre las funciones, para más detalles se puede revisar esa entrada o las referencias que dejé en esa entrada.

Algo más

Usando dplyr y reshape2, se pueden explorar otro tipo de relaciones entre los datos. Ejemplo, la información de los vuelos cuenta con 16 variables de las cuales 3 de ellas indican el año, mes y día. Estas tres variables pueden permitir construir la variable de fecha la cual nos permite preguntarnos sobre la información desde un punto de vista de series de tiempo.

Una serie de tiempo es analizar el comportamiento de un indicador con respecto a la variación del tiempo, esto si bien no es claro lo muestro en el ejemplo. Formalmente es un proceso estocástico que tienen un parámetro t de tiempo.

La intención con ciertas transformaciones y agrupamientos de los datos es permitir explorar la información de manera diferente a lo que podemos hacer directamente desde la tabla de datos.

Tratamiento 1

#Se cargan los datos y las librerías requeridas para el análisis
#Se cargan las librerías
library(nycflights13)
library(dplyr)
library(reshape2)
library(ggplot2)

#Revisión de los datos
data("flights")
dim(flights)
#[1] 336776 16
#missing value
#Se eliminan las filas sin datos completos

sum(is.na(flights))
#[1] 60593
flights2=na.omit(flights)
head(flights2,10)
str(flights2)
dim(flights2)
#[1] 327346 16

#Agrupamos por fechas y aeropuerto de origen
a1<-flights2%>%
 group_by(year,month,day,origin)%>%
 summarise(count=n(),retrasos=sum(dep_delay+arr_delay),distancia=mean(distance))

#Se revisan el tipo de datos con los siguientes comenados
head(a1)
dim(a1)
#[1] 1095 7
str(a1)

#Agrego una variable para la fecha

a2=mutate(a1,fechas=as.Date(paste(month,day,year,sep="/"),"%m/%d/%Y"))
dim(a2)
#1095 8
str(a2)

head(a2)
#Source: local data frame [6 x 8]
#Groups: year, month, day

 #year month day origin count retrasos distancia fechas
#1 2013 1 1 EWR 300 11455 1039.8033 2013-01-01
#2 2013 1 1 JFK 295 5944 1297.1424 2013-01-01
#3 2013 1 1 LGA 236 2617 843.6695 2013-01-01
#4 2013 1 2 EWR 341 17242 1010.3900 2013-01-02
#5 2013 1 2 JFK 317 3641 1281.3502 2013-01-02
#6 2013 1 2 LGA 270 3589 841.9519 2013-01-02

Lo que se hizo con el código anterior fue construir desde los datos originales un data.frame o una tabla de datos con la información ordenada de otro modo y lista para poder explorarla como serie de tiempo. Lo que se hizo fue contabilizar la cantidad de vuelos por día, sumar los minutos de retraso al despegar y al arribar y calcular la distancia media recorrida por día.

Se pueden hacer algunas cosas con esta información, primero visualizar el comportamiento de los datos por cada aeropuerto de origen y después explorar si existe alguna relación entre las 3 series de tiempo.

 #Serie de tiempo 1
ggplot(a2,aes(x=fechas,y=count))+geom_line(aes(colour=factor(origin)))+scale_x_date()+facet_grid(.~origin)+geom_smooth(col="black")+
 ggtitle("Número de vuelos por año")+theme(plot.title=element_text(lineheight = 1,face='bold'))+ylab("Cantidad de vuelo")+xlab("Fecha de Vuelos")

La gráfica es la siguiente:

Serie_T1

#Serie de tiempo 2
ggplot(a2,aes(x=fechas,y=retrasos))+geom_line(aes(colour=factor(origin)))+scale_x_date()+facet_grid(.~origin)+geom_smooth(col="black")+
 ggtitle("Minutos de retrasos por vuelo en el año")+theme(plot.title=element_text(lineheight = 1,face='bold'))+ylab("Minutos de retraso")+xlab("Fecha de Vuelos")

La gráfica es la siguiente:

Serie_T2

#Serie de tiempo 3
ggplot(a2,aes(x=fechas,y=distancia))+geom_line(aes(colour=factor(origin)))+scale_x_date()+facet_grid(.~origin)+geom_smooth(col="black")+
 ggtitle("Media de la distancia de vuelo por año")+theme(plot.title=element_text(lineheight = 1,face='bold'))+ylab("Media de la distancia de vuelo")+xlab("Fecha de Vuelos")

La gráfica es la siguiente:

Serie_T3

Revisando el “cruce” entre la cantidad de vuelos, la distancia media y el total de minutos de retraso se obtienen las siguientes gráficas.

#Cruce 1
ggplot(a2,aes(count,retrasos))+geom_point(aes(colour=origin))+facet_grid(.~origin)+geom_smooth(col="black")
 ggtitle("Scatter plot de Cantidad de vuelos vs retrasos ")+theme(plot.title=element_text(lineheight = 1,face='bold'))+ylab("Cantidad de Vuelos")+xlab("Retrasos por Aeropuertos")

La gráfica que se obtiene es la siguiente:

SC1

#Cruce 2
ggplot(a2,aes(count, distancia))+geom_point(aes(colour=origin))+facet_grid(.~origin)+geom_smooth(col="black")+
 ggtitle("Scatter plot de cantidad de vuelos vs distancia media por año")+theme(plot.title=element_text(lineheight = 1,face='bold'))+ylab("Cantidad de Vuelos")+xlab("Distancia media por Aeropuertos")

La gráfica es la siguiente:

SP2

#Cruce 3
ggplot(a2,aes(y=distancia,x=retrasos))+geom_point(aes(colour=origin))+facet_grid(.~origin)+geom_smooth(col="black")+
 ggtitle("Scatter plot entre la distancia media y los retrasos por año")+theme(plot.title=element_text(lineheight = 1,face='bold'))+ylab("Media de la distancia de vuelo")+xlab("Retrasos por Aeropuertos")

La gráfica es la siguiente:

SC3

El único fin de las 6 gráficas anteriores es tratar de detectar alguna información o relación entre la información que con los datos originales no sea visible.

Analizando someramente las gráficas se puede observar en la primera gráfica que en los dos primeros aeropuertos el comportamiento es “similar”, ya que incremente la cantidad de vuelos de primavera hasta verano y posteriormente decrece cuando el año llega al invierno.La diferencia con el aeropuerto LGA que muestra un incremento conforme el año se acerca al invierno.

En la segunda gráfica se aprecia que para los 3 aeropuerto los retrasos incrementan en verano, con un poco de cuidado se aprecia que entre los meses de Febrero e inicio de Marzo se muestra un pico en los retrasos, el cual puede deberse a la temporada de cierre de nevadas.

El tercero resulta claro cual aeropuerto es requerido para vuelos con distancias mayores, JFK, y cual parece para vuelos cercanos. Revisando esto con la gráfica 1 se puede tener sospechas de que el incremento en la cantidad de vuelos en el aeropuerto LGA al cierre de año está relacionado con que se realizan vuelos cortos, lo cual puede estar relacionado con los festejos de cierre de año y periodos vacacionales.

Las últimas 3 gráficas son más sutiles, ya que se requeriría revisar la correlación y el coeficiente de información mutua para explorar si existe alguna relación interesante en los cruces. Por ejemplo la gráfica 5 muestra ciertas relaciones lineales negativa entre la cantidad de vuelos y la media de distancia en el aeropuerto LGA, pero en los otros dos no se muestra que exista dicha relación.

La gráfica 4 y 6 no muestran relaciones lineales, pero es posible que si tengan relaciones no lineales para lo cual se podría usar el índice de información mutual para tener algún indicio.

 Lo que debe de estar claro es que con unas agrupaciones y transformaciones se pudo explorar la información desde otra perspectiva y permite pensar en hacer cierto tipo de análisis, ejemplo las tres primeras gráficas permiten incitar en realizar un análisis de series de tiempo de los datos y con ellos crear un modelo predictivo de la cantidad de vuelos por mes o el comportamiento de los retrasos por mes.

Tratamiento 2

El siguiente ejemplo son dos partes, la primera es elegir el número de cola del avión con mayor cantidad de vuelos realizados y analizar sus comportamiento con respecto al año. Por último tomo un ejemplo que resulta similar al presentado en la página de referencia de la librería dplyr, con algunas ligeras variaciones.

Para elegir la cola de avión con mayor cantidad de vuelos, es decir; la variable tailnum agrupo la información con respecto a esa variable y resumo la información para conocer el registro con mayor número de vuelos.

#Agrupación
plane1<-flights2%>%
   group_by(tailnum,origin)%>%
   summarise(count=n(),retrasos=sum(arr_delay),med_distancia=mean(distance))

head(plane1)
#Source: local data frame [6 x 5]
#Groups: tailnum

 #tailnum origin count retrasos med_distancia
#1 D942DN JFK 1 2 944.0000
#2 D942DN LGA 3 124 824.6667
#3 N0EGMQ EWR 48 531 719.0000
#4 N0EGMQ JFK 27 361 517.6667
#5 N0EGMQ LGA 277 2622 688.2816
#6 N10156 EWR 144 1706 758.6458

max(plane1$count)
#536
filter(plane1,count==536)
#Source: local data frame [1 x 5]
#Groups: tailnum

#tailnum origin count retrasos med_distancia
#1 N725MQ LGA 536 2508 561.4515

De este modo tenemos que N725MQ es el registro que tienen mayor cantidad de vuelos. Entonces filtro la información original para analizar solo los datos correspondientes a este registro.

#Filtro
Plane_N725MQ=filter(flights2,tailnum=="N725MQ")
#Observo en breve la información
head(Plane_N725MQ)
#Source: local data frame [6 x 16]

#year month day dep_time dep_delay arr_time arr_delay carrier tailnum
#1 2013 1 1 832 -8 1006 -24 MQ N725MQ
#2 2013 1 1 1305 -10 1523 3 MQ N725MQ
#3 2013 1 1 1840 -5 2055 25 MQ N725MQ
#4 2013 1 2 1205 0 1349 4 MQ N725MQ
#5 2013 1 2 1805 -5 1946 1 MQ N725MQ
#6 2013 1 3 1131 -4 1314 -16 MQ N725MQ
#Variables not shown: flight (int), origin (chr), dest (chr), air_time (dbl),
# distance (dbl), hour (dbl), minute (dbl)

#Agrego dos variables nuevas
Plane_N725MQ1=mutate(Plane_N725MQ,velocidad=distance/air_time*60,fechas=as.Date(paste(month,day,year,sep="/"),"%m/%d/%Y"))

Con lo anterior puedo seleccionar las columnas que deseo analizar y por medio de la librería Reshape2 puedo cambiar su estructura para apreciar su comportamiento con respecto al tiempo, es decir; como serie de tiempo.

#Explorando el comportamiento de su velocidad media y la media de retrasos al arribar

#Se procesan los datos para construir otra estructura de datos 

B1<-Plane_N725MQ1%>%
 group_by(fechas,origin, dest)%>%
 summarise(count=n(),retrasos=mean(arr_delay),vel_media=mean(velocidad))

head(B1)
Source: local data frame [6 x 6]
Groups: fechas, origin

 fechas origin dest count retrasos vel_media
1 2013-01-01 LGA CRW 1 25 277.5000
2 2013-01-01 LGA DTW 1 3 295.2941
3 2013-01-01 LGA RDU 1 -24 335.8442
4 2013-01-02 LGA BNA 1 1 342.0896
5 2013-01-02 LGA RDU 1 4 300.6977
6 2013-01-03 LGA CLE 1 -16 318.2278

#Gráficas para explorar el comportamiento en el año

ggplot(B1,aes(x=fechas,y=vel_media,colour=dest))+geom_line()+scale_x_date()+geom_smooth(col="black")+
 ggtitle("Comportamiento de la velocidad en el año")+theme(plot.title=element_text(lineheight = 1,face='bold'))+ylab("Velocidad media de vuelo")+xlab("Fecha de Vuelos")

ggplot(B1,aes(x=fechas,y=retrasos,colour=dest))+geom_line()+scale_x_date()+geom_smooth(col="black")+
 ggtitle("Comportamiento de los retrasos en el año")+theme(plot.title=element_text(lineheight = 1,face='bold'))+ylab("Media de los retraos en los vuelos")+xlab("Fecha de Vuelos")

Las gráficas que se obtienen son las siguientes:

N725MQ

N725MQ_retrasos

En las dos gráficas se aprecian comportamiento peculiares para ciertos destinos de vuelo, cada sigue un comportamiento. En la primera gráficas resalta el comportamiento de BNA y CLE, el primero está sobre la línea de tendencia y el segundo tiene bajas de media de velocidad considerables para el periodo de Abril a Julio.

En la segunda gráfica se aprecia un comportamiento de mayor media en los retrasos entre los meses de Abril-Mayo y en otro momento en los meses de Junio-Julio.

En resumen, esta exploración permite centralizar el análisis en los destinos que resultan con comportamiento más extraños, con esta sospecha se puede procesar la información para analizar como se comportan el resto de vuelos con respecto a esos destinos.

Otro modo de procesar los datos  y realizar el análisis es revisar como se comportan los retrasos pero dividida por aeropuerto de origen de vuelo, esto se puede hacer fácil haciendo uso de la función melt de la librería Reshape2.

#Separación del análisis

#Se procesan lo datos originales
B2<-Plane_N725MQ1%>%
 group_by(fechas,origin, dest)%>%
 summarise(retrasos_arr=mean(arr_delay),retrasos_dep=mean(dep_delay))

#Construyo una table que permita el análisis de los datos de manera sencilla

datos=melt(B2,id=c("fechas","origin","dest"))

#Con el comando head() se visualiza como se ordena la información
head(datos)
#     fechas origin dest variable value
#1 2013-01-01 LGA CRW retrasos_arr 25
#2 2013-01-01 LGA DTW retrasos_arr 3
#3 2013-01-01 LGA RDU retrasos_arr -24
#4 2013-01-02 LGA BNA retrasos_arr 1
#5 2013-01-02 LGA RDU retrasos_arr 4
#6 2013-01-03 LGA CLE retrasos_arr -16

#Se analizan los datos gráficamente

ggplot(data=datos,aes(x=fechas,y=value,colour=dest,group=variable))+geom_line()+scale_x_date()+
 ggtitle("Comportamiento de los retrasos en el año")+theme(plot.title=element_text(lineheight = 1,face='bold'))+ylab("Media de los retraos en los vuelos")+xlab("Fecha de Vuelos")+
 facet_grid(.~origin)

La gráfica es la siguiente:

Retrasos_origen-destino

En este último ejemplo queda claro cual es el aeropuerto que se usa para despegar y el comportamiento de las medias de retrasos muestra ciertos picos, una ligera alza en los meses de Junio-Julio-Agosto. Este último gráfico se puede obtener desde la primera estructura de datos analizada, la única finalidad de hacer este último ejemplo era mostrar como melt() modifica la información no regresa algo parecido a una tabla pivot de Excel.

La última es prácticamente equivalente a la que se puede consultar en la página de referencia de dplyr, lo único que agregué fue mostrar la referencia de origen de vuelo.

#Última clasificación
#Se agrupan por 2 columnas
by_tailnum <- group_by(flights, tailnum,origin)
head(by_tailnum)
#Se construye el resumen de los datos
 
delay <- summarise(by_tailnum,
 count = n(),
 dist = mean(distance, na.rm = TRUE),
 delay = mean(arr_delay, na.rm = TRUE))

delay <- filter(delay, count > 20, dist < 2000)

#Se genera la gráfica para ver el cruce

ggplot(delay, aes(dist, delay)) +
 geom_point(aes(size = count,colour=origin), alpha = 1/2) +
 geom_smooth() +
 scale_size_area()+
 xlab('Media de la Distancia')+
 ylab('Media de los Retrasos')+
 ggtitle('Cruce de la media de las distancia vs media de los retrasos')+
 theme(plot.title=element_text(lineheight = 1,face='bold'))

La gráfica que se obtiene es la siguiente:

Cruce_dist_vs_retrasos

El patrón que se muestra posiblemente pueda ser estudiado con alguna técnica de clasificación la cual indicaría un aproximado de como se aglomera la información, pensando burdamente en usar cluster supongo que jugando con 5 cluster podría obtenerse un modelo simple sobre los datos.

En Python

PENDIENTE

Split-Apply-Combine-Parte 1

Sobre Split-Apply-Combine

La idea de esta entrada es aplicar la visión de dividir-aplicar-combinar (Split-Apply-Combine), en las referencias [1,2] se explican más detalles sobre el modelo o metodología guía, la cual fue inspirada por Map-Reduce. Esta última es una metodología apropiada cuando se tienen un cluster de computadoras, en este caso la versión Split-Apply-Combina se centra en el caso del análisis de datos mediante una sola máquina para diseñar una buena estrategia de análisis.

En la entrada desarrollo el ejemplo mediante R project y Python. La intención es mostrar como los dos lenguajes cuentan con las herramientas para realizar el mismo tipo de análisis. Trato de hacer un ejemplo amplio, por lo cual elegí una base más o menos grande de R project, la cual si bien no tiene muchas variables ( solo 16) permite hacer varias cosas con ellas. Por supuesto que el análisis puede tener otro tipo de enfoque, lo único que hago es usar esa base para emplear varias herramientas de R y de Python, las cuales en general uno las consulta de mera independiente.

Divido esta entrada en 3 partes.

Parte 1.-Exploración directa de los datos.

Parte 2.-Procesando los datos y analizando otras posibles relaciones.

Parte 3.-Definiendo técnicas para aplicar y concentrar esto para realizar un proceso Split-Apply-Combine.

Decidí dividir en tres partes para compartir algunos ejemplos sencillos sobre como hacer algunos gráficos y también un poco de análisis. Todos los gráficos usados entran perfectamente en los temas de exploración de datos o análisis exploratorio, en ocasiones solo uso algunas funciones para mostrar como funcionan y el valor que agrega a la exploración puede ser pobre y sin relevancia mayor.

En la referencia [3] se pueden ver detalles sobre las gráficas en R y en la referencia [4] sobre las gráficas en Python.

Parte 1.-Exploración directa de los datos.

En el código considero instaladas las librerías requeridas, en caso de contar con ellas  en caso de usar RStudio se puede hacer uso de Tools para instalar paquetes, lo cual es rápido y fácil.

Los datos corresponden a los tiempos de vuelo del año 2013, contienen información de los vuelos de 3 aeropuertos JFK, LGA y EWR.

#Código para cargar las librerías y datos
library(nycflights13)
library(ggplot2)
library(grid)
library(GGally)
library(gridExtra)#Funciones para la librería Grid

Los datos en general primero son explorados para conocer cuantos datos se tienen como missing values y las dimensiones de estos datos, es decir; cuantas columnas y cuantas filas tiene el archivo.

#Código para revisar las propiedades básicas de los datos
data(flights)
#Para conocer información de los datos 
#?flights

dim(flights)
#[1] 336776     16

#missing value
sum(is.na(flights))
#[1] 60593
flights2=na.omit(flights)
dim(flights2)
#[1] 327346     16

#tipo de datos
str(flights2)
#Aspectos básicos de los datos
head(flights2)
#También se puede ver la parte final de la tabla de datos
#tail(flights2)

#Para revisar el resumen estadístico de los datos
#Se usa el siguiente comando

summary(flights2)

Debido a que se tienen varias variables que son categóricas, se puede revisar su comportamiento considerándolas como “factores”. Un modo de revisar como se comportan ese tipo de datos se puede hacer con el siguiente comando.

#Explorando unas variables categóricas

sort(summary(factor(flights2$dest)),decreasing = TRUE)
sort(summary(factor(flights2$origin)),decreasing = TRUE)
#EWR    JFK     LGA
#117127 109079  101140

Cuando se revisan los datos de manera numérica uno no aprecia mucho la relevancia del comportamiento de los datos. Entonces siempre se buscan algunas gráficas sencillas que representen la información de manera clara y fácil.

Todas las gráficas siguientes pueden ser mejoradas, al final depende del interés de cada persona para definir el tipo de gráfica que se hace.

#Primeras Gráficas

p1=ggplot(data=flights2,aes(x=factor(origin),y=month))+geom_boxplot(col="black",alpha=0.7,aes(fill=factor(origin)))+
xlab('Origen de los Vuelos')+ylab('Mes de vuelo')+ggtitle('Comportamiento de los lugares de salida por mes')+theme(plot.title=element_text(lineheight = 2,face='bold'))

p2=ggplot(data=flights2,aes(x=factor(month)))+geom_bar(col="black",aes(fill=factor(origin)),position="dodge")+
  xlab('Mes de Vuelos')+ylab('Número de vuelos')+ggtitle('Comportamiento de vuelos por mes con respecto al aeropuerto de salida')+theme(plot.title=element_text(lineheight = 2,face='bold'))

grid.arrange(p1,p2,ncol=2,main="Gráficas del comportamiento de vuelos por mes")

La gráfica que se obtiene es la siguiente:Gráficas_vuelos_por_mesLa intensión es revisar la relación que existe entre la cantidad de vuelos por mes respecto al aeropuerto de salida u origen. La primera gráfica muestra como los 3 aeropuertos tienen un mes en el cual es la media de vuelos y por otro lado el JFK muestra que cubre un mayor número de meses con actividad. En las barras se observa que en EWR se realiza un mayor número de vuelos, además muestra las altas en vuelos en primavera y verano.

##Comportamiento de los tiempos de retraso de despegue y de arrive

q1=ggplot(data=flights2,aes(x=dep_time))+geom_histogram(color="black",aes(fill=factor(origin)))+
  xlab('Hora de despegue')+ylab('Número de Vuelos')+ggtitle('Hora de Despegue por aeropuerto de salida')+theme(plot.title=element_text(lineheight = 2,face='bold'))

q2=ggplot(data=flights2,aes(x=arr_time))+geom_histogram(colour="black",aes(fill=factor(origin)))+
  xlab('Hora de arrivo')+ylab('Número de vuelos')+ggtitle('Hora de Arrivo por aeropuerto de salida')+theme(plot.title=element_text(lineheight = 2,face='bold'))

grid.arrange(q1,q2,nrow=1,ncol=2,main="Gráficas del comportamiento de vuelos por mes")

El gráfico que se obtiene es:

Despegue_y_arrivo_por_aeropuertoSe muestran los histogramas de la cantidad de vuelos por hora de despegue y arribo divididos por el aeropuerto de origen. Se aprecian las horas en las cuales se realizan el mayor número de vuelos , los cuales son sumamente notorios. Por otro lado se aprecia la hora de llegada, las cuales muestran cuatro horas pico o el horario con mayor número de arribos.

En la base de datos se tienen registrados los minutos de retraso en despegues y en arribos a su destino, una idea sencilla es graficar el cruce de ambos datos y explorarlos.

#Gráfica de retrasos en despeque y arrivo
ggplot(data=flights2,aes(x=dep_delay,y=arr_delay))+geom_point(aes(colour=factor(origin)))+
  geom_smooth(method='lm',colour='black')+
  xlab('Retraso en Despegar')+ylab('Retraso en arrivo')+ggtitle('Relación entre retraso de despegue y arrivo por aeropuerto')+theme(plot.title=element_text(lineheight = 2,face='bold'))

La gráfica que se genera es la siguiente:

Relación_despegue_vs_arrivoAgregué la línea de tendencias entre los datos calculada con la regresión lineal entre los datos. Esta recta muestra una buena descripción gráficamente o eso aparenta, pero con un poco de cuidado se ve que los mejores datos descritos son aquellos que tienen mayor tiempo de retraso, mientras que la mayoría de datos se muestran como una mancha en los primeros 500 minutos.

Lo que se me ocurrió revisar el comportamiento de la distribución de los datos, haciendo uso de un histograma, después de observar que solo se apreciaba unas cuantas barras estimadas probé transformar los dato y calcular su logaritmo y luego volver a revisar su histograma. Lo que se observa de los datos transformados es la siguiente gráfica:

Hist_Dist_despegue_y_arrivosLa cual se puede obtener con el siguiente código:

#Tranformación de los datos e Histograma
l1=ggplot(data=flights2,aes(x=log(87+arr_delay)))+geom_histogram(color="black",aes(fill=factor(origin)))+
xlab('Log(87+arr_delay)')+ylab('Conteo')+ggtitle('Histograma de los retrasos en tiempo de arrivo')+theme(plot.title=element_text(lineheight = 2,face='bold'))

l2=ggplot(data=flights2,aes(x=log(43+dep_delay)))+geom_histogram(color="black",aes(fill=factor(origin)))+
  xlab('Log(43+dep_delay)')+ylab('Conteo')+ggtitle('Histograma de los retrasos de tiempo de despegues')+theme(plot.title=element_text(lineheight = 2,face='bold'))

grid.arrange(l1,l2,nrow=1,ncol=2,main="Gráficas del comportamiento de vuelos por mes")

Se observa entonces que los datos muestran comportamientos de distribuciones de cols pesadas, lo cual requiere una revisión más meticulosa. Haciendo el gráfico de cruce entre los dos tiempos de retrasos, pero considerando sus valores transformados se obtienen el siguiente gráfico log-log.

Log_vs_log_distancia-arrivosEl otro par de variables que sugiere que tienen alguna relación son la distancia del vuelo y el tiempo de vuelo. La gráfica se puede realizar con el siguiente código:

ggplot(data=flights2,aes(x=air_time,y=distance))+geom_point(aes(colour=factor(origin)))+
  xlab('Tiempo de vuelo')+ylab('Distancia recorrida')+ggtitle('Cruce entre la distancia y el tiempo de vuelo')+theme(plot.title=element_text(lineheight = 2,face='bold'))

La gráfica que se obtiene es la siguiente:

Cruce_Dis_vs_Tiempo_vueloLa gráfica a primera vista sugiere una relación lineal, pero también se puede explorar si no resulta favorable considerar el log-log de los datos, esto por aquellos datos que muestran desconexión en la parte superior de la gráfica.

ggplot(data=flights2,aes(x=log(distance),y=log(air_time)))+geom_point(aes(colour=factor(origin)))+
  xlab('Logaritmo del Tiempo de vuelo')+ylab('Logaritmo de la Distancia recorrida')+ggtitle('Cruce entre la distancia y el tiempo de vuelo')+theme(plot.title=element_text(lineheight = 2,face='bold'))

La gráfica que se obtiene es la siguiente:

Cruce_Log_de_Dist_vs_Timepo_VueloTambién muestra una relación lineal pero es posible que revisando un ajuste lineal entre los datos sea favorable considerar el log-log para cubrir una mayor  población de los datos.

Los siguientes gráficos solo son ejemplos para explorar la información en busca de alguna posible relación entre los datos.

JitterHoras_de_despegueEl código de los dos gráficos es el siguiente:

ggplot(data=flights2,aes(x=factor(carrier),y=factor(origin)))+geom_jitter()
#ggplot(data=flights2,aes(x=hour))+geom_bar(aes(colour=origin,fill=origin),position = "dodge")+coord_flip()

ggplot(data=flights2,aes(x=hour))+geom_bar(color="black",aes(fill=origin))+
  xlab('hora del día')+ylab('Conteo')+ggtitle('Histograma de las horas de despegue o salida')

Estos dos gráficos, sobre todo el primero; son informativos por ejemplo las barras negras muestran por aeropuerto qué aerolíneas tienen actividad. Es inmediatamente es claro cuales aerolíneas tienen presencia y la diferencia entre la actividad por aeropuerto.

La segunda gráfica muestra 3 momentos donde no se registran vuelos, lo cual resulta raro pero los datos de muestran esto, lo cual implica a ser explorado con mayor detalle.

Un gráfico que suele hacerse entre varias variables o entre una matriz de datos es el scatter plot, para explorar la posible relación entre cada par de variables. Para realizarlo en R con ggplot, se ejecuta el siguiente código:

#Scatter Plot
ggpairs(data=flights2[,c("dep_delay","arr_delay","air_time","distance")],columns=1:4,title="Datos")

La gráfica que se obtiene es la siguiente:

Scatterplot-flighLa gráfica de arriba solo muestra la relación entre varias variables de los datos no de todas. La información de la gráfica incluye la correlación, la cual muestra que es casi nula entre algunas de las variables y muy próxima a 1 en dos pares de variables. La otra opción es hacer una revisión entre las variables pero no solo usando la correlación, sino la el índice de información Mutua. Este último no solo considera la posible relación lineal, sino también una posible relación no lineal. Para realizar el cálculo en R project se puede usar la librería “entropy”.

Ejemplo en R la información mutua entre arr_delay y distance, es de 0.0406. Lo cual indica que no existe alguna relación que analizar.

Con la breve exploración anterior se puede ir pensando en el tipo de cosas que pueden resultar interesantes o que podían ser modeladas, ejemplo pronosticar conociendo el lugar de despegue la hora y la aerolínea el lugar de arribo, o construir un modelo para clasificar los datos de vuelos por el aeropuerto de origen. Pronosticar la cantidad de vuelos por aeropuerto para el siguiente año, o pronosticar los tiempo de atraso en arribar un vuelo conociendo el origen, el horario, el mes de vuelo y el lugar al cual arribará.

Por su puesto que el tipo de preguntar a responder depende de la persona que está explorando los datos o el negocio que requiere este análisis. Por lo cual en esta entrada solo deseaba mostrar algunos gráficos que se pueden ayudar a explorar los datos, al final la primera parte de realizar algún análisis es explorar aspectos estadísticos básicos y apoyarse de alguna gráfica ayuda a tener alguna visualización de la información.

Análisis similar con Python

 El procesamiento de datos en Python en general es realizado siempre en el módulo pandas, las estructuras de datos para procesar información es similar a R project. Se cuenta con estructuras de datos como Data.Frame y Series, como su nombre lo indica funcionan equivalentemente a las de R.

Se pueden usar en Python algunos de los datos de R project haciendo uso de algún módulo, ejemplo el módulo Stats Models tienen un api para conectarse con los datos de R project, lo cuales se encuentran su lista en Rdatasets.

En el ejemplo lo que hice fue exportar los datos de la librería nycflights13 a un archivo csv y cargarlo en una sección en Pandas de Python.

El código es el siguiente:

#Carga de datos
import pandas as pd
#Cargar datos
flights=pd.read_csv('flights.csv')

Para procesar o explorar la descripción básica de los datos se hacer con los siguientes código.

#Revisión de los datos
#Se ve un resumen de lso datos
flights.head()
flights.shape()

Unas de las cosas estándar al hacer análisis es procesar los missing data, lo cual se puede hacer en python con el siguiente código.

#Revisión
#Un modo burdo de ver que tienen missin values
flights[flights.isnull()!=True].count()
#Nuevo data.frame para los datos sin missin value
flights2=flights.dropna()
flights2[flights2.isnull()!=True].count()

Otro aspecto básico es revisar la descripción de la información de manera general, es decir; revisar aquella variables que son categóricas como origin o dest y también revisar la descripción estadística de variables que pueden ser consideradas continuas, como dep_delay o distance.

#Revisión de algunas variables
#Para conocer la información en general de los datos
flights2.info()
flights2.describe()
#Se exploran algunas de las variables que no son númericas y otras númericas
#Variables continua
flights2['dep_delay'].describe()
flights2['distance'].describe()

#Variables categoricas
flights2['origin'].astype('category').value_counts()
flights2['origin'].astype('category').describe()
flights2['dest'].astype('category').value_counts()
flights2['dest'].astype('category').describe()

Ahora para visar algunas gráficas básicas sobre las variables se usa en general matplotlib, pese a que pandas cuanta con algunas gráficas básicas como histogramas,gráficas de barras y boxplot.

Otra opción de reciente desarrollo es ggplot en python, la cual no cuenta con la misma flexibilidad que en R, pero está en desarrollo. La diferencia entre matplot y ggplot es en esencia que este último fue desarrollado con la perspectiva de “Grammar of Graphics”, mientras que matplot fue pensado en el marco de desarrollo de gráficos  que se usan en matlab. El concepto de ggplot, en recientes años a tomado relevancia ya que se pueden construir gráficos de manera rápida que involucran de fondo una transformación de los datos.

Un ejemplo es el siguiente, el cual se hace con los gráficos de base de pandas.

#Boxplot
df=flights2[['month','origin']]
df.boxplot(by='origin')

Boxplot_por_Mes

El gráfico es similar al que se obtienen en ggplot con R, pero en su defecto en R no se requiere seleccionar los datos y directamente indicar que se agrupen para hacer el gráfico.

Entonces para hacer un análisis gráfico rápido me resulta más fácil pensarlo en R que en Python, pero se obtienen ciertas ventajas en Python que no se obtienen en R.

Otra opción que existe es usar los gráficos de R desde pandas con el módulo rpy2, por medio del cual de cierto modo se usa todo lo que se quiera de R y se utiliza en python.

Un ejemplo de ggplot en Python, usando el módulo para python y no desde R, es el siguiente:

#Gráfica de barras
from ggplot import *


print ggplot(flights2,aes(x='month',fill='factor(origin)'))+\
 geom_bar(color="black")+\
 ggtitle('Barras por cada Mes')+\
 theme(plot_title=element_text(lineheight=2))+\
 ylab('Cantidad de Vuelos por Mes')+\
 xlab('Meses')

Barras_por_mes

El gráfico es parecido al primero que mostré en R. El boxplot de los vuelos en matplotlib se puede obtener del siguiente modo:

#Boxplot Maplotlib
import numpy as np
import matplotlib.pyplot as plt

plt.style.use('ggplot')

fig,ax1=plt.subplots()
ax1.boxplot(Mes,showmeans=True)
ax1.set_title('Boxplot de vuelos por mes')
ax1.set_xlabel('Valores')
ax1.set_ylabel('Meses')
plt.show()

Boxplot de los meses

Los histogramas para explorar el comportamiento de los tiempo de retraso y de despegue se obtienen en ggplot para python se obtienen con el siguiente código.

#Histogramas
from ggplot import *
print ggplot(flights2,aes(x='dep_time',fill='origin'))+\
 geom_histogram(color="black")+\
 ggtitle('Histograma de la hora despegue')+\
 theme(plot_title=element_text(lineheight=2))+\
 ylab('Cantidad de vuelos')+\
 xlab('Horario')

print ggplot(flights2,aes(x='arr_time',fill='origin'))+\
 geom_histogram(color="black")+\
 ggtitle('Histograma de la hora de arribo')+\
 theme(plot_title=element_text(lineheight=2))+\
 ylab('Cantidad de vuelos')+\
 xlab('Horario')

Histograma_despegues

Histograma_arribos

La gráfica del cruce entre los retrasos de arribo y despegue se obtienen de manera similar a la que se realiza en ggplot para R, lo único que se hace notar, es que la librería ggplot para python es más lenta al crear una gráfica con muchos datos.

#Scatter plot
from ggplot import *
print ggplot(flights2,aes(x='dep_delay',y='arr_delay'))+\
 geom_point()+\
 geom_smooth(color="red")+\
 ggtitle('Cruce de retrasos')+\
 theme(plot_title=element_text(lineheight=2))+\
 ylab('Retrasos en arribos')+\
 xlab('Retrasos en despegue')

Cruce_de_retrasos

Para hacer una gráfica similar en matplotlib se observa que el código es más elaborado, se obtiene una gráfica similar pero con mayor cantidad de líneas de programación.

#Ejemplo en Matplotlib
import matplotlib.pyplot as plt

for aero in ['LGA', 'JFK', 'EWR']:
 x =flights2.loc[flights2.loc[:,'origin']==aero,'dep_delay']
 y = flights2.loc[flights2.loc[:,'origin']==aero,'arr_delay'] 
 
 if aero=='EWR':
 color='red'
 plt.scatter(np.asarray(x.values)[:len(x)-1],np.asarray(y.values)[:len(y)-1],
 c=color, label=aero)
 
 if aero =='JFK':
 color='blue'
 plt.scatter(np.asarray(x.values)[:len(x)-1],np.asarray(y.values)[:len(y)-1],
 c=color, label=aero)

 if aero=='LGA':
 color='green'
 plt.scatter(np.asarray(x.values)[:len(x)-1],np.asarray(y.values)[:len(y)-1],
 c=color, label=aero) 
 

plt.legend(loc=4)
plt.grid(True)
plt.title('Scatter plot de los retrasos')
plt.ylabel('Retrasos en despegues')
plt.xlabel('Retrasos en arribos')
plt.show()

La gráfica es la siguiente:

Scatter_plot_retrasos

Otra gráfica fácil de hacer es la del comportamiento de los histogramas de los retrasos, tanto de despegues como de arribos.

El código es el siguiente:

#Histogramas
import matplotlib.pyplot as plt

plt.style.use('bmh')

x =flights2['dep_delay']
y = flights2['arr_delay'] 

x1=np.asarray(x.values)[:len(x)-1]
y1=np.asarray(y.values)[:len(y)-1]
plt.hist(x1, histtype="stepfilled",
 bins=20, alpha=0.8, normed=True)
plt.hist(y1, histtype="stepfilled",
 bins=20, alpha=0.8, normed=True)
plt.title('Histogramas de los retrasos de despeque y de llegada')
plt.xlabel('Valores de los retrasos')
plt.ylab('Cantidad de retrasos')
plt.show()

La gráfica es la siguiente:

Histogramas_de_los_retrasos

Para hacer un panel de gráficos como el de la exploración en R, se tiene que diseñar el panel con el número de gráficos que se desea. Esto en el fondo también sucede con R, solo que en matplotlib uno tiene que hacer el código con mayor detalle.

Espero que estas cuantas gráficas hechas en python de una imagen de lo versátil que es, se requiere en general procesar más los datos pero se puede tener mayor control sobre el tipo de gráficos que se pueden hacer.En la referencia [4] se pueden consultar muchas gráficas de muestra.

Referencias

1.-http://static.googleusercontent.com/media/research.google.com/es//archive/mapreduce-osdi04.pdf

2.-http://www.jstatsoft.org/v40/i01/paper

3.-http://ggplot2.org/

4.-http://matplotlib.org/

5.-http://pandas.pydata.org/pandas-docs/stable/groupby.html

Algo de procesamiento de datos

Algunas herramientas de manipulación de datos

En general la intención de esta categoría es  comentar sobre las técnicas de exploración, pero revisando y diseñando las entradas me di cuenta que dejé de lado una cosa que creo es fundamental, el procesamiento de datos o pre-procesamiento.

Creo que es central el tema, ya que con un poco de suerte uno tienen la información ordenada y mal que bien lista para poderla usar y explorar.

Ejemplo, ingresando a la siguiente liga yahoo uno puede extraer datos de yahoo fincance sobre el el indicador S&P, tiene varios formatos y la información está lista para analizarse, todo se encontrará bien ordenado. Pero no pasa lo mismo si uno quiere analizar algún texto, ejemplo; si uno ingresa la http://www.gutenberg.org puede descargar algún texto en formato HTML, el cual requiere procesarse para posteriormente poder analizarse con técnicas de Procesamiento de Textos o Text Mining, pero peor aún si uno construye un corpus que es un conjunto de textos, procesarlos requiere todo un labor.

Por otro lado, si se tienen la información ordenada lo que uno en ocasiones requiere es ver la relación entre algunas variables o estimar la media o conocer los máximos o mínimos. Para esto se requiere saber como manipular los datos.

 La exploración y construcción de algún modelo según sea el interés u objetivo, requiere de aprender o conocer los datos. Si bien para esto están los procesadores de bases de datos, lo que uno trata con herramientas como R o Python, es hacer uso de otro tipo de algoritmos que posiblemente requieren manipular bases considerablemente grandes antes de implementar el algoritmo.

Tratando de limitar la entrada comentaré sobre algunas librerías para procesamiento de datos en R y sobre el procesamiento de datos en Python. Por último agrego algunos comentarios sobre el procesamiento con Julia.

Librerías en R.

Es posible que existan más librerías que hacen manipulación de datos en R, pero las recomendadas y sobre todo reconocidas en varios cursos como en DataCamp, son reshape2, dplyr y data.table.

Comento sobre las tres y dejo los artículos donde pueden consultar información al respecto en la lista de  referencias.

Reshape2….después de reshape

La librería fue desarrollada por Hadley Wickham, la última versión se llama reshape2 la cual sufrió algunos cambios con respecto a la primera versión, como toda librería actualizada. La utilidad de la librería la trato de explicar con los ejemplos y estos son tomados del artículo [1].

Los datos para los ejemplos se cargan dentro de la misma librería. Este código lo corrí en RStudio.

#Cargamos las librerías para el ejemplo

library(reshape2)
library(ggplot2)

#En caso de que se quiera ver la información de la libreria se puede poner el comando
#library(help='reshape2')

#Primer ejemplo
data("smiths")

#Exploramos los datos

head(smiths)
dim(smiths)

#Se observa que solo cuenta con dos filas y 5 columnas
#La primera función importante "melt"
#para ver información de la función se puede usar
#?melt o help(melt)

#Se definen los id y las medidas

melt(smiths,id=c("subjects","time"),measured=c("age","weight","height")
#R regresa algo como esto:
# subject time variable value
#1 John Smith 1 age 33.00
#2 Mary Smith 1 age NA
#3 John Smith 1 weight 90.00
#4 Mary Smith 1 weight NA
#5 John Smith 1 height 1.87
#6 Mary Smith 1 height 1.54

#Obtenemos el mismo resultado de la siguiente forma

melt(smiths,id=1:2,measured=c("age","weight","height")

#Se hace la misma petición pero evitando los datos NA

melt(smths,id=1:2,measured=c("age","weight","height"),na.rm=TRUE)

# subject time variable value
#1 John Smith 1 age 33.00
#2 John Smith 1 weight 90.00
#3 John Smith 1 height 1.87
#4 Mary Smith 1 height 1.54

La función melt lo que hace es elegir de la pequeña tabla de datos las variables que considera id y el resto las considera como variables sobre las cuales muestra sus valores para cada id. Si se solicita la dimensión o tamaño de la tabla al usar melt se observa que hay un cambio en la estructura de la tabla. Es decir, en una petición por melt, pedir dim(melt(….)) deja ver que se cambio la forma de la tabla, de allí el origen del nombre.

La razón de su utilidad es que se pueden elegir varias variables que se desean explorar como una sola, en este caso subject y time y para cada una de sus posibles combinaciones se observan los valores de las variables age, weight y height.

Esto a primera vista no tienen relevancia o no se aprecia, pero con una tabla de datos poder modificar su estructura facilida un  rápido el análisis tanto gráfico como exploratorio. En el último ejemplo se usa una base de tamaño mayor. Pero antes comento sobre otra función.

#Función cast
Se observa a forma orginal del data.frame smiths

# subject time age weight height
#1 John Smith 1 33 90 1.87
#2 Mary Smith 1 NA NA 1.54

#Se construye una nuevo data.frame con la función melt

smithsm=melt(smiths,id=1:2,measured=C("age","weight","height"),na.rm = TRUE,value.name="Valores")

#Se observa el data.frame construido

smithsm 

#subject time variable value
#1 John Smith 1 age 33.00
#2 John Smith 1 weight 90.00
#3 John Smith 1 height 1.87
#4 Mary Smith 1 height 1.54

dcast(smithm,subject+time~variable)
#Se recupera el data.frame orgianal
 subject time age weight height
1 John Smith 1 33 90 1.87
2 Mary Smith 1 NA NA 1.54

Lo que se observa del código anterior es que se recupera el arreglo original, pero no solo eso hace la función dcast. Con algo de cuidado se observa que dentro de la función se ingresó una expresión Col1+Col2~Fila1+Fila2+….+Filan=variable.

Esta expresión indica la relación entre las columnas y las filas de los datos, que es lo que mostrará esta función. Se puede agregar un elemento más a la función dcast y hace parecer que melt y esta función funcionan de manera inversa. El siguiente ejemplo muestra como usarla para estimar cierto valores.

Para el ejemplo hice uso de una base que está en la librería reshape2, french_fries.

#Cargamos los datos french_fries 
data(french_fries)
ffm<-melt(french_fries,id=1:4,na.rm = TRUE)
#Se eligen solo los datos que no contienen NA
#Se explora la información
head(ffm)
dim(ffm)
#Se calcula con la función cast la media de cada variable, en este caso los valores de variable son:potato,buttery,grassy,rancid y painty

dcast(ffm,rep~variable,mean)
# Se obtiene una table como la siguiente:
# rep potato buttery grassy rancid painty
#1 1 6.956196 1.863663 0.5752161 3.916427 2.361383
#2 2 6.948851 1.784195 0.7528736 3.788218 2.682133

dcast(ffm,.~variable,mean)
# . potato buttery grassy rancid painty
#1 . 6.952518 1.823699 0.6641727 3.85223 2.521758

Lo que hace el código anterior es primero cargar los datos, después modificar sus estructura original para obtener una tabla con las primeras 4 variables como id, esto mediante la función melt. Usando la función dcast, se estima primero desde ffm la relación entre rep y las variables, se calcula su media. En la segunda petición de la función dcast se estima la media en general de las variables, sin separarse por rep como en el primer ejemplo.

El último ejemplo es equivalente al que se muestra en las referencias [1,2].

#Caso de estudio de la base de datos french_fries

data(french_fries)
dim(french_fries)
#[1] 696   9

#Se modifica la estructura de los datos por medio de melt
ffm<-melt(french_fries,id=1:4,na.rm=TRUE)
head(ffm)
# time treatment subject rep variable value
#1 1 1 3 1 potato 2.9
#2 1 1 3 2 potato 14.0
#3 1 1 10 1 potato 11.0
#4 1 1 10 2 potato 9.9
#5 1 1 15 1 potato 1.2
#6 1 1 15 2 potato 8.8

dim(ffm)
#[1] 3471    6

#Se revisan los valores de subject y de time
summary(ffm$time)
# 1 2 3 4 5 6 7 8 9 10 
#360 360 360 360 355 360 359 357 300 300 
summary(ffm$subject)
# 3 10 15 16 19 31 51 52 63 78 79 86 
#270 300 295 299 300 270 300 300 300 300 267 270

#Se hace un cruce entre subject y time para determinar analizar entre la pareja de datos entre los cuales no se tienen información

#Se agregan inicialmente los margenes y en el segundo comando no.
dcast(ffm,subject~time,length,margins=T)
dcast(ffm,subject~time,length)

#Se pueden analizar los valores máximos de cada una de las variables por medio de dcast

dcast(ffm,variable~.,max)

#Construye la gráfica de los valores de las variables
library(ggplot2)
qplot(value,data=ffm,geom="histogram",facets=.~variable,bindwidth=1,main="Histogramas de las variables")

La gráfica que se obtienen es la siguiente:

Histo_ejem_da

La idea general del ejemplo anterior es, construir un data.frame , explorar el comportamiento de los campos sin valor entre dos variables y revisar sus valores máximos. Por último construir una herramienta gráfica, en este caso un histograma; para revisar el comportamiento de las variables. En las referencias se explican más detalles sobre el ejemplo.

La librería reshape2 cuenta con más funciones, sobre las cuales no menciono nada, pero melt y cast son las más destacadas. Para conocer como usar el resto de funciones se puede revisar el manual de uso, referencia[3].

Librería dplyr…antes plyr

Esta librería también fue diseñada por Hadley Wickham, es fundamental debido a que los data.frame son usados como objetos y en general se usa para grandes tablas de datos.

La importancia de este paquete está ligada a la metodología de trabajo “split-apply-combine”, que fue planteada como una metodología similar a mapReduce que fue la solución de manipulación y procesamiento de grandes volúmenes de datos en Google, sobre este tema se puede revisar la referencia [4,5].

Para ser breve comparto un ejemplo publicado en parte en  [6] y trato de comentar de modo breve lo que se hace con esta librería.

Tiene alta relevancia con el procesamiento de datos, por lo cual invito a que se consulte en otros blog con algún tutorial amplio o en algunas de las referencias que dejo.

Cuando se manejan datos uno espera dos cosas, que sea rápido y que sea intuitivo, es decir; que sea fácil para uno pueda recordar los nombres de algunas funciones o la lógica al programar con los comandos. Para replicar el código se requiere una librería hflights que contiene solo una base considerablemente grande para mostrar como se usan las funciones principales de dplyr, o en su defecto se puede usar la librería nycflights13.  Este código lo corrí en RStudio.

#Se se cargan las librerías.

library(dplyr)
library(hflights)

#Cargamos los datos del ejemplo
#Se exploran el tipo de datos y las dimensiones de la tabla.
data(hflights)
head(hflights)
dim(hflights)
#[1] 227496 21

#Data frame local, mediante la función tbl_df()

flights<-tbl_df(hflights)
flights

#Para ver más registros de los dados por default 

print(flights,n=15)

#Convertir en un data frame normal y ver toda la información de manera normal en R

data.frame(head(flights))

###################################################
#Uso de la función filter()

#Método clásico en R
#Eligiendo vuelos de Enero 1

flights[flights$Month==1 & flights$DayofMonth==1,]

#Usando filter para elegir valores del dataframe

filter(flights,Month==1,DayofMonth==1)

#Usando filter y la pipe para condiciones "or"

filter(flights,UniqueCarrier=="AA"|UniqueCarrier=="UA")
#Uso del operador %in%

filter(flights,UniqueCarrier %in% c("AA","UA"))

###################################################
#Función slice
slice(flights,1:13)

###################################################
#Función arrange()
arrange(flights,Year,Month,DayofMonth)
arrange(flights,desc(FlightNum))

###################################################
#Selección de tres columnas
#Método clásico
flights[,c("DepTime","ArrTime","FlightNum")]
#Por medio de dplyr
select(flights,DepTime,ArrTime, FlightNum)

#Otras opciones para seleccinoar
select(flights, Year:DayofMonth,contains("Taxi"),contains("Delay"))

#Asignando títulos a los valores de algunas columnas
select(flights,Mes=Month,Vuelo=FlightNum)

El código anterior, si uno va revisando las salidas en R no resulta nada extraño, pero se pueden observar ciertas comodidades al seleccionar o al filtrar la información, que al hacerlo de manera normal en R resultaría más engorroso.

La siguiente parte del código se hace uso de comandos u operadores propios de la librería.

#Usos de operadores
#Anidación

#Método 1-por filtro
filter(select(flights, UniqueCarrier, DepDelay), DepDelay > 60)

#Método 2-por cadena
flights%>%
 select(UniqueCarrier, DepDelay)%>% 
 filter(DepDelay > 60)

#Distancia Euclideana entre dos vectores

x1<-1:5;x2<-2:6
sqrt(sum((x1-x2)^2))

#Método de cadena para para la distancia Euclideana
(x1-x2)^2%>%sum()%>%sqrt()

##################################################
#Reordernar filas-arrange
#Método estándas en R

flights[order(flights$DepDelay),c("UniqueCarrier","DepDelay")]

#Método en dplyr

flights%>%
 select(UniqueCarrier,DepDelay)%>%
 arrange(DepDelay)

#Método en dplyr en orden descendente

flights%>%
 select(UniqueCarrier,DepDelay)%>%
 arrange(desc(DepDelay))

####################################################
#Mutate
#Creando nuevas variables en los datos
#Método clásico por R

flights$Speed<-flights$Distance/flights$AirTime*60
#Visualización de la nueva variable
flights[,c("Distance","AirTime","Speed")]

#Método 2-Aproximación por dplyr

flights%>%
 select(Distance,AirTime)%>%
 mutate(Speed=Distance/AirTime*60)

flights<-flights%>%mutate(speed=Distance/AirTime*60)

####################################################
#summarise: Reduce variables a valores

#Método estándar en R

head(with(flights, tapply(ArrDelay, Dest, mean, na.rm=TRUE)))
head(aggregate(ArrDelay ~ Dest, flights, mean))

#Método por dplyr

flights %>%
 group_by(Dest) %>%
 summarise(avg_delay = mean(ArrDelay, na.rm=TRUE))

#Se calcula el porcentaje de vuelos cancelados o desviados para cada aerolínea

flights %>%
 group_by(UniqueCarrier) %>%
 summarise_each(funs(mean), Cancelled, Diverted)

Se observa que los operadores permiten tener mejor legibilidad sobre las operaciones pedidas en el data.frame. Por otro lado al contar con nombres intuitivos permite que sea claro al leer el código.

Por último, solo agrego un ejemplo de la referencia [6] de construir subconjuntos y graficar el comportamiento medio de retrasos.

#Procesamiento y gráfica de datos procesados
#Cargamos las librerías necesarias

library(ggplot2)

by_tailnum <- group_by(flights, TailNum)
delay <- summarise(by_tailnum,
 count = n(),
 dist = mean(Distance, na.rm = TRUE),
 delay = mean(ArrDelay, na.rm = TRUE))
delay <- filter(delay, count > 20, dist < 2000)

#Se visualiza el data.frame

delay
dim(delay)

#Distancia media por vuelo

ggplot(delay, aes(dist,delay)) +
 geom_point(aes(size = count), alpha = 1/2,colour='purple') +
 geom_smooth() +ggtitle("Variación de los retrasos contra distancia")+xlab("Distancia")+ylab("Retraso")


La gráfica que se obtiene es la siguiente:

Retrasos_vuelos

 

A primera vista no parece del todo relevante la librería, pero algunas de las cosas que no comento es la conexión con bases de datos, esto permite hacer un análisis con herramientas más ricas de estadística y mejores gráficas propias de R, que las que cuenta un gestos de bases de datos.

Al final, cuando se hace un análisis de datos dos piezas son importantes, el enfoque de la investigación que se realice y las herramientas con las cuales se haga. En este caso, explotar el uso de esta librería dependerá del proyecto que se aborde.

Nota: Para procesar algunas tablas también se puede hacer mediante la librería plyr, la cual es un poco más lenta para tablas grandes pero es igual de amigable que dplyr.

Librería data.table

La librería es de años recientes y poco a poco gana popularidad, se ha mejorado para que sea aún más rápida. Prácticamente sus uso es para procesar datos en R. En buena medida la sintaxis es similar a SQL.

En un base como las que se encuentran en las librerías de R no se notará su utilidad, pero en particular para bases grandes permite hacer manipulación de los datos. Existen varias comparaciones en cuanto al rendimiento contra tiempo al comparar  operaciones similares entre dplyr y data.table, la ventaja de la primera es la facilidad de operar con sus funciones y la segunda resulta ser más eficiente ante volúmenes considerablemente muy grandes.

Solo comento algo breve sobre la creación de tablas en esta librería, la definición de key’s y un poco sobre la agrupación. Se puede ver más detalles en las referencias [7,8]

#Código de data.table
#Se carga la librería

library(data.table)

#Se compara la creación de un data.frame y una tabla 
DF<-data.frame(x=c("b","b","b","a","a"),v=rnorm(5))
DF

class(DF)

#Generado con data.frame de data.table

DT<-data.table(x=c("b","b","b","a","a"),v=rnorm(5))

DT
class(DT)

#Importan datos de  R project

data(cars)
head(cars)

#Se transforman a una data.table

CARS=data.table(cars)
CARS

#La observación es sobre el indicador MB
#Muestra las tablas creadas con data.table
tables()

#Se se pueden aplicar las funciones estándar sobre el data.table

sapply(DT,class)


El código anterior muestra solo como crear tablas con data.table, al observarlas resultan prácticamente iguales a un data.frame estándar y heredas el tipo de funcionalidades o aplicaciones de funciones comunes en R, como sapply.

#Key en data.table

#Se usa DT del código pasada

setkey(DT,x)
DT

#Se observa en tables() que ahora key tienen "x" respecto a DT

#Ejemplo para seleccionar filas
DT[2,]
DT[x=='b',]
#Despues de definir la key se tienen
DT["b",]

Las key puedes ser definidas por varias columnas, al observar las tablas definidas en data.table resultan tener comportamientos similares a los data.frame normales.

Un modo de comparar los procesos normales de R con respecto a los procesos de data.table es comparando los tiempos requeridos para seleccionar partes de una tabla de tamaño, la tabla contendrá 10 millones de filas y 3 columnas.

El siguiente código son solo 3 ejemplos comparativos del tiempo que requiere R y del tiempo que se requiere en data.table para seleccionar o procesar la misma muestra de datos.

#Ejemplos comparativos

#Comparación de busqueda vector scan y binary search

grpsize = ceiling(1e7/26^2)
#Valor
grpsize

#Definición y tiempo de creación d eun data.frame
tt=system.time( DF <- data.frame(
 x=rep(LETTERS,each=26*grpsize),
 y=rep(letters,each=grpsize),
 v=runif(grpsize*26^2),
 stringsAsFactors=FALSE))
#Se revisan las propiedades generales del data.frame
head(DF)
tail(DF)
dim(DF)

tt=system.time(ans1 <- DF[DF$x=="R" & DF$y=="h",]) # 'vector scan'
tt
#Se revisan los datos obtenidos de la selección de valores
head(ans1)
tail(ans1)
dim(ans1)

#Se crea una data.table
DT=as.data.table(DF)
system.time(setkey(DT,x,y))

ss=system.time(ans2 <- DT[list("R","h")])
ss
head(ans2)
tail(ans2)
dim(ans2)

############################
#Otro ejemplo de velocidad entre R y data.table

#Por data.table
system.time(ans1 <- DT[x=="R" & y=="h",])

#Por R 
system.time(ans2 <- DF[DF$x=="R" & DF$y=="h",])

#Validamos que son la misma tabla
mapply(identical,ans1,ans2)

#########################################
#Ejemplo de aplicación de una función sobre las columnas

#Por medio de R
ttt=system.time(tt <- tapply(DT$v,DT$x,sum)); ttt

head(tt)

#Por medio de data.table
sss=system.time(ss <- DT[,sum(v),by=x]); sss

head(ss)

Los resultados que se obtendrán al replicar al código dan una idea de las mejoras que se obtienen al manipular datos con data.table, pero hay recientes estudios donde se compara con dplyr y pandas( de Python). Pero aún así se está mejorando parte de los algoritmos para que sea aún más eficiente. Otra opción es revisar los ejemplos que cuenta la librería , basta con poner en R el comando example(data.table).

Un ejemplo de las ventajas se puede ver en la siguiente liga:

http://stackoverflow.com/questions/11054208/lapply-and-do-call-running-very-slowly

Espero las 3 librerías den una idea en general de lo que se puede hacer en R para procesar datos, unos cursos breves pero que dan una idea global se encuentran en DataCamp, de los cuales la plataforma permite realizar algunos de manera libre y otros tienen costo.

Pandas-Python

La librería para realizar análisis de datos en Python es Pandas, la cual en casi todas las distribuciones de paquetes de librerías o módulos, como Anaconda; viene por default.

Existe mucha referencia respecto al manejo de dicho módulo. Los objetos principales son Series y DataFrame, los cuales son nombres similares a los usados en R project. 

El módulo fue creado por Wes McKinney y la documentación, casos de estudio y todo lo relacionado con pandas se puede consultar en la referencia[9], la cual es sumamente amplia y cuanta con ligas a otras referencias de manuales breves, ejemplos de uso y actualizaciones. El texto clásico y base para aprender sobre este módulo fue escrito por el mismo Wes, es la referencia [10]. Por supuesto es altamente recomendable y los datos para realizar los ejemplos, se encuentran en el repositorio:

https://github.com/pydata/pydata-book

Para el uso del código siguiente considero que se tiene instalado Python 2.7.9 y de preferencia con el paquete de módulos de Anaconda.

El ejemplo para verificar que funciona adecuadamente la el conjunto de módulos con panda, se corre el siguiente código en ipython.

#Prueba 
%pylab
plot(arange(10))

 El código debe de arrojar la siguiente imagen.

Pandas_1

Lo que se hizo fue probar que %pylab, las librerías básicas de desarrollo en ipython se encuentran funcionando. Los módulos base son numpy, matplotlib y pandas.

En los ejemplos solo muestra el uso de pandas para manipulación de datos, con la finalidad de mostrar ejemplos similares a los de R project. En el código de ejemplo considero en general la estructura de datos DataFrame y solo al final doy un breve ejemplo del uso de la estructura Series.

Creación de DataFrame y manipulación básicas.

#Uso de Pandas
#Considerando que se inicio ipython 

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

from pandas import DataFrame

#Se definen datos para creas un DataFrame

data={'Estado':['Jalisco','Sonora','Guanajuato',Veracruz'],'Población':['7838010','2892464','5769524','7985893'],'clave':['09','13','14','23']}

Estados=DataFrame(data)

En Ipython se obtendrá una tabla de datos, sobre la cual se puede manipular la columnas, las filas, agregar nuevas columnas y filas. Si se desea conocer aspectos básicos o métodos que  la estructura de datos DataFrame tienen por default, se puede hacer lo siguiente para conocerlos:

#Datos
Estados.<tabs>

Se presentará una lista de todos los métodos disponibles, ejemplo:

#Métodos
Estado.shape
#Dará el número de columnas y filas que tienen el DataFrame

Estado.head()
#Mostrará solo la cabecera, es decir; una muestra de los primeros 4 registros o filas del DataFrame

Estado.describe()
#En caso de tener una tabla con muchos datos, da un resumen estadístico de la información 

#Si se desea sabes para que sirve cierto método se hace uso de los siguiente.
help(Estados.*)
#Donde en lugar de * se pone el método que se desea conocer

Varios aspectos no menciono, pero en caso de no estar familiarizado con estructuras de datos, recomiendo la referencia [11]. Lo que se oculta al construir el DataFrame por este método es que se asigna data={}, estos paréntesis indican que es un diccionario y [] indica que es una lista.

En R project se hace uso de ciertas librerías para manipular los datos y poder procesarlos de mejor manera, como reshape o dplyr, en python todas esa operaciones se hacen solo con pandas.

Más aún, una postura ante el análisis de datos es la metodología “Split-Apply-Combine”, la cual fue descrita y ejemplificada por Hadley Wickham en la referencia [12]. Esta metodología fue considerada por Wes para el desarrollo de algunos aspectos de pandas, sobre todo en este tema comparto una entrada en esta categoría.

Ejemplo es sencillo y solo muestro algunas operaciones o procesos simples, para mayor detalle y sobre todo para revisar ejemplos de como hacer usos de pandas se puede consultar en la referencia [9].

En el siguiente código se extraen datos desde yahoo finance de los indicadores AAPL e IBM. Se obtiene  un DataFrame de datos, los cuales para mostrar algunas cosas básicas se agregan algunas columnas, de piden algunos datos básicos para explorar la información y se utilizan los dos DataFrame solo para mostrar una operación simple de “Muchos a Muchos”.Se hacen algunas operaciones para agrupar y transformar los datos y por ultimo una gráfica simple sobre el comportamiento de los cierres de los dos indicadores.

#Ejemplo Pandas
#Se cargan las librerias necesarias

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from pandas import DataFrame
import pandas.io.data
import datetime

#Se cargan datos desde yahoo para AAPL y IBM
aapl=pd.io.data.get_data_yahoo('AAPL',start=datetime.datetime(2010,1,1),end=datetime.datetme(2015,1,1))
ibm=pd.io.data.get_data_yahoo('IBM',start=datetime.datetime(2010,1,1),end=datetime.datetme(2015,1,1))
 
#Se exploran los datos, aspectos básicos solo para aapl

type(aapl)
aapl.describe()
aapl.columns
aapl.shape
aapl.head()

#Las dimensiones de los datos

aapl.shape
ibm.shape

Un podo de Missing Values.

#Ejemplo Pandas
#Para trabajar con los Missin Value se toman datos de aapl

aapl_Close=aapl.Close
type(aapl_Close)
#Se oservará que es una estructura de datos serie
#Se toman 200 enteros de entre 0 y 1258 como ejemplo de Missing Value

ind=np.random.randint(1258,size=(200))
aapl_Close[ind]=np.nan

#Para identificar que se cuentan con Missin Value
aapl_Close.isnull()
aapl_Close[aapl_Close.isnull()]
aapl_Close[aapl_Close.notnull()]

#Para elegir solo los datos sin nan
aapl_Close.dropna()
#Para asigar algun valor a los nan, ejemplo se asigna el valor cero "0"
aapl_Close.fillna(0)

Algunas operaciones como la librería Reshape en R project.

#Ejemplo Pandas
#Se define otra columna para elegir el indice como columna

aapl['ind']=aapl.index 

#Se define un nuevo DataFrame por medio de melt

df=pn.melt(aapl,id_vars=['ind'])
df.head()
#Se modifica la estructura de los datos
#Como ejemplo usa df para construir una tabla pivot

pivoted=df.pivot('ind','variables','value')

pivot.head()
#Tienen la misma estructura de los datos originales appl

Ejemplos simple de muchos a muchos y agrupación de los datos.

#Ejemplo Pandas
ibm['index']=ibm.index 

df1=pd.merge(aapl,ibm,left_on='ind',right_on='index',suffixes=('_aapl','_ibm'))

df1.head()
df1.shape
#Se obtienen una tabla con todos los datos de los dos indicadores

#####group by######

#df que proviene de melt
df2=df.groupby('variable')

#Se genera un DataFrame con estructura distinta al origina, el cual muestra agrupamientos por los valores de "variable"
df2.head()

#Algunas funciones sobre el DataFrame agrupado

df2.agg([np.max,np.min,np.mean,np.std])

#Ejemplo de groupby and transform

#Se contruye una Series de los datos de AAPL
aapl.head()
ts=aapl.Open

#Se revisa el tipo de datos y la forma
type(ts)
ts.head()

kye=lambda x:x.year
zscore:lambda x:(x-x.mean())/x.std()

#Caso 1

ts_trans=ts.groupby(key).transform(zscore)

TS_grp=ts_trans.groupby(key)
TS_grp.mean()
TS_grp.std()

#Caso 2

grp=ts.groupby(key)
grp.mean()
grp.std()

#Gráfica de los datos de cierre de los dos inficadores

Indicadores=DataFrame({'AAPL_Close':df1.Close_aapl,'IBM_Close':df1.Close_ibm})
Indicadores.plot()
plt.show()

La gráfica que se obtiene es la siguiente:

Pandas_2

Espero que con este sencillo ejemplo se aprecie que las herramientas de Python para hacer data analysis son bastante buenas, rápidas y dentro de Ipython se tienen prácticamente un entorno de desarrollo solo basta agregarle un block de notas para escribir el código. En la referencia [9] se pueden consultar todos los detalles sobre pandas, muchas cosas al respecto no las menciono pero creo que el mejor lugar para aprender es allí.

Data.Frame-Julia

PENDIENTE DE TERMINAR

 

Tutoriales:

1.-https://www.youtube.com/watch?v=8SGif63VW6E

2.-https://www.youtube.com/watch?v=qLrdYhizEMg

3.-https://www.youtube.com/watch?v=MvH1eTdsekA

Referencias:

1.-http://www.jstatsoft.org/v21/i12/paper

2.-http://had.co.nz/thesis/practical-tools-hadley-wickham.pdf

3.-http://cran.r-project.org/web/packages/reshape2/reshape2.pdf

4.-http://www.jstatsoft.org/v40/i01/paper

5.-http://static.googleusercontent.com/media/research.google.com/es//archive/mapreduce-osdi04.pdf

6.http://cran.rstudio.com/web/packages/dplyr/vignettes/introduction.html

7.-http://cran.r-project.org/web/packages/data.table/vignettes/datatable-intro.pdf

8.-http://blog.yhathq.com/posts/fast-summary-statistics-with-data-dot-table.html

9.-http://pandas.pydata.org/pandas-docs/stable/index.html

10.-http://shop.oreilly.com/product/0636920023784.do

11.-http://interactivepython.org/runestone/static/pythonds/index.html

Algo sobre reducción de dimensiones y datos de sistemas no-lineales.

Debido a la facilidad con la que se pueden adquirir datos y con la cantidad que se generan a diario, también empezaron aparecer técnicas para manipularlos , para explorarlos y procesarlos.Más aún hasta se han bautizado nuevos capos como “data sciences” o “big data“.

Sin entrar en detalles técnicos sobre las “tecnologías” que se tienen y pensando más en los aspectos entre teóricos y prácticos de los métodos para hacer análisis, resulta abrumador la cantidad de herramientas con las que se cuenta ahora. Por inclinación e interés personal tenía curiosidad de los aspectos teóricos de “variedades de aprendizaje”  (manifold learning), principalmente lo desarrollado por Mikhaile Belkin y Partha Niyogi, siendo este último uno de los que a mi parecer pusieron muy claro como obtener aspectos sumamente geométricos  y topológicos de las variedades, el trabajo realizado por ellos fue principalmente la técnica “Laplacian EigenMaps” , pero sin duda vale la pena no solo revisar ese trabajo sino todo lo posible de lo hecho por P. Niyogi. En otro momento hablo sobre esos trabajos.

Pero la intensión en este caso no es hablar de la técnica antes mencionada, sino de una “herramienta” hecha por un alumno del que  quizás  es el gurú a seguir hoy día sobre trabajos de procesamientos de datos o técnicas no lineales ,  J. Tenenmbaum del MIT. Su alumno D.N. Reshef  desarrollado algo que hoy podría llamarse , según Terry Speed, “la correlación del siglo 21”. En un artículo de la revista Sciences Terry hace una explicación de la diferencia entre la correlación o índice de correlación de Pearson que de cierto modo captura la relación lineal entre variables aleatorias y la cual fue y ha sido usada de manera general por años. Pero la propuesta de Reshef, es algo llamado  “coeficiente de información máxima“, MIC por sus siglas en ingles. Este coeficiente es un tipo de generalización de la correlación, es principalmente mucho mas sofisticado en cuanto a terminología y su origen vienen del concepto de “información mutua” ( mutual information) concepto desarrollado por C.Shannon en Teoría de la Información, pero más aún este coeficiente resulta natural para hacer comparación o análisis cuando la relación entre las variables es no-lineal.

Por lo cual, vale la pena ponerlo a prueba y ver como va madurando tanto el software para hacer su cálculo como sus aspectos teóricos, que aún son recientes y están en proceso. Siguiendo a K Murphy, no estaría mal hacer un comparativo y ver en que tipo de situación arroja información y como entender o interpretar dicho coeficiente.

Así que aprenderlo y usarlo.