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