Instalación e Interfaz

• Lo primero que debe hacer para empezar a utilizar R es instalarlo en su ordenador. R funciona en casi todas las plataformas disponibles, incluidos los sistemas Windows, Mac OS X y Linux. Links de descarga para R y RStudio.

• RStudio es un entorno de desarrollo integrado (IDE) disponible para R, el cual tiene un buen editor con resaltado de sintaxis, un visor de objetos de R y un gran número de características agradables que están integradas.

Data Frames

• Un data frame es una estructura de datos clave en estadística y en R. La estructura básica de un data frame es que hay una observación por fila y cada columna representa una variable, medida, rasgo o característica de esa observación. R tiene una implementación interna de los data frame que es probablemente la que se utiliza con más frecuencia.

El paquete dplyr

• El paquete dplyr fue desarrollado por Hadley Wickham de RStudio y es una versión optimizada y destilada de su paquete plyr. Una importante contribución del paquete dplyr es que proporciona una “gramática” (en particular, verbos) para la manipulación de datos y para operar con data frames. Con esta gramática, se puede comunicar de forma comprensible lo que se está haciendo a un data frame. Esto es útil porque proporciona una abstracción para la manipulación de datos que antes no existía. Otra contribución útil es que las funciones de dplyr son muy rápidas, ya que muchas operaciones clave están codificadas en C++.

Gramática dplyr

• Algunos de los “verbos” clave proporcionados por el paquete dplyr son

  • select: devuelve un subconjunto de las columnas de un dataframe, utilizando una notación flexible
  • filter: extraer un subconjunto de filas de un dataframe basándose en condiciones lógicas
  • arrange: reordenar las filas de un dataframe
  • rename: renombrar las variables de un dataframe
  • mutate: añadir nuevas variables/columnas o transformar variables existentes
  • summarise: generar estadísticas de resumen de diferentes variables en el dataframe, posiblemente dentro de los estratos
  • %>%: el operador “pipe” se utiliza para conectar varias acciones verbales en una pipeline

Propiedades comunes de la función dplyr

• Todas las funciones que trataremos en este capítulo tendrán algunas características comunes. En particular,

  1. El primer argumento es un dataframe.
  2. Los argumentos siguientes describen qué hacer con el dataframe especificado en el primer argumento, y puede referirse a las columnas del dataframe directamente sin usar el operador $ (sólo use los nombres de las columnas).
  3. El resultado de retorno de una función es un nuevo dataframe
  4. Los dataframe deben estar correctamente formateados y anotados para que todo esto sea útil. En resumen, debe haber una observación por fila y cada columna debe representar un rasgo o característica de esa observación.

Instalación del paquete dplyr

• El paquete dplyr se puede instalar desde CRAN o desde GitHub utilizando el paquete devtools y la función install_github(). El repositorio de GitHub suele contener las últimas actualizaciones del paquete y la versión de desarrollo.

• Para instalarlo desde CRAN, basta con ejecutar

install.packages("dplyr")

• Para instalar desde GitHub puedes ejecutar

install_github("hadley/dplyr")

• Después de instalar el paquete es importante que lo cargues en su sesión de R con la función library()

library(dplyr)

• Si la instalación de dplyr no fue realizada correctamente, un mensaje de instalación requerida aparecerá resaltado en amarillo en la parte superior de RStudio. Haga click en Install y espere un momento mientras se lleva a cabo la instalación. Puede usar también esta estrategia de instalación en lo que sigue

• NOTA: Si alguna vez se encuentra con un problema en el que R confunde la función que quiere llamar, puede especificar el nombre completo de una función utilizando el operador ::. Por ejemplo, la función de filter del paquete dplyr tiene el nombre completo dplyr::filter. Llamar a las funciones con su nombre completo resolverá cualquier confusión sobre la función a la que se quiere llamar.

select()

Para los ejemplos de esta sección utilizaremos un dataset que contiene información sobre la contaminación del aire y la temperatura de la ciudad de Chicago (ver Air pollution dataset). Todos los datos usados en este curso los podrás descargar en: RDataSets

• Después de descomprimir el archivo, puede cargar los datos en R utilizando la función readRDS(). Esta función proporciona los medios para guardar y restaurar un único objeto R en una conexión. Esta difiere de save y load, que guardan y restauran uno o más objetos con nombre en un entorno. Son ampliamente utilizadas por el propio R, por ejemplo para almacenar los metadatos de un paquete y para almacenar las bases de datos de help.search: la extensión de archivo ".rds" es la más utilizada.

chicago <- readRDS("C:/Data/chicago.rds")

• Puede ver algunas características básicas del dataset con las funciones dim() y str().

dim(chicago)

## [1] 6940    8

str(chicago)

## 'data.frame':    6940 obs. of  8 variables:
##  $ city      : chr  "chic" "chic" "chic" "chic" ...
##  $ tmpd      : num  31.5 33 33 29 32 40 34.5 29 26.5 32.5 ...
##  $ dptp      : num  31.5 29.9 27.4 28.6 28.9 ...
##  $ date      : Date, format: "1987-01-01" "1987-01-02" ...
##  $ pm25tmean2: num  NA NA NA NA NA NA NA NA NA NA ...
##  $ pm10tmean2: num  34 NA 34.2 47 NA ...
##  $ o3tmean2  : num  4.25 3.3 3.33 4.38 4.75 ...
##  $ no2tmean2 : num  20 23.2 23.8 30.4 30.3 ...

head(chicago)

##   city tmpd   dptp       date pm25tmean2 pm10tmean2 o3tmean2 no2tmean2
## 1 chic 31.5 31.500 1987-01-01         NA   34.00000 4.250000  19.98810
## 2 chic 33.0 29.875 1987-01-02         NA         NA 3.304348  23.19099
## 3 chic 33.0 27.375 1987-01-03         NA   34.16667 3.333333  23.81548
## 4 chic 29.0 28.625 1987-01-04         NA   47.00000 4.375000  30.43452
## 5 chic 32.0 28.875 1987-01-05         NA         NA 4.750000  30.33333
## 6 chic 40.0 35.125 1987-01-06         NA   48.00000 5.833333  25.77233

• La función select()puede utilizarse para seleccionar las columnas de un dataframe en las que desee centrarse. A menudo se tiene un dataframe grande que contiene “todos” los datos, pero un análisis determinado puede utilizar solo un subconjunto de variables u observaciones. La función select() le permite obtener las pocas columnas que pueda necesitar.

• Supongamos que queremos tomar sólo las 3 primeras columnas. Hay varias formas de hacerlo. Podríamos, por ejemplo, utilizar índices numéricos. Pero también podemos utilizar los nombres directamente.

names(chicago)[1:3]

## [1] "city" "tmpd" "dptp"

subset <- select(chicago, names(chicago)[1:3])
head(subset)

##   city tmpd   dptp
## 1 chic 31.5 31.500
## 2 chic 33.0 29.875
## 3 chic 33.0 27.375
## 4 chic 29.0 28.625
## 5 chic 32.0 28.875
## 6 chic 40.0 35.125

subset <- select(chicago, city:dptp)
head(subset)

##   city tmpd   dptp
## 1 chic 31.5 31.500
## 2 chic 33.0 29.875
## 3 chic 33.0 27.375
## 4 chic 29.0 28.625
## 5 chic 32.0 28.875
## 6 chic 40.0 35.125

• Tenga en cuenta que el signo : normalmente no se puede utilizar con nombres o cadenas, pero dentro de la función select() puede utilizarlo para especificar un rango de nombres de variables. También puede omitir variables con la función select() utilizando el signo negativo

subset <- select(chicago, -(city:dptp))
head(subset)

##         date pm25tmean2 pm10tmean2 o3tmean2 no2tmean2
## 1 1987-01-01         NA   34.00000 4.250000  19.98810
## 2 1987-01-02         NA         NA 3.304348  23.19099
## 3 1987-01-03         NA   34.16667 3.333333  23.81548
## 4 1987-01-04         NA   47.00000 4.375000  30.43452
## 5 1987-01-05         NA         NA 4.750000  30.33333
## 6 1987-01-06         NA   48.00000 5.833333  25.77233

• El código equivalente en R base sería el siguiente, el cual es menos intuitivo

i <- match("city", names(chicago))
j <- match("dptp", names(chicago))
head(chicago[, -(i:j)])

##         date pm25tmean2 pm10tmean2 o3tmean2 no2tmean2
## 1 1987-01-01         NA   34.00000 4.250000  19.98810
## 2 1987-01-02         NA         NA 3.304348  23.19099
## 3 1987-01-03         NA   34.16667 3.333333  23.81548
## 4 1987-01-04         NA   47.00000 4.375000  30.43452
## 5 1987-01-05         NA         NA 4.750000  30.33333
## 6 1987-01-06         NA   48.00000 5.833333  25.77233

• La función select() también permite una sintaxis especial que permite especificar nombres de variables basados en patrones. Así, por ejemplo, si quisiéramos conservar todas las variables que terminan con un “2”, podríamos hacer

subset <- select(chicago, ends_with("2"))
str(subset)

## 'data.frame':    6940 obs. of  4 variables:
##  $ pm25tmean2: num  NA NA NA NA NA NA NA NA NA NA ...
##  $ pm10tmean2: num  34 NA 34.2 47 NA ...
##  $ o3tmean2  : num  4.25 3.3 3.33 4.38 4.75 ...
##  $ no2tmean2 : num  20 23.2 23.8 30.4 30.3 ...

• O si quisiéramos mantener todas las variables que empiezan por “d”, podríamos hacer

subset <- select(chicago, starts_with("d"))
str(subset)

## 'data.frame':    6940 obs. of  2 variables:
##  $ dptp: num  31.5 29.9 27.4 28.6 28.9 ...
##  $ date: Date, format: "1987-01-01" "1987-01-02" ...

• También puede utilizar expresiones regulares más generales si es necesario. Consulte la página de ayuda escribiendo en consola (??select) para más detalles.

filter()

• La función filter() se utiliza para extraer subconjuntos de filas de un dataframe. Esta función es similar a la función subset() existente en R, pero es bastante más rápida. Supongamos que queremos extraer las filas del dataframe PM2.5 (la materia particulada 2.5, incluye sustancias químicas orgánicas, polvo, hollín y metales.) que sean superiores a 30 (que es un nivel razonablemente alto), podríamos hacer

chic.f <- filter(chicago, pm25tmean2 > 30)
str(chic.f)

## 'data.frame':    194 obs. of  8 variables:
##  $ city      : chr  "chic" "chic" "chic" "chic" ...
##  $ tmpd      : num  23 28 55 59 57 57 75 61 73 78 ...
##  $ dptp      : num  21.9 25.8 51.3 53.7 52 56 65.8 59 60.3 67.1 ...
##  $ date      : Date, format: "1998-01-17" "1998-01-23" ...
##  $ pm25tmean2: num  38.1 34 39.4 35.4 33.3 ...
##  $ pm10tmean2: num  32.5 38.7 34 28.5 35 ...
##  $ o3tmean2  : num  3.18 1.75 10.79 14.3 20.66 ...
##  $ no2tmean2 : num  25.3 29.4 25.3 31.4 26.8 ...

• Puede ver que ahora sólo hay 194 filas en el dataframe y la distribución de los valores de pm25tmean2 es

summary(chic.f$pm25tmean2)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   30.05   32.12   35.04   36.63   39.53   61.50

• Podemos colocar una secuencia lógica arbitrariamente compleja dentro de filter(), por lo que podríamos, por ejemplo, extraer las filas en las que PM2.5 es mayor que 30 y la temperatura, tmpd es mayor que 80 grados Fahrenheit

chic.f <- filter(chicago, pm25tmean2 > 30 & tmpd > 80)
select(chic.f, date, tmpd, pm25tmean2)

##          date tmpd pm25tmean2
## 1  1998-08-23   81   39.60000
## 2  1998-09-06   81   31.50000
## 3  2001-07-20   82   32.30000
## 4  2001-08-01   84   43.70000
## 5  2001-08-08   85   38.83750
## 6  2001-08-09   84   38.20000
## 7  2002-06-20   82   33.00000
## 8  2002-06-23   82   42.50000
## 9  2002-07-08   81   33.10000
## 10 2002-07-18   82   38.85000
## 11 2003-06-25   82   33.90000
## 12 2003-07-04   84   32.90000
## 13 2005-06-24   86   31.85714
## 14 2005-06-27   82   51.53750
## 15 2005-06-28   85   31.20000
## 16 2005-07-17   84   32.70000
## 17 2005-08-03   84   37.90000

arrange()

• La función arrange() se utiliza para reordenar las filas de un dataframe según una de las variables/columnas. Reordenar las filas de un dataframe (conservando el orden correspondiente de otras columnas) suele ser una tarea ardua en R. La función arrange() simplifica la tarea. Aquí podemos ordenar las filas del dataframe chicago por fecha, de modo que la primera fila sea la observación más antigua y la última la más reciente

chicago <- arrange(chicago, date)

• Ahora podemos visualizar las primeras y últimas filas de las siguiente manera

head(select(chicago, date, pm25tmean2), 3)

##         date pm25tmean2
## 1 1987-01-01         NA
## 2 1987-01-02         NA
## 3 1987-01-03         NA

tail(select(chicago, date, pm25tmean2), 3)

##            date pm25tmean2
## 6938 2005-12-29    7.45000
## 6939 2005-12-30   15.05714
## 6940 2005-12-31   15.00000

• Las columnas también se pueden ordenar de forma descendente utilizando el operador especial desc()

chicago <- arrange(chicago, desc(date))
head(chicago)

##   city tmpd dptp       date pm25tmean2 pm10tmean2  o3tmean2 no2tmean2
## 1 chic   35 30.1 2005-12-31   15.00000       23.5  2.531250  13.25000
## 2 chic   36 31.0 2005-12-30   15.05714       19.2  3.034420  22.80556
## 3 chic   35 29.4 2005-12-29    7.45000       23.5  6.794837  19.97222
## 4 chic   37 34.5 2005-12-28   17.75000       27.5  3.260417  19.28563
## 5 chic   40 33.6 2005-12-27   23.56000       27.0  4.468750  23.50000
## 6 chic   35 29.6 2005-12-26    8.40000        8.5 14.041667  16.81944

• Si se observan las tres primeras y las tres últimas filas, las fechas aparecen en orden descendente

head(select(chicago, date, pm25tmean2), 3)

##         date pm25tmean2
## 1 2005-12-31   15.00000
## 2 2005-12-30   15.05714
## 3 2005-12-29    7.45000

tail(select(chicago, date, pm25tmean2), 3)

##            date pm25tmean2
## 6938 1987-01-03         NA
## 6939 1987-01-02         NA
## 6940 1987-01-01         NA

rename()

• Renombrar una variable en un dataframe en R es sorprendentemente difícil de hacer. La función rename() está diseñada para facilitar este proceso. Aquí puedes ver los nombres de las cinco primeras variables del dataframe chicago.

head(chicago[, 1:5], 3)

##   city tmpd dptp       date pm25tmean2
## 1 chic   35 30.1 2005-12-31   15.00000
## 2 chic   36 31.0 2005-12-30   15.05714
## 3 chic   35 29.4 2005-12-29    7.45000

• La columna dptp se supone que representa la temperatura del punto de rocío y la columna pm25tmean2 proporciona los datos PM2.5. Sin embargo, estos nombres son incómodos y probablemente deban ser renombrados a algo más sencillos. Dew point: temperatura más alta a la que empieza a condensarse el vapor de agua contenido en el aire, produciendo rocío, neblina, cualquier tipo de nube o, en caso de que la temperatura sea lo suficientemente baja, escarcha.

chicago <- rename(chicago, dewpoint = dptp, pm25 = pm25tmean2)
head(chicago[, 1:5], 3)

##   city tmpd dewpoint       date     pm25
## 1 chic   35     30.1 2005-12-31 15.00000
## 2 chic   36     31.0 2005-12-30 15.05714
## 3 chic   35     29.4 2005-12-29  7.45000

• La sintaxis dentro de la función rename() es tener el nuevo nombre en el lado izquierdo del signo = y el nombre antiguo en el lado derecho.

mutate()

• La función mutate() existe para calcular las transformaciones de las variables de un dataframe. A menudo, se desea crear nuevas variables derivadas de variables existentes y mutate() proporciona una interfaz limpia para hacerlo. Por ejemplo, con los datos de contaminación atmosférica, a menudo queremos eliminar la tendencia de los datos restando la media de los mismos

• De este modo, podemos ver si el nivel de contaminación atmosférica de un día determinado es superior o inferior a la media (en lugar de observar su nivel absoluto). Aquí creamos una variable pm25detrend que resta la media de la variable pm25.

chicago <- mutate(chicago, pm25detrend = pm25 - mean(pm25, na.rm = TRUE))
head(chicago)

##   city tmpd dewpoint       date     pm25 pm10tmean2  o3tmean2 no2tmean2
## 1 chic   35     30.1 2005-12-31 15.00000       23.5  2.531250  13.25000
## 2 chic   36     31.0 2005-12-30 15.05714       19.2  3.034420  22.80556
## 3 chic   35     29.4 2005-12-29  7.45000       23.5  6.794837  19.97222
## 4 chic   37     34.5 2005-12-28 17.75000       27.5  3.260417  19.28563
## 5 chic   40     33.6 2005-12-27 23.56000       27.0  4.468750  23.50000
## 6 chic   35     29.6 2005-12-26  8.40000        8.5 14.041667  16.81944
##   pm25detrend
## 1   -1.230958
## 2   -1.173815
## 3   -8.780958
## 4    1.519042
## 5    7.329042
## 6   -7.830958

• También existe la función relacionada llamada transmute(), que hace lo mismo que mutate() pero luego elimina todas las variables no transformadas. En este caso, se remueve la tendencia de las variables PM10 y ozono (O3)

head(transmute(chicago, 
               pm10detrend = pm10tmean2 - mean(pm10tmean2, na.rm = TRUE),
               o3detrend = o3tmean2 - mean(o3tmean2, na.rm = TRUE)))

##   pm10detrend  o3detrend
## 1  -10.395206 -16.904263
## 2  -14.695206 -16.401093
## 3  -10.395206 -12.640676
## 4   -6.395206 -16.175096
## 5   -6.895206 -14.966763
## 6  -25.395206  -5.393846

• Observe que sólo hay dos columnas en el dataframe transmutado.

group_by()

• La función group_by() se utiliza para generar estadísticas de resumen del dataframe dentro de los estatus definidos por una variable. Por ejemplo, en este dataset de contaminación atmosférica, podría querer saber cuál es el nivel medio anual de PM2,5. Entonces el estatus es el año, y eso es algo que podemos derivar de la variable fecha. Junto con la función group_by() se suele utilizar la función summarize() (o summarise()` para algunas partes del mundo).

• La operación general aquí es una combinación de dividir un dataframe en piezas separadas definidas por una variable o grupo de variables (group_by()), y luego aplicar una función de resumen en esos subconjuntos (summarize()). En primer lugar, podemos crear una variable de año utilizando as.POSIXlt(). Las opciones para días y meses respectivamente, están disponibles usando $mday y mon+1 (ver as.POSIXlt). Los valores anuales se almacenan utilizando un valor de índice base de 1900. Así, 2015 se almacena como 115 ($year = 115).

chicago <- mutate(chicago, year = as.POSIXlt(date)$year + 1900)
head(chicago)

##   city tmpd dewpoint       date     pm25 pm10tmean2  o3tmean2 no2tmean2
## 1 chic   35     30.1 2005-12-31 15.00000       23.5  2.531250  13.25000
## 2 chic   36     31.0 2005-12-30 15.05714       19.2  3.034420  22.80556
## 3 chic   35     29.4 2005-12-29  7.45000       23.5  6.794837  19.97222
## 4 chic   37     34.5 2005-12-28 17.75000       27.5  3.260417  19.28563
## 5 chic   40     33.6 2005-12-27 23.56000       27.0  4.468750  23.50000
## 6 chic   35     29.6 2005-12-26  8.40000        8.5 14.041667  16.81944
##   pm25detrend year
## 1   -1.230958 2005
## 2   -1.173815 2005
## 3   -8.780958 2005
## 4    1.519042 2005
## 5    7.329042 2005
## 6   -7.830958 2005

• Ahora podemos crear un dataframe separado que divida el dataset original por año

years <- group_by(chicago, year)
years

## # A tibble: 6,940 × 10
## # Groups:   year [19]
##    city   tmpd dewpoint date        pm25 pm10tmean2 o3tmean2 no2tmean2
##    <chr> <dbl>    <dbl> <date>     <dbl>      <dbl>    <dbl>     <dbl>
##  1 chic     35     30.1 2005-12-31 15          23.5     2.53      13.2
##  2 chic     36     31   2005-12-30 15.1        19.2     3.03      22.8
##  3 chic     35     29.4 2005-12-29  7.45       23.5     6.79      20.0
##  4 chic     37     34.5 2005-12-28 17.8        27.5     3.26      19.3
##  5 chic     40     33.6 2005-12-27 23.6        27       4.47      23.5
##  6 chic     35     29.6 2005-12-26  8.4         8.5    14.0       16.8
##  7 chic     35     32.1 2005-12-25  6.7         8      14.4       13.8
##  8 chic     37     35.2 2005-12-24 30.8        25.2     1.77      32.0
##  9 chic     41     32.6 2005-12-23 32.9        34.5     6.91      29.1
## 10 chic     22     23.3 2005-12-22 36.6        42.5     5.39      33.7
## # ℹ 6,930 more rows
## # ℹ 2 more variables: pm25detrend <dbl>, year <dbl>

• Por último, calculamos los estadísticos de resumen para cada año del dataframe con la función summarize(). La función summarize() devuelve un dataframe con year como primera columna, y luego los promedios anuales de pm25, o3 y no2.

summarize(years, pm25 = mean(pm25, na.rm = TRUE),
          o3 = max(o3tmean2, na.rm = TRUE),
          no2 = median(no2tmean2, na.rm = TRUE))

## # A tibble: 19 × 4
##     year  pm25    o3   no2
##    <dbl> <dbl> <dbl> <dbl>
##  1  1987 NaN    63.0  23.5
##  2  1988 NaN    61.7  24.5
##  3  1989 NaN    59.7  26.1
##  4  1990 NaN    52.2  22.6
##  5  1991 NaN    63.1  21.4
##  6  1992 NaN    50.8  24.8
##  7  1993 NaN    44.3  25.8
##  8  1994 NaN    52.2  28.5
##  9  1995 NaN    66.6  27.3
## 10  1996 NaN    58.4  26.4
## 11  1997 NaN    56.5  25.5
## 12  1998  18.3  50.7  24.6
## 13  1999  18.5  57.5  24.7
## 14  2000  16.9  55.8  23.5
## 15  2001  16.9  51.8  25.1
## 16  2002  15.3  54.9  22.7
## 17  2003  15.2  56.2  24.6
## 18  2004  14.6  44.5  23.4
## 19  2005  16.2  58.8  22.6

• En un ejemplo algo más trabajado, podríamos querer saber cuáles son los niveles medios de de ozono (o3) y dióxido de nitrógeno (no2) dentro de los quintiles de pm25. Una forma más sencilla de de hacerlo sería a través de un modelo de regresión, pero en realidad podemos hacerlo rápidamente con group_by() y summarize().

• Primero, podemos crear una variable categórica de pm25 dividida en quintiles La función cut divide el rango en intervalos y codifica los valores según el intervalo al que pertenecen.

qq <- quantile(chicago$pm25, seq(0, 1, 0.20), na.rm = TRUE)
qq

##     0%    20%    40%    60%    80%   100% 
##  1.700  8.700 12.375 16.700 22.610 61.500

chicago <- mutate(chicago, pm25.quint = cut(pm25, qq))
head(chicago)

##   city tmpd dewpoint       date     pm25 pm10tmean2  o3tmean2 no2tmean2
## 1 chic   35     30.1 2005-12-31 15.00000       23.5  2.531250  13.25000
## 2 chic   36     31.0 2005-12-30 15.05714       19.2  3.034420  22.80556
## 3 chic   35     29.4 2005-12-29  7.45000       23.5  6.794837  19.97222
## 4 chic   37     34.5 2005-12-28 17.75000       27.5  3.260417  19.28563
## 5 chic   40     33.6 2005-12-27 23.56000       27.0  4.468750  23.50000
## 6 chic   35     29.6 2005-12-26  8.40000        8.5 14.041667  16.81944
##   pm25detrend year  pm25.quint
## 1   -1.230958 2005 (12.4,16.7]
## 2   -1.173815 2005 (12.4,16.7]
## 3   -8.780958 2005   (1.7,8.7]
## 4    1.519042 2005 (16.7,22.6]
## 5    7.329042 2005 (22.6,61.5]
## 6   -7.830958 2005   (1.7,8.7]

• Ahora podemos agrupar el dataframe por la variable pm25.quint.

quint <- group_by(chicago, pm25.quint)
quint

## # A tibble: 6,940 × 11
## # Groups:   pm25.quint [6]
##    city   tmpd dewpoint date        pm25 pm10tmean2 o3tmean2 no2tmean2
##    <chr> <dbl>    <dbl> <date>     <dbl>      <dbl>    <dbl>     <dbl>
##  1 chic     35     30.1 2005-12-31 15          23.5     2.53      13.2
##  2 chic     36     31   2005-12-30 15.1        19.2     3.03      22.8
##  3 chic     35     29.4 2005-12-29  7.45       23.5     6.79      20.0
##  4 chic     37     34.5 2005-12-28 17.8        27.5     3.26      19.3
##  5 chic     40     33.6 2005-12-27 23.6        27       4.47      23.5
##  6 chic     35     29.6 2005-12-26  8.4         8.5    14.0       16.8
##  7 chic     35     32.1 2005-12-25  6.7         8      14.4       13.8
##  8 chic     37     35.2 2005-12-24 30.8        25.2     1.77      32.0
##  9 chic     41     32.6 2005-12-23 32.9        34.5     6.91      29.1
## 10 chic     22     23.3 2005-12-22 36.6        42.5     5.39      33.7
## # ℹ 6,930 more rows
## # ℹ 3 more variables: pm25detrend <dbl>, year <dbl>, pm25.quint <fct>

• Por último, podemos calcular la media de o3 y no2 dentro de los quintiles de pm25

na.omit(summarize(quint, o3 = mean(o3tmean2, na.rm = TRUE), no2 = mean(no2tmean2, na.rm = TRUE)))

## # A tibble: 5 × 3
##   pm25.quint     o3   no2
##   <fct>       <dbl> <dbl>
## 1 (1.7,8.7]    21.7  18.0
## 2 (8.7,12.4]   20.4  22.1
## 3 (12.4,16.7]  20.7  24.4
## 4 (16.7,22.6]  19.9  27.3
## 5 (22.6,61.5]  20.3  29.6

• Con base en la tabla, no hay una relación fuerte entre pm25 y o3, pero parece haber una correlación positiva entre pm25 y no2. Un modelo estadístico más sofisticado puede ayudar a dar respuestas precisas a estas preguntas, pero una simple aplicación de las funciones dplyr puede llevarle a la mayoría de los caminos.

• Nótese además que si omitimos el uso de la función na.omit, podremos observar una fila cuya clase corresponde a <NA>. Esto está asociado a observaciones para las cuales existen valores de o3 y no2, pero, pm25 es faltante. Estudiaremos mas adelante como realizar imputación de este tipo de datos.

summarize(quint, o3 = mean(o3tmean2, na.rm = TRUE), no2 = mean(no2tmean2, na.rm = TRUE))

## # A tibble: 6 × 3
##   pm25.quint     o3   no2
##   <fct>       <dbl> <dbl>
## 1 (1.7,8.7]    21.7  18.0
## 2 (8.7,12.4]   20.4  22.1
## 3 (12.4,16.7]  20.7  24.4
## 4 (16.7,22.6]  19.9  27.3
## 5 (22.6,61.5]  20.3  29.6
## 6 <NA>         18.8  25.8

%>%

• El operador pipeline %>% es muy útil para encadenar múltiples funciones dplyr en una secuencia de operaciones. Nótese que hasta ahora, cada vez que se desea aplicar más de una función, la secuencia quedaba dentro de una secuencia de llamadas a funciones anidadas que es difícil de leer, por ejemplo

third(second(first(x)))

• Este anidamiento no es una forma natural de pensar en una secuencia de operaciones. El operador %>% permite encadenar las operaciones de izquierda a derecha, es decir

first(x) %>% 
  second %>% 
  third

• Tomemos el ejemplo que acabamos de hacer en la última sección, en el que hemos calculado la media de o3 y no2 dentro de los quintiles de pm25. Allí tuvimos que

  1. Crear una nueva variable pm25.quint
  2. Agrupar el dataframe por esa nueva variable
  3. Calcular la media de o3 y no2 en los subgrupos definidos por pm25.quint

• Eso se puede hacer con la siguiente secuencia en una sola expresión de R

library(tidyr)

mutate(chicago, pm25.quint = cut(pm25, qq)) %>%
  group_by(pm25.quint) %>%
  summarize(o3 = mean(o3tmean2, na.rm = TRUE), no2 = mean(no2tmean2, na.rm = TRUE)) %>%
  drop_na()

## # A tibble: 5 × 3
##   pm25.quint     o3   no2
##   <fct>       <dbl> <dbl>
## 1 (1.7,8.7]    21.7  18.0
## 2 (8.7,12.4]   20.4  22.1
## 3 (12.4,16.7]  20.7  24.4
## 4 (16.7,22.6]  19.9  27.3
## 5 (22.6,61.5]  20.3  29.6

• De esta manera no tenemos que crear un conjunto de variables temporales en el camino o crear una secuencia masiva anidada de llamadas a funciones. Observe en el código anterior que el dataframe chicago pasa a la primera llamada a mutate(), pero después no tiene que pasar el primer argumento a group_by() o a summarize(). Una vez que se recorre el pipeline con %>%, el primer argumento se toma como la salida del elemento anterior en la pipeline.

• Otro ejemplo podría ser el cálculo del nivel medio de contaminantes pm25, o3, no2 por mes. Esto podría ser útil para ver si hay alguna tendencia estacional en los datos.

mutate(chicago, month = as.POSIXlt(date)$mon + 1) %>%
  group_by(month) %>%
  summarize(pm25 = mean(pm25, na.rm = TRUE),
            o3 = max(o3tmean2, na.rm = TRUE),
            no2 = median(no2tmean2, na.rm = TRUE))

## # A tibble: 12 × 4
##    month  pm25    o3   no2
##    <dbl> <dbl> <dbl> <dbl>
##  1     1  17.8  28.2  25.4
##  2     2  20.4  37.4  26.8
##  3     3  17.4  39.0  26.8
##  4     4  13.9  47.9  25.0
##  5     5  14.1  52.8  24.2
##  6     6  15.9  66.6  25.0
##  7     7  16.6  59.5  22.4
##  8     8  16.9  54.0  23.0
##  9     9  15.9  57.5  24.5
## 10    10  14.2  47.1  24.2
## 11    11  15.2  29.5  23.6
## 12    12  17.5  27.7  24.5

• Aquí podemos ver que el o3 tiende a ser bajo en los meses de invierno (11, 12, 1, 2, 3) y alto en el verano (4,5,…,10) mientras que no2 es más alto en invierno y más bajo en verano.

Graficar un dataset

• El paquete ggplot2 puede utilizarse para crear varios tipos de gráficos y diagramas a partir de un dataframe. La función ggplot() se utiliza para definir gráficos, y se le pueden pasar una serie de parámetros y unirse a otras funciones para producir finalmente un gráfico de salida. El primer parámetro que se pasa a ggplot() es el dataframe que se quiere representar. Normalmente será un objeto data.frame, pero también puede ser un subconjunto de un dataframe definido con la función subset().

library(readr)
dfCrime = read_csv("C:/Data/Crime_Data.csv", col_names = TRUE)

• Puede obtener un recuento del número de registros con la función nrow(). La función View() se puede utilizar para ver los datos en un formato tabular.

nrow(dfCrime)

## [1] 481376

head(dfCrime)

## # A tibble: 6 × 11
##   `Report Number` `Occurred Date` `Occurred Time` `Reported Date`
##             <dbl> <chr>                     <dbl> <chr>          
## 1         2.01e13 12/13/1908                 2114 12/13/2008     
## 2         2.01e13 06/15/1964                    0 06/15/2010     
## 3         2.01e12 01/01/1973                    0 01/25/2012     
## 4         2.01e13 06/01/1974                    0 09/09/2013     
## 5         2.02e13 01/01/1975                    0 08/11/2016     
## 6         1.98e12 12/16/1975                  900 12/16/1975     
## # ℹ 7 more variables: `Reported Time` <dbl>, `Crime Subcategory` <chr>,
## #   `Primary Offense Description` <chr>, Precinct <chr>, Sector <chr>,
## #   Beat <chr>, Neighborhood <chr>

• Un tibble, o tbl_df, es una reimaginación moderna del data.frame, manteniendo lo que el tiempo ha demostrado que es eficaz, y desechando lo que no lo es. En este caso limitaremos las columnas a lo siguiente: Reported Date, Crime Subcategory, Primary Offense Description, Precinct, Sector, Beat y Neighborhood

• tibble tiene un método de imprimir refinado que muestra sólo las 10 primeras filas y todas las columnas que caben en la pantalla. Esto hace que sea mucho más fácil trabajar con datos de gran tamaño. A diferencia de los data.frame, los tibble no admiten operaciones aritméticas en todas las columnas, mientras que con un data.frame.

dfCrime = select(dfCrime,
                 "Reported Date",
                 "Crime Subcategory",
                 "Primary Offense Description",
                 "Precinct",
                 "Sector",
                 "Beat",
                 "Neighborhood")

• También es posible que quiera cambiar el nombre de las columnas para que sean más fáciles de leer o simplificar los nombres. La función select() se puede utilizar para hacer esto también.

dfCrime = select(dfCrime, 
                 "CrimeDate" = "Reported Date", 
                 "Category" = "Crime Subcategory", 
                 "Description" = "Primary Offense Description", 
                 "Precinct",
                 "Sector", 
                 "Beat", 
                 "Neighborhood")
head(dfCrime)

## # A tibble: 6 × 7
##   CrimeDate  Category             Description Precinct Sector Beat  Neighborhood
##   <chr>      <chr>                <chr>       <chr>    <chr>  <chr> <chr>       
## 1 12/13/2008 DUI                  DUI-LIQUOR  EAST     G      G2    CENTRAL ARE…
## 2 06/15/2010 FAMILY OFFENSE-NONV… CHILD-OTHER WEST     Q      Q2    QUEEN ANNE  
## 3 01/25/2012 SEX OFFENSE-OTHER    SEXOFF-OTH… NORTH    N      N2    NORTHGATE   
## 4 09/09/2013 SEX OFFENSE-OTHER    SEXOFF-OTH… UNKNOWN  <NA>   <NA>  UNKNOWN     
## 5 08/11/2016 SEX OFFENSE-OTHER    SEXOFF-OTH… UNKNOWN  <NA>   <NA>  UNKNOWN     
## 6 12/16/1975 BURGLARY-RESIDENTIAL BURGLARY-F… SOUTH    R      R3    LAKEWOOD/SE…

• El filtrado del dataset le permite centrarse en un subconjunto de filas en lugar de en todo el dataset. El paquete dplyr incluye una función filter() que soporta esta capacidad

dfCrime2 = filter(dfCrime, Neighborhood == "QUEEN ANNE")
nrow(dfCrime2)

## [1] 25172

• También puede incluir varias expresiones en una función filter(), como sigue, vecindario QUEEN ANNE y categoría BURGLARY-RESIDENTIAL. Examinaremos expresiones de filtro más complejas en una sección posterior

dfCrime3 = filter(dfCrime, 
                  Neighborhood == "QUEEN ANNE", 
                  Category == "BURGLARY-RESIDENTIAL")
head(dfCrime3)

## # A tibble: 6 × 7
##   CrimeDate  Category             Description Precinct Sector Beat  Neighborhood
##   <chr>      <chr>                <chr>       <chr>    <chr>  <chr> <chr>       
## 1 04/22/2008 BURGLARY-RESIDENTIAL BURGLARY-N… WEST     Q      Q3    QUEEN ANNE  
## 2 01/02/2008 BURGLARY-RESIDENTIAL BURGLARY-F… WEST     Q      Q3    QUEEN ANNE  
## 3 01/02/2008 BURGLARY-RESIDENTIAL BURGLARY-F… WEST     D      D2    QUEEN ANNE  
## 4 01/07/2008 BURGLARY-RESIDENTIAL BURGLARY-F… WEST     Q      Q3    QUEEN ANNE  
## 5 01/09/2008 BURGLARY-RESIDENTIAL BURGLARY-N… WEST     Q      Q2    QUEEN ANNE  
## 6 01/10/2008 BURGLARY-RESIDENTIAL BURGLARY-F… WEST     Q      Q3    QUEEN ANNE

• En la ventana de la consola de RStudio añada el código que ve a continuación para agrupar los delitos por ronda policial

dfCrime2 = group_by(dfCrime2, Beat)

• La función n() se utiliza para obtener un recuento del número de registros de cada grupo. Añada a el código que se ve a continuación. La columna Beat corresponde a el límite del sector policial designado donde se produjeron las infracciones (ver Seattle Police Department)

dfCrime2 = summarise(dfCrime2, n = n())
head(dfCrime2)

## # A tibble: 4 × 2
##   Beat      n
##   <chr> <int>
## 1 D2     4373
## 2 Q1       88
## 3 Q2    10851
## 4 Q3     9860

Gráficos con ggplot()

• En la consola de RStudio añada el código que ve a continuación. La función ggplot() en este caso recibe el dataframe dfCrime. La función geom_col() se utiliza para definir la geometría del gráfico (diagrama de barras) y se le pasa un parámetro de asignación que se define llamando a la función aes() y pasando las columnas para el eje x (Beat), y el eje y (n = count)

library("ggplot2")

ggplot(data=dfCrime2) + 
  geom_col(mapping = aes(x = Beat, y = n), fill = "red")

Gráficos de robos por mes y año

• En este ejercicio crearemos un par de gráficos de barras que muestren el número de robos por año y por mes para el barrio de Queen Anne. Además de los paquetes dplyr y ggplot2 que que usamos anteriormente, también usaremos el paquete lubridate para manipular la información de las fechas.

• En el panel de paquetes de RStudio, cargue el paquete lubridate. El paquete lubridate forma parte de tidyverse y se utiliza para trabajar con fechas y horas. Además, asegúrese de que los paquetes readr, dplyr y ggplot2 están cargados

library(readr)
dfCrime = read_csv("C:/Data/Crime_Data.csv", col_names = TRUE)

• Especifique las columnas y los nombres de las columnas

dfCrime = select(dfCrime, 
                 "CrimeDate" = "Reported Date",
                 "Category" = "Crime Subcategory",
                 "Description" = "Primary Offense Description", 
                 "Precinct",
                 "Sector", 
                 "Beat", 
                 "Neighborhood")

• Filtrar los registros para que sólo se conserven los robos residenciales en el barrio de Queen Anne

dfCrime2 = filter(dfCrime, 
                  Neighborhood == "QUEEN ANNE", 
                  Category == "BURGLARY-RESIDENTIAL")

• El paquete dplyr incluye la capacidad de crear dinámicamente nuevas columnas en un dataframe mediante la manipulación de los datos de las columnas existentes en el dataframe. La función mutate() se utiliza para crear las nuevas columnas. Aquí la función mutate() se utilizará para extraer el año de la columna CrimeDate.

• Agregue el siguiente código para ver esto en acción. El segundo parámetro crea una nueva columna llamada YEAR y la rellena utilizando la función year() del paquete lubridate. Dentro de la función year() la columna CrimeDate, que es una columna de caracteres, se convierte en una fecha con base en el formato format="%m/%d/%Y".

library(lubridate)
dfCrime3 = mutate(dfCrime2, YEAR = year(as.Date(dfCrime2$CrimeDate,format="%m/%d/%Y")))

• Observe la columna YEAR al final del dataframe. La función mutate() siempre añade nuevas columnas al final del dataframe

head(dfCrime3)

## # A tibble: 6 × 8
##   CrimeDate  Category       Description Precinct Sector Beat  Neighborhood  YEAR
##   <chr>      <chr>          <chr>       <chr>    <chr>  <chr> <chr>        <dbl>
## 1 04/22/2008 BURGLARY-RESI… BURGLARY-N… WEST     Q      Q3    QUEEN ANNE    2008
## 2 01/02/2008 BURGLARY-RESI… BURGLARY-F… WEST     Q      Q3    QUEEN ANNE    2008
## 3 01/02/2008 BURGLARY-RESI… BURGLARY-F… WEST     D      D2    QUEEN ANNE    2008
## 4 01/07/2008 BURGLARY-RESI… BURGLARY-F… WEST     Q      Q3    QUEEN ANNE    2008
## 5 01/09/2008 BURGLARY-RESI… BURGLARY-N… WEST     Q      Q2    QUEEN ANNE    2008
## 6 01/10/2008 BURGLARY-RESI… BURGLARY-F… WEST     Q      Q3    QUEEN ANNE    2008

• Ahora agruparemos los datos por año y resumiremos obteniendo un recuento del número de delitos por año

dfCrime4 = group_by(dfCrime3, YEAR)
dfCrime4 = summarise(dfCrime4, n = n())
head(dfCrime4)

## # A tibble: 6 × 2
##    YEAR     n
##   <dbl> <int>
## 1  2008   205
## 2  2009   194
## 3  2010   181
## 4  2011   203
## 5  2012   170
## 6  2013   209

• Cree un gráfico de barras llamando a las funciones ggplot() y geom_col() como se ve abajo. Defina YEAR como columna para el eje x y el número de delitos para el eje el eje y. Esto debería producir el gráfico que se ve a continuación

ggplot(data=dfCrime4) + geom_col(mapping = aes(x=YEAR, y=n), fill="red")

• Ahora crearemos otro gráfico de barras que muestre el número de delitos por por mes en lugar de por año. En primer lugar, cree una columna MONTH utilizando la función mutate()

dfCrime3 = mutate(dfCrime2, MONTH = month(as.Date(dfCrime2$CrimeDate, format="%m/%d/%Y")))

• Agrupamos y resumimos los datos por mes

dfCrime4 = group_by(dfCrime3, MONTH) 
dfCrime4 = summarise(dfCrime4, n = n())
dfCrime4

## # A tibble: 12 × 2
##    MONTH     n
##    <dbl> <int>
##  1     1   211
##  2     2   162
##  3     3   173
##  4     4   188
##  5     5   206
##  6     6   203
##  7     7   196
##  8     8   144
##  9     9   187
## 10    10   154
## 11    11   168
## 12    12   201

• Para obtener una mejor visualización, podemos por ejemplo utilizar la función month.abb[] para cambiar los id de los meses a las iniciales correspondientes

dfCrime4 = mutate(dfCrime4, MONTH_S = month.abb[MONTH])
dfCrime4

## # A tibble: 12 × 3
##    MONTH     n MONTH_S
##    <dbl> <int> <chr>  
##  1     1   211 Jan    
##  2     2   162 Feb    
##  3     3   173 Mar    
##  4     4   188 Apr    
##  5     5   206 May    
##  6     6   203 Jun    
##  7     7   196 Jul    
##  8     8   144 Aug    
##  9     9   187 Sep    
## 10    10   154 Oct    
## 11    11   168 Nov    
## 12    12   201 Dec

• A continuación, se crea el gráfico de barras

library(scales)

dfCrime4$MONTH_S <- factor(dfCrime4$MONTH_S, levels=unique(dfCrime4$MONTH_S))

ggplot(data=dfCrime4) + geom_col(mapping = aes(x=MONTH_S, y=n), fill="blue") +
  xlab("Month") + ylab("Frequency") +
  ggtitle("Number of crimes per month")

Conclusiones

• En este capítulo se han abordado algunas técnicas básicas para la exploración y visualización utilizando el paquete tidyverse y su ecosistema de subpaquetes. Después de instalar y cargar el paquete mediante RStudio, se realizaron una serie de tareas usando el lenguaje de programación R con una serie de subpaquetes tidyverse. Se cargó un dataset desde un archivo CSV utilizando readr.
• Después, se manipuló los datos de varias maneras utilizando el paquete dplyr. Se utilizó la función select() para incluir y renombrar columnas. El contenido del dataframe se filtró con la función filter(). A continuación, los datos se agruparon y finalmente se produjeron varios gráficos utilizando ggplot2. En la siguiente sección se estudiará más sobre cómo usar el paquete readr para cargar datos de fuentes de datos externas.

Lectura de archivos csv con read.table()

• La primera función que examinaremos es read.table(). La función read.table() es una función función integrada en R que puede utilizarse para leer varios formatos de archivo en un dataframe, esta es probablemente la función interna más común utilizada para leer archivos simples en R. Sin embargo, como veremos más adelante en esta sección, tidyverse incluye funciones similares que son más eficientes en la lectura de datos externos en R. La sintaxis de read.table() es aceptar un nombre de archivo, que será la ruta junto con un indicador TRUE|FALSE para la cabecera. Si se establece como TRUE, se asume que los nombres de las columnas están en la línea de cabecera del archivo. La ruta no es no es necesaria si ya se ha establecido el directorio de trabajo. La salida de la función read.table() es un objeto dataframe. La función file.choose() es una función práctica que puede utilizar para seleccionar interactivamente el archivo que desea importar en lugar de tener que codificar la ruta del dataset. En este ejercicio aprenderá a utilizar la función read.table() para cargar un archivo con formato csv

• Abra RStudio y busque el panel de la consola. Si es necesario, establezca el directorio de trabajo escribiendo el código que ve a continuación en el panel de la consola o desde las opcciones Sesión | Establecer directorio de trabajo | Elegir Directorio en el menú de RStudio.

setwd(<installation directory for exercise data>)

• La carpeta RDataSets contiene un archivo llamado StudyArea.csv, que es un archivo separado por comas que contiene datos de incendios forestales de los años 1980-2016 para los estados de California, Oregón, Washington, Idaho, Montana, Wyoming, Colorado, Utah, Nevada, Arizona y Nuevo México. Hay un poco más de 439.000 registros en este archivo y hay 37 columnas de información que describen cada incendio durante este periodo.

• Utilice la función read.table() para cargar estos datos en un nuevo objeto dataframe

df = read.table("C:/Data/StudyArea.csv", header = TRUE)

• Recibirá un mensaje de error cuando intente ejecutar esta línea de código. El mensaje de error mensaje de error debería aparecer como se ve a continuación

• La razón por la que se generó un mensaje de error en este caso es que la función read.table() utiliza espacios como delimitador entre registros y nuestro archivo utiliza comas como delimitador

• Actualice su llamada a read.table() como se ve a continuación para incluir el argumento sep, que debería ser una coma

df = read.table("C:/Data/StudyArea.csv", sep=",", header = TRUE)

• Cuando ejecute esta línea de código verá un nuevo error

• La función read.table() NO rellenará automáticamente los valores que falten con un valor por defecto como NA, así que como algunas de las columnas están vacías en nuestras filas obtenemos un mensaje de error que indica que una línea en particular no tiene las 14 columnas de información. Podemos arreglar esto añadiendo el parámetro fill y poniéndolo igual a TRUE. Actualice su código como se ve a continuación para añadir el parámetro de fill

df = read.table("C:/Data/StudyArea.csv", header=TRUE, fill=TRUE, sep=",")
nrow(df)

## [1] 153095

• Cuando ejecute esta línea de código, importará el contenido del archivo a un objeto objeto dataframe. Sin embargo, si mira la pestaña Global Enviroment en RStudio verá que sólo se han cargado 153.095 registros y sin embargo sabemos que hay más de 400.000.

• Las comillas (simples o dobles) en un archivo csv pueden hacer que no se carguen los registros. Agreguemos un parámetro más para manejar los registros que fueron descartados debido a comillas.

df = read.table("C:/Data/StudyArea.csv", header=TRUE, fill=TRUE, quote="\"", sep=",")
nrow(df)

## [1] 439362

• Al ejecutar esta línea de código, se deberían importar 440.476 registros. Los datos se cargan en un objeto R data.frame que es una estructura que se asemeja a una tabla. Por ahora, puede pensar en ellos como tablas que contienen columnas y filas.

• La función read.table() se utiliza normalmente para cargar archivos de texto delimitados por tabulaciones, pero muchos intentarán utilizar la función read.table() con archivos de formato csv sin entender todos los parámetros que pueden necesitar ser incluidos. En su lugar, también puede utilizar read.csv() como haremos en el siguiente paso.

Lectura de archivos csv con read.csv()

• La función read.csv() también es una función incorporada en R que es casi idéntica a read.table(), con la excepción de que los argumentos de cabecera y relleno se establecen en TRUE por defecto. En este paso se observa que es mucho más fácil cargar un archivo csv usando read.csv().

• La función read.csv() maneja automáticamente la mayoría de las situaciones que deben identificarse cuando se utiliza read.table() para cargar un archivo csv. Introduzca y ejecute el código que se ve a continuación para ver cuánto más fácil es esto con read.csv(). Esto cargará correctamente los más de 400.000 registros del archivo csv. En general esta función es mucho más fácil de usar que read.table().

df = read.csv("C:/Data/StudyArea.csv")
nrow(df)

## [1] 439362

Lectura de archivo delimitado por tabulaciones con read.table()

• La función read.table() se utiliza más a menudo para leer el contenido de un archivo delimitado por tabulaciones. En este paso aprenderás a hacerlo

• La carpeta RDataSets incluye un archivo llamado all_genes_pombase.txt, que está delimitado por texto. Abra este archivo con Excel o alguna otra aplicación para ver la estructura de los campos y los delimitadores.

• En la ventana de la consola de R introduzca y ejecute el código que ve a continuación para importar el archivo.

df2 = read.table("C:/Data/all_genes_pombase.txt", header=TRUE, sep="\t", quote="\"")
nrow(df2)

## [1] 7019

head(df2)

##     ensembl_id        name chromosome
## 1  SPAC1002.01 SPAC1002.01          I
## 2  SPAC1002.02       pom34          I
## 3 SPAC1002.03c        gls2          I
## 4 SPAC1002.04c       taf11          I
## 5 SPAC1002.05c        jmj2          I
## 6 SPAC1002.06c        bqt2          I
##                                                     description   feature_type
## 1                                     conserved fungal protein  protein_coding
## 2                                            nucleoporin Pom34  protein_coding
## 3                            glucosidase II alpha subunit Gls2  protein_coding
## 4 transcription factor TFIID complex subunit Taf11 (predicted)  protein_coding
## 5                                     histone demethylase Jmj2  protein_coding
## 6                               bouquet formation protein Bqt2  protein_coding
##   strand   start     end
## 1      1 1798347 1799015
## 2      1 1799061 1800053
## 3     -1 1799915 1803141
## 4     -1 1803624 1804491
## 5     -1 1804548 1806797
## 6     -1 1807270 1807781

• Esto debería cargar 7019 registros en el data.frame. Notará que muchos de los parámetros todavía tienen que ser utilizados cuando se carga el dataset, por lo que no es tan fácil de como se podría esperar incluso en este caso.

Uso de readr para cargar datos

• Hasta ahora en este capítulo hemos visto varias funciones incorporadas en R para leer archivos externos como dataframe. El paquete tidyverse incluye un subpaquete llamado readr que también puede utilizarse para cargar datos externos. El paquete readr incluye una función read_csv() que carga los datos mucho más rápido que la función interna read.csv().

• En general, la función read_csv() del paquete readr es preferible a las funciones que se encuentran en la instalación básica de R. El paquete readr también incluye algunas otras funciones para cargar varios formatos de archivo, incluyendo read_delim(), read_csv2(), y read_tsv(). Cada una estas funciones aceptan los mismos parámetros, por lo que una vez que haya aprendido a utilizar cualquiera de las funciones de R para cargar datos, puede fácilmente utilizar cualquiera de las otras.

• En este paso vas a utilizar la función read_csv() que se encuentra en el paquete readr para cargar los datos en un dataframe.

1. Cargar la librería readr

library(readr)

2. La función read_csv() del paquete readr puede utilizarse para cargar archivos csv. En comparación con las funciones de carga base que vimos anteriormente, las funciones readr son significativamente más rápidas (10 veces), incluyen una útil barra de progreso para proporcionar información sobre el progreso de la carga de archivos grandes, y todas las funciones trabajan exactamente de la misma manera.

• Añada y ejecute el código a continuación. Observe que los datos se cargan mucho más rápido en el objeto dataframe. El argumento col_types fue utilizado en este caso para cargar todas las columnas como un tipo de datos de caracteres para simplificar. De lo contrario, tendríamos que hacer un preprocesamiento adicional de los datos para tener en cuenta los distintos tipos de datos de las columnas.

dfReadr = read_csv("C:/Data/StudyArea.csv", col_types = cols(.default = "c"), col_names = TRUE)

3. Ahora vamos a ejecutar esta función de nuevo, pero esta vez quitando el argumento col_types para poder ver un ejemplo de algunos de los posibles errores de carga que pueden ocurrir. Lo primero que verá es una lista de las columnas que se importarán junto con el tipo de datos de la columna. A continuación aparecerá un mensaje de advertencia indicando que hubo errores de análisis en la carga.

dfReadr = read_csv("C:/Data/StudyArea.csv", col_names = TRUE)

4. Puede utilizar la función problems() para obtener una lista de los errores de análisis. Añada y ejecute el código que ves a continuación.

problems(dfReadr)

## # A tibble: 0 × 5
## # ℹ 5 variables: row <int>, col <int>, expected <chr>, actual <chr>, file <chr>

5. Por el aspecto de los mensajes de error parece que hay un problema con la columna UNIT. Si vuelve a mirar la lista de columnas y tipos de datos, observará que la columna UNIT fue creada como un tipo de datos entero. Sin embargo, si abre el archivo StudyArea.csv en Excel o en otra aplicación, verá rápidamente verá que no todos los valores son numéricos, algunos incluyen letras. Esto explica los errores de análisis en el dataset.

• Actualice su código como se ve a continuación y ejecútelo de nuevo. Esto establece la columna UNIT a un tipo de datos de carácter (texto).

dfReadr = read_csv("C:/Data/StudyArea.csv", col_types = list(UNIT = col_character()),
col_names = TRUE)

• Esta vez debería obtener una carga limpia del dataset. Eso no significa que los datos no necesiten alguna preparación y limpieza adicional. Por ejemplo, hay algunos campos de fecha, como STARTDATED, que se cargaron como caracteres pero es mejor que se carguen como campos de fecha. Podemos dejar este trabajo de preparación adicional para un ejercicio posterior.

6. Puede examinar las primeras líneas del dataframe introduciendo la función head() como se ve a continuación.

head(dfReadr)

## # A tibble: 6 × 14
##     FID ORGANIZATI UNIT  SUBUNIT SUBUNIT2        FIRENAME CAUSE YEAR_ STARTDATED
##   <dbl> <chr>      <chr> <chr>   <chr>           <chr>    <chr> <dbl> <chr>     
## 1     0 FWS        81682 USCADBR San Diego Bay … PUMP HO… Human  2001 1/1/01 0:…
## 2     1 FWS        81682 USCADBR San Diego Bay … I5       Human  2002 5/3/02 0:…
## 3     2 FWS        81682 USCADBR San Diego Bay … SOUTHBAY Human  2002 6/1/02 0:…
## 4     3 FWS        81682 USCADBR San Diego Bay … MARINA   Human  2001 7/12/01 0…
## 5     4 FWS        81682 USCADBR San Diego Bay … HILL     Human  1994 9/13/94 0…
## 6     5 FWS        81682 USCADBR San Diego Bay … IRRIGAT… Human  1994 4/22/94 0…
## # ℹ 5 more variables: CONTRDATED <chr>, OUTDATED <chr>, STATE <chr>,
## #   STATE_FIPS <dbl>, TOTALACRES <dbl>

• En esta sección se estudiaron varias funciones para cargar un archivo de datos externo incluyendo las funciones incorporadas en R read.table() y read.csv(). Aunque estas funciones funciones pueden hacer el trabajo, la función read_csv() que se encuentra en el paquete readr es una función mucho más eficiente para cargar datos externos. En la siguiente sección se aprenderá a transformar datasets utilizando el paquete dplyr. Se estudiarán técnicas para filtrar el contenido de un dataframe, seleccionar las columnas específicas que se utilizarán, organizar las filas en orden ascendente o descendente y resumir y agrupar un dataset.

Filtrar registros para crear un subconjunto

• La primera función de dplyr que examinaremos es filter(). La función filter() se utiliza para crear un subconjunto de registros basado en algún valor. Por ejemplo, puede querer crear un dataframe de incendios forestales que contenga incidentes que hayan quemado más de 25.000 acres. Siempre que tenga un dataframe existente que incluya una columna que mida el número de acres quemadas, puede crear este subconjunto utilizando la función filter().

• Como será el caso de todas las funciones de dplyr que examinemos, el primer argumento que se pasa a la función filter() es un objeto de dataframe. Cada parámetro adicional pasado a la función es una expresión condicional utilizada para filtrar el dataframe.

• Por ejemplo, observe la siguiente línea de código. Esta declaración llama a la función filter() para crear una nueva variable llamada df25k, que contendrá sólo las filas en las que la columna ACRES contenga un valor superior a 25000.

df25k = filter(df, ACRES >= 25000)

• En el siguiente ejemplo de código, se pasan dos expresiones condicionales. La primera se utiliza para filtrar los registros de manera que el número de ACRES sea mayor mayor o igual a 25.000, y la segunda filtra los registros para los que la columna YEAR contenga un valor de 2016.

df25k = filter(df, ACRES >= 25000, YEAR == 2016)

• En este caso, la variable df25k incluirá los registros en los que coincidan ambas condiciones: la superficie quemada es superior a 25000 y el año del incendio fue 2016. Este también puede reescribirse como un único parámetro que utiliza el operador & para combinar como se ve a continuación.

df25k = filter(df, ACRES >= 25000 & YEAR == 2016)

• En este ejercicio aprenderá a utilizar la función filter() para crear un subconjunto de registros basado en algún valor. Los ejercicios de este capítulo requieren los siguientes paquetes: readr, dplyr, ggplot2. Pueden cargarse desde el panel de paquetes, el panel de la consola o un script. Si es necesario, establezca el directorio de trabajo escribiendo el código: setwd().

• El repositorio RDataSets contiene un archivo llamado StudyArea.csv, que es un archivo separado por comas que contiene datos de incendios forestales de los años 1980-2016 para los estados de California, Oregón, Washington, Idaho, Montana, Wyoming, Colorado, Utah, Nevada, Arizona y Nuevo México. Hay algo más de 439.000 registros en este archivo y hay 37 columnas de información que describen cada incendio durante este periodo.

• Utilice la función read_csv() para cargar el conjunto de datos en un dataframe.

dfFires = read_csv("C:/Data/StudyArea.csv", col_types = list(UNIT = col_character()), col_names = TRUE)

• Utilice la función nrow() para asegurarse de que los aproximadamente 439.000 registros se hayan cargado.

nrow(dfFires)

## [1] 439362

• Inicialmente utilizaremos una única expresión condicional con la función filter() para crear un subconjunto de registros que contenga sólo los incendios forestales de más de 25.000 acres.

df25k = filter(dfFires, TOTALACRES >= 25000)

• Obtenga un recuento del número de registros que coinciden con el filtro. Debería haber 655 filas. También puede utilizar la función View(df25k) para ver los datos en un formato tabular.

nrow(df25k)

## [1] 655

head(df25k)

## # A tibble: 6 × 14
##     FID ORGANIZATI UNIT  SUBUNIT SUBUNIT2        FIRENAME CAUSE YEAR_ STARTDATED
##   <dbl> <chr>      <chr> <chr>   <chr>           <chr>    <chr> <dbl> <chr>     
## 1  1875 FWS        22570 USAZKGR Kofa National … KING VA… Human  2005 10/1/05 0…
## 2  1968 FWS        81674 USCAHPR Hopper Mountai… PIRU     Human  2003 10/23/03 …
## 3  2476 FWS        81720 USCATNR San Diego Nati… CEDAR    Human  2003 10/25/03 …
## 4  2477 FWS        81720 USCATNR San Diego Nati… MINE     Human  2003 10/26/03 …
## 5  6929 FWS        61520 USMTCMR Charles M. Rus… SODA CR… Natu…  2006 7/16/06 0…
## 6  9161 FWS        14621 USORSLR Sheldon Nation… BADGER … Natu…  1999 8/24/99 0…
## # ℹ 5 more variables: CONTRDATED <chr>, OUTDATED <chr>, STATE <chr>,
## #   STATE_FIPS <dbl>, TOTALACRES <dbl>

• También puede incluir múltiples expresiones condicionales como parte del filtro.

df1k = filter(dfFires, TOTALACRES >= 1000, YEAR_ == 2016)
nrow(df1k)

## [1] 152

• También puede combinar las expresiones en una sola expresión con varias condiciones, como se ve a continuación. El carácter & es el operador “y”. Tendrías que usar el carácter | para incluir un operador “o”.

df1k = filter(dfFires, TOTALACRES >= 1000 & YEAR_ == 2016)

• Por último, cuando se tiene una lista de valores potenciales que se desea incluir por el filtro se puede utilizar la sentencia %in%. Esta línea de código en particular crearía un dataframe que contiene los incendios ocurridos en los años 2010, 2011 o 2012.

dfYear = filter(dfFires, YEAR_ %in% c(2010, 2011, 2012))
head(dfYear)

## # A tibble: 6 × 14
##     FID ORGANIZATI UNIT  SUBUNIT SUBUNIT2        FIRENAME CAUSE YEAR_ STARTDATED
##   <dbl> <chr>      <chr> <chr>   <chr>           <chr>    <chr> <dbl> <chr>     
## 1 12831 FWS        81627 USCASRR Sacramento Riv… SWR-SR-… Human  2010 5/10/10 0…
## 2 12839 FWS        84555 USNVDSR Desert Nationa… MORMON … Human  2010 6/28/10 0…
## 3 12840 FWS        13561 USWALPR Little Pend Or… WINSLOW… Natu…  2010 7/31/10 0…
## 4 12841 FWS        13700 USWAHFR Hanford Reach … SADDLE … Natu…  2010 7/31/10 0…
## 5 12847 FWS        84555 USNVDSR Desert Nationa… DEAD HO… Natu…  2010 7/25/10 0…
## 6 12921 FWS        22530 USAZBAR Buenos Aires N… M M 14.… Unde…  2010 11/2/10 0…
## # ℹ 5 more variables: CONTRDATED <chr>, OUTDATED <chr>, STATE <chr>,
## #   STATE_FIPS <dbl>, TOTALACRES <dbl>

Acotar la lista de columnas con select()

• El archivo StudyArea.csv con el que ha estado trabajando en los ejercicios incluye 37 columnas de información. En la mayoría de los casos no necesitará todas las columnas. La función select() puede utilizarse para reducir la lista de columnas para incluir sólo las necesarias para una tarea.

• Utilice la función read_csv() para cargar el conjunto de datos en un dataframe. Nota: En aras de la exhaustividad, se cargarán los datos externos del archivo StudyArea.csv en el dataframe dfFires, pero este paso no es absolutamente necesario si está haciendo los ejercicios en secuencia en la misma sesión de R Studio.

dfFires = read_csv("C:/Data/StudyArea.csv", col_types = list(UNIT = col_character()), col_names = TRUE)

• En una nueva línea, añada una llamada a la función select() como se ve a continuación para limitar las columnas que se devuelven.

dfFires2 = select(dfFires, FIRENAME, TOTALACRES, YEAR_)

• Visualice las primeras filas y observe que ahora sólo tenemos tres columnas.

head(dfFires2)

## # A tibble: 6 × 3
##   FIRENAME   TOTALACRES YEAR_
##   <chr>           <dbl> <dbl>
## 1 PUMP HOUSE        0.1  2001
## 2 I5                3    2002
## 3 SOUTHBAY          0.5  2002
## 4 MARINA            0.1  2001
## 5 HILL              1    1994
## 6 IRRIGATION        0.1  1994

• Muchos de los nombres de las columnas que importes no serán muy fáciles de leer, por lo que no es raro que quiera renombrar las columnas también. Esto puede lograrse usando la función select().

dfFires2 = select(dfFires, "FIRE" = "FIRENAME", "ACRES" ="TOTALACRES", "YR" = "YEAR_")
head(dfFires2)

## # A tibble: 6 × 3
##   FIRE       ACRES    YR
##   <chr>      <dbl> <dbl>
## 1 PUMP HOUSE   0.1  2001
## 2 I5           3    2002
## 3 SOUTHBAY     0.5  2002
## 4 MARINA       0.1  2001
## 5 HILL         1    1994
## 6 IRRIGATION   0.1  1994

• También hay una serie de útiles funciones de ayuda que puede utilizar con la función select() para filtrar las columnas devueltas. Estas incluyen starts_with(), ends_with(), contains(), matches() y num_range().

dfFires3 = select(dfFires, contains("DATE"))
head(dfFires3)

## # A tibble: 6 × 3
##   STARTDATED   CONTRDATED   OUTDATED
##   <chr>        <chr>        <chr>   
## 1 1/1/01 0:00  1/1/01 0:00  <NA>    
## 2 5/3/02 0:00  5/3/02 0:00  <NA>    
## 3 6/1/02 0:00  6/1/02 0:00  <NA>    
## 4 7/12/01 0:00 7/12/01 0:00 <NA>    
## 5 9/13/94 0:00 9/13/94 0:00 <NA>    
## 6 4/22/94 0:00 4/22/94 0:00 <NA>

• También puede hacer varias llamadas a estas funciones de ayuda

dfFires3 = select(dfFires, contains("DATE"), starts_with("TOTAL"))
head(dfFires3)

## # A tibble: 6 × 4
##   STARTDATED   CONTRDATED   OUTDATED TOTALACRES
##   <chr>        <chr>        <chr>         <dbl>
## 1 1/1/01 0:00  1/1/01 0:00  <NA>            0.1
## 2 5/3/02 0:00  5/3/02 0:00  <NA>            3  
## 3 6/1/02 0:00  6/1/02 0:00  <NA>            0.5
## 4 7/12/01 0:00 7/12/01 0:00 <NA>            0.1
## 5 9/13/94 0:00 9/13/94 0:00 <NA>            1  
## 6 4/22/94 0:00 4/22/94 0:00 <NA>            0.1

Organizar las filas

• La función arrange() del paquete dplyr puede utilizarse para ordenar las filas de un dataframe. Esta función acepta un conjunto de columnas para ordenarlas, siendo el orden por defecto ordenación de las filas en orden ascendente. Sin embargo, puede pasar la función desc() helper para ordenar las filas en orden descendente.

• Utilice la función read_csv() para cargar el conjunto de datos en un dataframe.

dfFires = read_csv("C:/Data/StudyArea.csv", col_types = list(UNIT = col_character()), col_names = TRUE)

• Filtrar el conjunto de datos para que contenga solo los incendios de más de 1.000 acres quemados del año 2016.

df1k = filter(dfFires, TOTALACRES >= 1000, YEAR_ == 2016)

• Añade y ejecuta el código que ves a continuación para crear un subconjunto de columnas y renombrarlas

df1k = select(df1k, "NAME" = "FIRENAME", "ACRES" = "TOTALACRES", "YR" = "YEAR_")

• Ordena las filas para que estén en orden ascendente

arrange(df1k, ACRES)

## # A tibble: 152 × 3
##    NAME        ACRES    YR
##    <chr>       <dbl> <dbl>
##  1 Crackerbox  1000   2016
##  2 Lakes       1000   2016
##  3 Choulic 2   1008   2016
##  4 Amigo Wash  1020.  2016
##  5 Granite     1030   2016
##  6 Tie         1031   2016
##  7 Black       1040   2016
##  8 Bybee Creek 1072   2016
##  9 MARSHES     1080   2016
## 10 Bug Creek   1089   2016
## # ℹ 142 more rows

• Utilice la función de ayuda desc() para ordenar las filas en orden descendente.

arrange(df1k, desc(ACRES))

## # A tibble: 152 × 3
##    NAME        ACRES    YR
##    <chr>       <dbl> <dbl>
##  1 PIONEER    188404  2016
##  2 Junkins    181320  2016
##  3 Range 12   171915  2016
##  4 Erskine     48007  2016
##  5 Cedar       45977  2016
##  6 Maple       45425  2016
##  7 Rail        43799  2016
##  8 North Fire  42102  2016
##  9 Laidlaw     39813  2016
## 10 BLUE CUT    36274  2016
## # ℹ 142 more rows

Añadir filas con mutate()

• La función mutate() se utiliza para añadir nuevas columnas a un dataframe que son resultado de operaciones entre otras columnas del dataframe. Las nuevas columnas creadas con la función mutate() se añadirán al final del dataframe. Para este ejercicio usará el paquete lubridate. El paquete lubridate forma parte de tidyverse y se utiliza para trabajar con fechas y horas.

• Utilice la función read_csv() para cargar el conjunto de datos en un dataframe

dfFires = read_csv("C:/Data/StudyArea.csv", 
                   col_types = list(UNIT = col_character()),
                   col_names = TRUE)

• Utilice la función select() para definir un conjunto de columnas para el dataframe.

df = select(dfFires, ORGANIZATI, STATE, YEAR_, TOTALACRES, CAUSE, STARTDATED)

• Realice un filtrado básico de los datos para que sólo se incluyan los incendios de más de 1.000 acres quemados y que tengan una causa humana o natural

df = filter(df, TOTALACRES >= 1000 & CAUSE %in% c("Human", "Natural"))

• Utilice la función mutate() para crear una nueva columna DOY que contenga el día del año en que se inició el incendio. La función yday() del paquete lubridate se utiliza para devolver el día del año utilizando una fecha de entrada formateada de la columna STARTDATED.

df = mutate(df, DOY = yday(as.Date(df$STARTDATED, format="%m/%d/%y%H:%M")))
knitr::kable(head(df))

Resumir y agrupar

• Las estadísticas de resumen de un dataframe pueden producirse con la función summarize(). La agrupación de datos en un dataframe facilita el paradigma de dividir-aplicar-combinar. Este paradigma primero divide los datos en grupos, utilizando la función group_by() en dplyr, luego aplica el análisis al grupo y, finalmente, combina los resultados.

• En este ejercicio utilizará las funciones mutate(), summarize() y group_by() para agrupar los incendios forestales por década y producir un resumen del tamaño medio de los incendios forestales para cada década.

• Utilice la función read_csv() para cargar el conjunto de datos en un dataframe.

dfFires = read_csv("C:/Data/StudyArea.csv", 
                   col_types = list(UNIT = col_character()), 
                   col_names = TRUE)

• Seleccione las columnas que se utilizarán en el ejercicio.

df = select(dfFires, ORGANIZATI, STATE, YEAR_, TOTALACRES, CAUSE)

• Filtrar los registros

df = filter(df, TOTALACRES >= 1000)

• Utilice la función mutate() para crear una nueva columna llamada DECADE que define la década en la que se produjo cada incendio. En este caso se llama a la función ifelse() para producir los valores de cada década.

df <- mutate(df, DECADE = ifelse(YEAR_ %in% 1980:1989, "1980-1989", 
                          ifelse(YEAR_ %in% 1990:1999, "1990-1999",
                          ifelse(YEAR_ %in% 2000:2009, "2000-2009", 
                          ifelse(YEAR_ %in% 2010:2016, "2010-2016", "-99")))))
knitr::kable(head(df))

• Utilice la función group_by() para agrupar el dataframe por década.

grp = group_by(df, DECADE)

• Resuma el tamaño medio de los incendios forestales por década utilizando la función summarize()

sm = summarize(grp, mean(TOTALACRES))
knitr::kable(sm)

• Pongamos en orden las cosas cambiando el nombre de la nueva columna producida por la función summarize()

names(sm) <- c("DECADE", "MEAN_ACRES_BURNED")

• Por último, vamos a crear un gráfico de barras con los resultados. Discutiremos la creación de diferentes tipos de tablas y gráficos a medida que avancemos en los capítulos posteriores

library(ggplot2)

ggplot(data=sm) + 
  geom_col(mapping = aes(x=DECADE, y=MEAN_ACRES_BURNED), fill="red")

Tuberías

• La canalización es una forma eficiente de enviar la salida de una función a otra función sin crear un conjunto de datos intermedio y es más útil cuando se tiene una serie de funciones para ejecutar. La sintaxis para la canalización consiste en utilizar los caracteres %>% al final de cada sentencia que desee canalizar. En este ejercicio aprenderá a utilizar la canalización para encadenar dataframe de entrada y salida. Resolveremos el ejercicio anterior usando las tuberias.

library(lubridate)

df = read_csv("C:/Data/StudyArea.csv", col_types = list(UNIT = col_character()), col_names = TRUE) %>%
  select(ORGANIZATI, STATE, YEAR_, TOTALACRES, CAUSE, STARTDATED) %>%
  filter(TOTALACRES >= 1000 & CAUSE %in% c("Human", "Natural")) %>%
  mutate(DOY = yday(as.Date(STARTDATED, format="%m/%d/%y %H:%M")))

knitr::kable(head(df))

• La primera línea de código lee el contenido del archivo externo StudyArea.csv en una variable de dataframe df como hemos hecho en todos los demás ejercicios de este capítulo. Sin embargo, observará la inclusión de la sentencia de canalización (%>%>) al final de la línea. Esto asegura que el contenido de la variable df será automáticamente automáticamente a la función select().

• Observe que la función select() no crea una variable como hemos hecho en ejercicios anteriores, y que hemos omitido el primer parámetro, que normalmente habría sido la variable del dataframe. Se da a entender que la variable df será pasada a la función select().

• Este mismo proceso de incluir la sentencia piping al final de cada línea y dejando fuera el primer parámetro se repite para todas las líneas de código adicionales en las que queramos pasar automáticamente la variable df a la siguiente función de dplyr.

• En esta sección se utilizó el paquete dplyr para realizar varias funciones de transformación de datos. Se aprendió a limitar columnas con la función select() filtrar un marco de datos basado en una o más expresiones, añadir columnas con mutate(), y a resumir y agrupar datos. Por último, ha aprendido a utilizar la canalización para hacer su código más eficiente. En el próximo capítulo se aprenderá a crear conjuntos de datos ordenados con el paquete tidyr.

Ejercicio para entregar

• Crear un nuevo dataframe que sea un subconjunto del dataframe original de dfFires. El subconjunto debe contener todos los incendios del Estado de Idaho y las columnas deben ser limitadas para que sólo estén presentes las columnas YEAR_, CAUSE y TOTALACRES. Cambie el nombre de las columnas. Agrupe los datos por CAUSE y YEAR_ y luego resuma por el total de acres quemados. Trazar los resultados.

• Resuelva la Tarea 1.1 de la Sección 1 de Python utilizando R.

Recopilación

• Un problema común en muchos conjuntos de datos es que los nombres de las columnas no son variables sino valores de una variable. En la figura siguiente, las columnas 1999 y 2000 son en realidad valores de la variable YEAR. Cada fila de la tabla existente representa en realidad dos observaciones. El paquete tidyr puede utilizarse para reunir estas columnas existentes en una nueva variable. En este caso, tenemos que crear una nueva columna columna llamada YEAR y luego reunir en esta los valores existentes en las columnas 1999 y 2000.

• La función gather() del paquete tidyr puede utilizarse para realizar la recopilación de datos.

gather("1999", "2000", key = "year", value = "cases")

• Hay tres parámetros de la función gather(). El primero es el conjunto de columnas que representan lo que deberían ser valores y no variables. Estas serían las columnas 1999 y 2000 en el ejemplo que hemos seguido.
• A continuación, tendrá que nombrar la variable de la nueva columna. Esto también se llama la clave key, y en este caso será la variable del año. Por último, tendrás que proporcionar el valor value, que es el nombre de la variable cuyos valores se reparten por las celdas.

• En este ejercicio aprenderás a utilizar la función gather() para resolver este tipo de problemas:

1. Descargue el archivo CountryPopulation.csv localizado en RDataSets. Abra este archivo, por ejemplo en Microsoft Excel, o algún otro tipo de software de hoja de cálculo. El archivo debería tener un aspecto similar al de la pantalla que se muestra a continuación. Esta hoja de cálculo incluye el año 2017.

• Las columnas de cada año representan valores, no variables. Estas columnas deben reunirse en un nuevo par de variables que representen el año y la población. En este ejercicio utilizarás la función gather() para realizar esta tarea de ordenación de datos.

2. Si es necesario, cargue los paquetes readr y tidyr haciendo clic en las casillas de verificación en el panel de paquetes o incluyendo la siguiente línea de código.

library(readr)
library(tidyr)

3. Cargue el archivo CountryPopulation.csv en RStudio escribiendo el código que ve abajo en el panel de la consola

dfPop = read_csv("C:/Data/CountryPopulation.csv", col_names = TRUE)
str(dfPop)

## spc_tbl_ [264 × 10] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ Country Name: chr [1:264] "Aruba" "Afghanistan" "Angola" "Albania" ...
##  $ Country Code: chr [1:264] "ABW" "AFG" "AGO" "ALB" ...
##  $ 2010        : num [1:264] 101669 28803167 23369131 2913021 84449 ...
##  $ 2011        : num [1:264] 102053 29708599 24218565 2905195 83751 ...
##  $ 2012        : num [1:264] 102577 30696958 25096150 2900401 82431 ...
##  $ 2013        : num [1:264] 103187 31731688 25998340 2895092 80788 ...
##  $ 2014        : num [1:264] 103795 32758020 26920466 2889104 79223 ...
##  $ 2015        : num [1:264] 104341 33736494 27859305 2880703 78014 ...
##  $ 2016        : num [1:264] 104822 34656032 28813463 2876101 77281 ...
##  $ 2017        : num [1:264] 105264 35530081 29784193 2873457 76965 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   `Country Name` = col_character(),
##   ..   `Country Code` = col_character(),
##   ..   `2010` = col_double(),
##   ..   `2011` = col_double(),
##   ..   `2012` = col_double(),
##   ..   `2013` = col_double(),
##   ..   `2014` = col_double(),
##   ..   `2015` = col_double(),
##   ..   `2016` = col_double(),
##   ..   `2017` = col_double()
##   .. )
##  - attr(*, "problems")=<externalptr>

4. Utilice la función View() para mostrar los datos en una estructura tabular.

View(dfPop)

5. Utilice la función gather() como se ve a continuación

dfPop2 = gather(dfPop, 
                "2010", "2011", "2012", "2013", "2014", "2015", "2016", "2017", 
                key = "YEAR", 
                value = "POPULATION")
knitr::kable(head(dfPop2, 10))

• Otra opción también sería la siguiente, cuando contamos con un gran número de columnas de este tipo

years <- colnames(dfPop)[grep("^\\d{4}$", colnames(dfPop))]

dfPop2 <- dfPop %>%
  gather(key = "YEAR", value = "POPULATION", all_of(years))

knitr::kable(head(dfPop2, 10))

• Donde, en expresiones regulares, "^\\d{4}$" tiene el siguiente significado:

  • ^: Representa el inicio de una cadena.
  • \\d: Representa un dígito.
  • {4}: Indica que el elemento anterior (\\d, en este caso) debe aparecer exactamente 4 veces.
  • $: Representa el final de una cadena.

• En resumen, "^\\d{4}$" coincide con cualquier cadena que contenga exactamente 4 dígitos y no contenga ningún otro carácter adicional antes o después de los dígitos.

Distribución

• La distribución es lo contrario de la reunión y se utiliza cuando una observación se distribuye en varias filas. En el diagrama siguiente, la tabla debe definir una observación de un país por año. Sin embargo, observará que está repartida en dos filas. Una fila para cases y otra para population.

• Podemos utilizar la función spread() para solucionar este problema. La función spread() toma dos parámetros: la columna que contiene los nombres de las variables repetidas, conocida como key y una columna que contiene valores de múltiples variables, que llamamos value

• Instale el paquete devtools y los conjuntos de datos DSR utilizando el código que ve a continuación escribiendo en el panel de la consola. Alternativamente, puede utilizar el panel de paquetes para instalar los librerías

install.packages("devtools")
devtools::install_github("garrettgman/DSR")

knitr::kable(head(table2))

• Utilice la función spread() para corregir este problema.

table2b = spread(table2, key = type, value = count)
knitr::kable(head(table2b))

Separación

• Otro caso común es el de dos variables que se colocan en la misma columna. Por ejemplo, la hoja de cálculo siguiente tiene una columna State-County Name que en realidad contiene dos variables separadas por una barra.

• La función separate() puede utilizarse para dividir una columna en varias columnas dividiendo por un separador. Por defecto, la función separate() buscará automáticamente cualquier carácter no alfanumérico o se puede definir un carácter específico (ver separate()).

1. En la carpeta de RDataSets hay un archivo llamado usco2005.csv. Abra este archivo, por ejemplo con Microsoft Excel, o algún otro tipo de software de hoja de cálculo. El archivo debería tener un aspecto similar al de la captura de pantalla de abajo.

2. Cargue el archivo usco2005.csv en RStudio escribiendo el código que ve a continuación en el panel de la consola

df = read_csv("C:/Data/usco2005.csv", col_names = TRUE)
knitr::kable(head(df, 10))

3. Utilice la función separate() para separar el contenido de la columna State-County Name en las columnas StateAbbrev y CountyName

df2 = separate(df, "State-County Name",into = c("StateAbbrev", "CountyName"))
knitr::kable(head(df2, 10))

Unión

• La función unite() es exactamente lo contrario de separate(), ya que combina varias columnas en una sola. Aunque no se utiliza tan a menudo como separate(), puede haber ocasiones en las que necesite la funcionalidad proporcionada por unite(). En este ejercicio unirás el marco de datos que fue separado en el último ejercicio.

• En el panel de la consola, añada el código que ve a continuación para unir las columnas StateAbbrev y CountyName en una sola columna. Por defecto unite() realiza unión con base en la expresión regular "_", "unite(data, col, ..., sep = "_", remove = TRUE)"

df3 = unite(df2, State_County_Name, StateAbbrev, CountyName)
knitr::kable(head(df3, 10))

Conclusión

• En este capítulo se ha presentado el paquete tidyr y su conjunto de funciones para crear conjuntos de datos ordenados. El próximo capítulo le enseñará los fundamentos de la exploración de datos datos utilizando R y tidyverse.

Ejercicio para entregar

• Considere el conjunto de datos us_state_population.tsv utilizado en la sección de Python, para la creación del mapa coroplético de Estados Unidos. Repita el procedimiento planteado en cada ítem de esta sección para obtener el nuevo dataframe con las nuevas columnas Year y Population. Realice unión y separación utilizando las columnas State y Code.

Técnicas básicas de exploración

• El Análisis Exploratorio de Datos (EDA) es un flujo de trabajo diseñado para obtener una mejor comprensión de sus datos. El flujo de trabajo consta de tres pasos.

  1. El primero consiste en generar preguntas sobre los datos. En este paso se quiere ser lo más amplio posible, ya que en este momento no se tiene una buena idea de los datos.
  2. A continuación, se buscan las respuestas a estas preguntas mediante la visualización, la transformación y el modelado de los datos.
  3. Por último, refine sus preguntas o genere nuevas preguntas.

• En R hay dos herramientas principales que apoyan el proceso de exploración de datos: los gráficos y los estadísticos de resumen. Los datos pueden dividirse generalmente en tipos categóricos o continuos.

• Las variables categóricas consisten en un pequeño conjunto de valores, mientras que las continuas tienen un un conjunto potencialmente infinito de valores ordenados.

• Las variables categóricas suelen visualizarse con gráficos de barras, y las continuas con histogramas. Tanto datos categóricos y continuos pueden representarse a través de varios gráficos creados con R.

• Al realizar la visualización básica de las variables, tendemos a medir la variación o la covariación. La variación es la tendencia de los valores de una variable a cambiar de una medición a otra. La covariación es la tendencia de los valores de dos o más variables a variar de una manera relacionada.

  • Medir la variación categórica con un gráfico de barras
  • Medición de la variación continua con un histograma
  • Medición de la covariación con gráficos de caja
  • Medición de la covariación con el tamaño de los símbolos
  • Creación de gráficos 2D y hexadecimales
  • Generación de estadísticas de resumen

Medición de la variación categórica con un gráfico de barras

• Un gráfico de barras es una buena manera de visualizar datos categóricos. Se separa cada categoría en una barra separada y la altura de cada barra se define por el número de ocurrencias en esa categoría.

1. Los ejercicios de esta sección requieren los siguientes paquetes: readr, dplyr, ggplot2. Pueden cargarse desde el panel de paquetes, el panel de la consola o un script

2. Utilice la función read_csv() para cargar el conjunto de datos en un marco dataframe.

library(dplyr)
library(readr)
library(ggplot2)

dfFires <- read_csv("C:/Data/StudyArea.csv", col_types = list(UNIT = col_character()), col_names = TRUE)
str(dfFires)

## spc_tbl_ [439,362 × 14] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ FID       : num [1:439362] 0 1 2 3 4 5 6 18 20 21 ...
##  $ ORGANIZATI: chr [1:439362] "FWS" "FWS" "FWS" "FWS" ...
##  $ UNIT      : chr [1:439362] "81682" "81682" "81682" "81682" ...
##  $ SUBUNIT   : chr [1:439362] "USCADBR" "USCADBR" "USCADBR" "USCADBR" ...
##  $ SUBUNIT2  : chr [1:439362] "San Diego Bay National Wildlife Refuge" "San Diego Bay National Wildlife Refuge" "San Diego Bay National Wildlife Refuge" "San Diego Bay National Wildlife Refuge" ...
##  $ FIRENAME  : chr [1:439362] "PUMP HOUSE" "I5" "SOUTHBAY" "MARINA" ...
##  $ CAUSE     : chr [1:439362] "Human" "Human" "Human" "Human" ...
##  $ YEAR_     : num [1:439362] 2001 2002 2002 2001 1994 ...
##  $ STARTDATED: chr [1:439362] "1/1/01 0:00" "5/3/02 0:00" "6/1/02 0:00" "7/12/01 0:00" ...
##  $ CONTRDATED: chr [1:439362] "1/1/01 0:00" "5/3/02 0:00" "6/1/02 0:00" "7/12/01 0:00" ...
##  $ OUTDATED  : chr [1:439362] NA NA NA NA ...
##  $ STATE     : chr [1:439362] "California" "California" "California" "California" ...
##  $ STATE_FIPS: num [1:439362] 6 6 6 6 6 6 6 6 6 6 ...
##  $ TOTALACRES: num [1:439362] 0.1 3 0.5 0.1 1 0.1 3 0.1 0.1 0.1 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   FID = col_double(),
##   ..   ORGANIZATI = col_character(),
##   ..   UNIT = col_character(),
##   ..   SUBUNIT = col_character(),
##   ..   SUBUNIT2 = col_character(),
##   ..   FIRENAME = col_character(),
##   ..   CAUSE = col_character(),
##   ..   YEAR_ = col_double(),
##   ..   STARTDATED = col_character(),
##   ..   CONTRDATED = col_character(),
##   ..   OUTDATED = col_character(),
##   ..   STATE = col_character(),
##   ..   STATE_FIPS = col_double(),
##   ..   TOTALACRES = col_double()
##   .. )
##  - attr(*, "problems")=<externalptr>

3. Para este análisis, filtraremos los datos de modo que solo estén representados los incendios que quemaron más de 1.000 acres (4046.86 metros cuadrados o 0.4047 hectáreas) en los años 2010 a 2016. Añada el código que que se ve a continuación para filtrar los datos y enviar los resultados a un gráfico de barras.

df <- filter(dfFires, TOTALACRES >= 1000, 
             YEAR_ %in% c(2010, 2011, 2012, 2013, 2014, 2015, 2016))
ggplot(data = df) + 
  geom_bar(mapping = aes(x = YEAR_)) +
  scale_x_continuous(breaks = c(2010, 2011, 2012, 2013, 2014, 2015, 2016)) +
  labs(x = "Year", y = "Count")

4. Utilice la función count() para obtener el recuento real de cada categoría.

count(df, YEAR_)

## # A tibble: 7 × 2
##   YEAR_     n
##   <dbl> <int>
## 1  2010   127
## 2  2011   240
## 3  2012   352
## 4  2013   170
## 5  2014   131
## 6  2015   205
## 7  2016   152

Medición de la variación continua con un histograma

• La distribución de una variable continua se puede medir con el uso de un histograma. En este ejercicio crearás un histograma de los acres quemados en incendios forestales

1. En una nueva línea, utilice la función read_csv() para cargar el archivo StudyArea.csv

dfFires <- read_csv("C:/Data/StudyArea.csv", col_types = list(UNIT = col_character()), col_names = TRUE)

2. Introduzca el marco de datos y utilice la función select() para limitar las columnas y filtrar las filas para que sólo se incluyan los incendios superiores a 1.000 acres. Dado que tenemos un gran número de incendios forestales que quemaron sólo un pequeño número de acres, nos centraremos en los incendios que son un poco más grandes en este caso.

df %>%
  select(ORGANIZATI, STATE, YEAR_, TOTALACRES, CAUSE) %>%
  filter(TOTALACRES >= 1000)

## # A tibble: 1,377 × 5
##    ORGANIZATI STATE      YEAR_ TOTALACRES CAUSE       
##    <chr>      <chr>      <dbl>      <dbl> <chr>       
##  1 FWS        Washington  2010       2755 Natural     
##  2 FWS        Idaho       2010       1801 Undetermined
##  3 FWS        Montana     2010      27898 Natural     
##  4 FWS        Montana     2012       9468 Natural     
##  5 FWS        Washington  2012       1141 Human       
##  6 FWS        Oregon      2010       3083 Human       
##  7 FWS        Oregon      2013      51015 Undetermined
##  8 FWS        Washington  2016     171915 Human       
##  9 FWS        Arizona     2015       6000 Natural     
## 10 FWS        Arizona     2016       1624 Natural     
## # ℹ 1,367 more rows

3. Cree el histograma utilizando ggplot() con geom_hist() y un tamaño de la celda de 500. Los datos están claramente sesgados hacia el extremo inferior del número de acres quemados.

df %>%
  select(ORGANIZATI, STATE, YEAR_, TOTALACRES, CAUSE) %>%
  filter(TOTALACRES >= 1000) %>%
  slice(1:300) %>%
  ggplot() + 
  geom_histogram(mapping = aes(x=TOTALACRES), binwidth=1000) +
  labs(x = "Total Acres", y = "Count")

4. También se puede obtener un recuento físico del número de incendios que han caído en cada casilla. Al ver el histograma y el recuento, se puede oservar que la gran mayoría de incendios son pequeños.

df %>%
  count(cut_width(TOTALACRES, 100))

## # A tibble: 338 × 2
##    `cut_width(TOTALACRES, 100)`     n
##    <fct>                        <int>
##  1 [950,1050]                      43
##  2 (1050,1150]                     73
##  3 (1150,1250]                     38
##  4 (1250,1350]                     34
##  5 (1350,1450]                     39
##  6 (1450,1550]                     36
##  7 (1550,1650]                     26
##  8 (1650,1750]                     43
##  9 (1750,1850]                     37
## 10 (1850,1950]                     24
## # ℹ 328 more rows

5. Reconstruya el histograma utilizando un tamaño de celda de 5000. ¿Cuál es el efecto en la salida?

Medición de la covariación con gráficos de caja

• Los gráficos de caja proporcionan una representación visual de la dispersión de los datos de una variable. Estos gráficos muestran el rango de valores de una variable junto con la mediana y los cuartiles. Siga las instrucciones proporcionadas a continuación para crear un gráfico de caja que mide la covariación entre la organización y la superficie total quemada.

1. Utilice la función read_csv() para cargar el archivo StudyArea.csv en un dataframe

dfFires <- read_csv("C:/Data/StudyArea.csv", 
                    col_types = list(UNIT = col_character()), 
                    col_names = TRUE)

2. Extraiga el marco de datos y filtre las filas para que sólo se incluyan los incendios de entre 5000 y 10000 acres. A continuación, agrupe los datos por organización. La columna ORGANIZATI del conjunto de datos contiene datos categóricos de los organismos del gobierno federal de EE.UU. que han tenido tierras afectadas por incendios forestales. Por último utilice ggplot() con geom_boxplot() para crear un boxplot que muestre la distribución de incendios forestales por organización

df %>%
  filter(TOTALACRES >= 5000 & TOTALACRES <= 10000) %>%
  group_by(ORGANIZATI) %>%
  ggplot(mapping = aes(x = ORGANIZATI, y = TOTALACRES)) + geom_boxplot() +
  labs(x = "Organization", y = "Total Acres")

3. Puede crear también un nuevo boxplot que mapee la covariación de CAUSE y TOTALACRES

Medición de la covariación con el tamaño del símbolo

• La función geom_count() puede utilizarse con ggplot() para medir la covariación entre variables utilizando diferentes tamaños de símbolos. Siga las instrucciones proporcionadas a continuación para medir la covariación entre la organización y la causa del incendio forestal utilizando el tamaño del símbolo

1. Use la función read_csv() para cargar el conjunto de datos en un dataframe

dfFires <- read_csv("C:/Data/StudyArea.csv", 
                    col_types = list(UNIT = col_character()), 
                    col_names = TRUE)

2. Extraiga el marco de datos y filtre las filas para que sólo se incluyan los incendios forestales que se originaron por causas naturales o humanas. Esto eliminará cualquier registro que son desconocidos o tienen valores perdidos. A continuación, utilice geom_count() para crear un gráfico de símbolos graduados basado en el número de incendios por organización.

df %>%
  filter(CAUSE == "Natural" | CAUSE == "Human") %>%
  group_by(ORGANIZATI) %>%
  ggplot() + geom_count(mapping = aes(x = ORGANIZATI, y = CAUSE)) +
  labs(x = "Organization", y = "Cause")

3. También puede obtener un recuento exacto del número de incendios por organización y causa

df %>%
  count(ORGANIZATI, CAUSE)

## # A tibble: 11 × 3
##    ORGANIZATI CAUSE            n
##    <chr>      <chr>        <int>
##  1 BIA        Human           49
##  2 BIA        Natural         91
##  3 BLM        Human          187
##  4 BLM        Natural        386
##  5 FS         Human          158
##  6 FS         Natural        431
##  7 FWS        Human           10
##  8 FWS        Natural          7
##  9 FWS        Undetermined     6
## 10 NPS        Human            6
## 11 NPS        Natural         46

Gráficos 2D de contenedores y hexágonos

• También puede utilizar gráficos 2D de contenedores y hexadecimales como una forma alternativa de ver la distribución de dos variables. Siga las instrucciones que se proporcionan a continuación para crear gráficos 2D de cubos y hexágonos que visualicen la relación entre el año y la superficie total quemada

1. Utilice la función read_csv() para cargar el conjunto de datos en un marco de datos

dfFires <- read_csv("C:/Data/StudyArea.csv", 
                    col_types = list(UNIT = col_character()),
                    col_names = TRUE)

2. Cree un mapa 2D con YEAR_ en el eje \(x\) y TOTALACRES en el eje \(y\)

ggplot(data = dfFires) + 
  geom_bin2d(mapping = aes(x=YEAR_, y=TOTALACRES)) +
  labs(x = "Year", y = "Total Acres")

Generación de estadísticas de resumen

• Otra técnica básica para realizar un análisis exploratorio de datos es generar varias estadísticas de resumen sobre un conjunto de datos. R incluye una serie de funciones individuales para generar estadísticas de resumen específicas o puede utilizar la función summary() para generar un conjunto de estadísticas de resumen

1. Recargar el archivo StudyArea.csv en un marco de datos

df <- read_csv("C:/Data/StudyArea.csv", 
               col_types = list(UNIT = col_character()),
               col_names = TRUE)

2. Restringir la lista de columnas

df <- select(df, ORGANIZATI, STATE, YEAR_, TOTALACRES, CAUSE)

3. Filtra la lista para incluir sólo los incendios forestales de más de 1000 acres

df <- filter(df, TOTALACRES >= 1000)

4. Llame a la función mean(), pasando una referencia al marco de datos y a la columna columna TOTALACRES.

mean(df$TOTALACRES)

## [1] 10813.06

5. Llama a la función median()

median(df$TOTALACRES)

## [1] 3240

6. En lugar de llamar a las funciones individuales de estadísticas de resumen, puede simplemente utilizar la función summary() para devolver una lista de estadísticas de resumen.

summary(df$TOTALACRES)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1000    1670    3240   10813    8282  590620

Conclusiones

• En esta sección ha aprendido algunas técnicas básicas de exploración de datos utilizando R. Ha aprendido a medir la variación categórica y continua con gráficos de barras e histogramas, y la covariación con gráficos de caja y diferentes tamaños de símbolos. Por último, finalmente, aprendió a generar estadísticas de resumen y a crear gráficos 2D y hexadecimales.

Técnicas básicas de visualización de datos

• El paquete ggplot2 es una biblioteca que permite la creación de muchos tipos de visualización de datos, incluyendo varios tipos de gráficos y diagramas. Esta biblioteca fue creada creada por Hadley Wickham en 2005 y es una implementación en R de la Wilkinson’s Grammar of Graphics. La idea de este paquete es especificar bloques de construcción de gráficos y luego combinarlos para crear una pantalla gráfica.

• En este capítulo cubriremos los siguientes temas:

  • Creación de un gráfico de dispersión
  • Añadir una línea de regresión a un gráfico de dispersión
  • Trazado de categorías
  • Etiquetado del gráfico
  • Diseños de leyendas
  • Creación de una faceta
  • Tematización
  • Creación de gráficos de barras
  • Creación de gráficos de violines
  • Creación de gráficos de densidad

Creación de un gráfico de dispersión

• Un diagrama de dispersión es un gráfico en el que los valores de dos variables se trazan a lo largo de dos ejes, y el patrón de los puntos resultantes revela cualquier correlación presente

1. Los ejercicios de este capítulo requieren los siguientes paquetes: readr, dplyr, ggplot2. Pueden cargarse desde el panel de paquetes, el panel de la consola o un script

2. Abra RStudio y busque el panel de la consola

3. Cargar el contenido del archivo StudyArea.csv en un marco de datos

dfWildfires <- read_csv("C:/Data/StudyArea.csv", 
                        col_types = list(UNIT = col_character()), 
                        col_names = TRUE)

4. Crear un subconjunto de columnas

df <- select(dfWildfires, ORGANIZATI, STATE, YEAR_, TOTALACRES, CAUSE)

5. Agrupa los registros por año

grp <- group_by(df, YEAR_)

6. Resuma los datos según el número total de hectáreas quemadas

sm <- summarize(grp, totalacres = sum(TOTALACRES))

7. Utilice ggplot() para crear un gráfico de dispersión con el año en el eje x y el total de acres quemados en el eje Y.

ggplot(data=sm) + geom_point(mapping = aes(x=YEAR_, y=totalacres))

8. Hay ocasiones en las que tiene sentido utilizar las escalas logarítmicas en los cuadros y gráficos. Una de las razones es responder a la asimetría hacia los valores grandes, es decir, los casos en los que uno o unos pocos puntos son mucho más grandes que el grueso de los datos. En el gráfico que acabamos de crear hay un par de puntos que entran en esta categoría en en el eje Y. Vuelva a crear el gráfico, pero esta vez utilice la función log() en la columna totalacres

ggplot(data=sm) + geom_point(mapping = aes(x=YEAR_, y=log(totalacres)))

Añadir una línea de regresión al gráfico de dispersión

• Los gráficos construidos con ggplot() pueden tener más de una geometría. Es habitual añadir una línea de predicción (regresión) al gráfico.

1. Hay varias maneras de añadir una línea de regresión al gráfico de dispersión, una de ellas es utilizar la función geom_smooth() con el método establecido como lm (línea recta) y el parámetro se establecido en FALSE. Añade la línea de código que ves abajo en la ventana de la consola

ggplot(data=sm, aes(x=YEAR_, y=log(totalacres))) + 
  geom_point() +
  geom_smooth(method=lm, se=FALSE)

2. Cambia el método a loess el efecto en la línea de regresión

ggplot(data=sm, aes(x=YEAR_, y=log(totalacres))) + 
  geom_point() +
  geom_smooth(method=loess, se=FALSE)

3. Puede añadir un intervalo de confianza alrededor de la línea de regresión estableciendo se = TRUE

ggplot(data=sm, aes(x=YEAR_, y=log(totalacres))) + 
  geom_point() +
  geom_smooth(method=loess, se=TRUE)

Trazado de categorías

• En lugar de representar gráficamente todo el conjunto de incendios forestales, es posible que se quiera comprender mejor las tendencias por estado. En este paso creará un nuevo gráfico de dispersión que visualice las tendencias de los incendios forestales a lo largo del tiempo por estado

1. Reagrupar el marco de datos de los incendios forestales por estado y año

grp <- group_by(df, STATE, YEAR_)

2. Resuma los grupos por el total de hectáreas quemadas

sm <- summarize(grp, totalacres = sum(TOTALACRES))

3. Añade un parámetro de parámetro de color a la función aes() para que los puntos y la línea de regresión se se mapeen de acuerdo con el estado en el que se produjeron

ggplot(data=sm, aes(x=YEAR_, y=totalacres, colour=STATE)) + 
  geom_point(aes(colour = STATE)) + 
  stat_smooth(method=lm, se=FALSE)

Etiquetar el gráfico

• Puede añadir etiquetas a su gráfico mediante la función geom_text() o la función geom_ label()

1. Etiquete cada uno de los puntos del gráfico de dispersión utilizando geom_text() con un tamaño de etiqueta de 3

ggplot(data=sm, aes(x=YEAR_, y=log(totalacres))) + 
  geom_point() +
  geom_smooth(method=loess, se=TRUE) + 
  geom_text(aes(label=STATE), size=3)

La pantalla está extremadamente desordenada así que vamos a ajustar algunos parámetros para hacer esto más fácil de leer.

2. Puede utilizar el parámetro check_overlap para eliminar cualquier etiqueta superpuesta. Actualice su código como se ve a continuación

ggplot(data=sm, aes(x=YEAR_, y=log(totalacres))) + 
  geom_point() +
  geom_smooth(method=loess, se=TRUE) + 
  geom_text(aes(label=STATE), size=3, check_overlap = TRUE)

3. Esto se ve un poco mejor, pero si se cambia el tamaño de la etiqueta a 2 se reducirá aún más reducir el desorden y la superposición, mientras que esperamos que siga siendo legible

4. Habrás observado que las etiquetas se sitúan directamente sobre los temas. Puede utilizar los parámetros nudge_x y nudge_y para mover las etiquetas en relación con el punto. Utilice nudge_x como se ve a continuación para ver cómo esto mueve las etiquetas horizontalmente

ggplot(data=sm, aes(x=YEAR_, y=log(totalacres))) + 
  geom_point() +
  geom_smooth(method=loess, se=TRUE) + 
  geom_text(aes(label=STATE), size=2, check_overlap = TRUE, nudge_x = 1.0)

5. También puede colorear las etiquetas por categoría añadiendo el parámetro de color a la función aes() para geom_text()

ggplot(data=sm, aes(x=YEAR_, y=log(totalacres))) + 
  geom_point() +
  geom_smooth(method=loess, se=TRUE) + 
  geom_text(aes(label=STATE, color=STATE), size=2, check_overlap = TRUE, nudge_x = 1.0)

6. También puedes añadir un subtítulo y un pie de foto con el código que ves a continuación

ggplot(data=sm, aes(x=YEAR_, y=log(totalacres))) + 
  geom_point() +
  geom_smooth(method=loess, se=TRUE) + 
  labs(title=paste("Acreage Burned by Wildfires Has Increased In the Past Few Decades"), subtitle=paste("1980-2016"), caption="Data from USGS")

7. También puede actualizar las etiquetas X e Y del gráfico. Actualice estas etiquetas en su gráfico usando el código que ve abajo

ggplot(data=sm, aes(x=YEAR_, y=log(totalacres))) + 
  geom_point() +
  geom_smooth(method=loess, se=TRUE) + 
  labs(title=paste("Acreage Burned by Wildfires Has Increased In the Past Few Decades"), subtitle=paste("1980-2016"), caption="Data from USGS") +
  scale_y_continuous(name="Log of Total Acres Burned") + 
  scale_x_continuous(name="Burn Year")

Diseños de leyendas

• La función theme() puede utilizarse para controlar la ubicación de la leyenda y la función guide() puede utilizarse para proporcionar un control adicional de la leyenda

1. La función theme() junto con el argumento legend.postion se utiliza para controlar la ubicación de la leyenda en el gráfico. Por defecto, la leyenda que hemos visto hasta ahora se ha colocado en el lado derecho del gráfico con una orientación vertical. Reposicione la leyenda en la parte inferior con el código siguiente.

ggplot(data=sm, aes(x=YEAR_, y=log(totalacres), color=STATE)) +
  geom_point() + 
  labs(title=paste("Acreage Burned by Wildfires Has Increased In the Past Few Decades"),
  subtitle=paste("1980-2016"), caption="Data from USGS") + 
  scale_y_continuous(name="Log of Total Acres Burned") +
  scale_x_continuous(name="Burn Year") + 
  theme(legend.position="bottom")

2. También puede eliminar explícitamente una leyenda estableciendo legend.position = "none". Pruebe eso ahora si lo desea

3. Otros aspectos de la leyenda, como el número de filas en la leyenda, así como el tamaño del símbolo pueden ser controlados a través de la función guides(). Utilice el código que se ve a continuación para actualizar la leyenda a dos filas y con cada símbolo de tamaño 4

ggplot(data=sm, aes(x=YEAR_, y=log(totalacres), color=STATE)) +
  geom_point() +
  labs(title=paste("Acreage Burned by Wildfires Has Increased In the Past Few
Decades"), subtitle=paste("1980-2016"), caption="Data from USGS") +
  scale_y_continuous(name="Log of Total Acres Burned") +
  scale_x_continuous(name="Burn Year") +
  theme(legend.position = "bottom") +
  guides(color=guide_legend(nrow=2,override.aes=list(size=4)))

Crear una faceta

• Una forma particularmente buena de graficar variables categóricas es dividir el gráfico en facetas, que son subparcelas que muestran cada una un subconjunto de los datos. Las funciones facet_wrap() y facet_grid() pueden utilizarse para crear facetas

1. Utilice la función facet_wrap() que se muestra en el código siguiente para crear un mapa de facetas que muestre el total de acres quemados por estado

ggplot(data=sm, mapping = aes(x=YEAR_, y=log(totalacres))) + 
  geom_point() + facet_wrap(~STATE) + geom_smooth(method=loess, se=TRUE)

Temática

• Existen ocho temas incorporados que pueden utilizarse para personalizar el estilo de ggplot2 de los elementos que no son datos de su gráfico.

1. Los ocho temas incluidos en ggplot2 son theme_bw, theme_ classic, theme_dark, theme_gray, theme_light, theme_ linedraw, theme_minimal, theme_void. Añade el código que ves a continuación para cambiar la faceta a theme_dark

ggplot(data=sm, mapping = aes(x=YEAR_, y=log(totalacres))) + 
  geom_point() + facet_wrap(~STATE) + geom_smooth(method=loess, se=TRUE) +
  theme_dark()

2. Experimenta con los temas para ver las diferencias de estilo.

Creación de gráficos de barras

• Puede utilizar geom_bar() o geom_chart() para crear gráficos de barras con ggplot2. Sin embargo, hay una diferencia significativa entre ambos. La función geom_bar() generará un recuento del número de instancias de una variable. En otras palabras, cambia la estadística que ya se ha generado para el grupo. La función geom_col() mantiene la variable ya generada para el grupo. Para ver la diferencia, complete los siguientes pasos

1. Cargue el archivo StudyArea.csv y obtenga un subconjunto de columnas

dfWildfires <- read_csv("C:/Data/StudyArea.csv", 
                        col_types = list(UNIT = col_character()), col_names = TRUE)
df <- select(dfWildfires, ORGANIZATI, STATE, YEAR_, TOTALACRES, CAUSE)

2. Filtrar el marco de datos para que sólo se incluyan los incendios forestales de California

df <- filter(df, STATE == "California")

3. Agrupar el marco de datos por YEAR_

grp <- group_by(df, YEAR_)

4. Trace los datos utilizando geom_bar() como se ve a continuación. Observe que el gráfico de barras que se producido es un recuento del número de incendios para cada año

ggplot(data=grp) + geom_bar(mapping = aes(x=YEAR_), fill="red")

5. Ahora utilice geom_col() para ver la diferencia. La variable TOTALACRES se mantiene en este caso

ggplot(data=grp) + geom_col(mapping = aes(x=YEAR_, y=TOTALACRES), fill="red")

Creación de tramas de violín

• Los gráficos de violín, que son similares a los gráficos de caja, también muestran la densidad de probabilidad en varios valores. Las áreas más gruesas del gráfico de violín indican una mayor probabilidad en ese valor. Normalmente, los gráficos de violín también incluyen un marcador para la mediana junto con el rango intercuartil (IQR). La función geom_violin() se utiliza para crear gráficos violín en ggplot2

1. Cargue el archivo StudyArea.csv y obtenga un subconjunto de columnas

dfWildfires <- read_csv("C:/Data/StudyArea.csv", 
                        col_types = list(UNIT = col_character()), col_names = TRUE)
df <- select(dfWildfires, ORGANIZATI, STATE, YEAR_, TOTALACRES, CAUSE)

2. Filtrar el marco de datos para que sólo se incluyan los incendios forestales de más de 5.000 acres se incluyan.

dfWildfires <- filter(dfWildfires, TOTALACRES >= 5000)

3. Agrupe los incendios forestales por organización

grpWildfires <- group_by(dfWildfires, ORGANIZATI)

4. Crear una gráfica básica de violín.

ggplot(data=grpWildfires, mapping = aes(x=ORGANIZATI, y=log(TOTALACRES))) +
  geom_violin()

5. Puede añadir las observaciones individuales utilizando geom_jitter(). Si necesita mantener los puntos y dispersarlos horizontalmente, puede utilizar geom_jitter(height = 0)

ggplot(data=grpWildfires, mapping = aes(x=ORGANIZATI, y=log(TOTALACRES))) +
  geom_violin() + 
  geom_jitter(height = 0, width = 0.1)

6. La media puede añadirse utilizando stat_summary() como se ve a continuación

ggplot(data=grpWildfires, mapping = aes(x=ORGANIZATI, y=log(TOTALACRES))) +
  geom_violin() + geom_jitter(height = 0, width = 0.1) + 
  stat_summary(fun.y=median, geom="point", size=3, color="red") +
  stat_summary(fun.y=mean, geom="point", size=3, color="green")

7. La función box_plot() puede utilizarse para añadir la mediana y el IQR

ggplot(data=grpWildfires, mapping = aes(x=ORGANIZATI, y=log(TOTALACRES))) +
  geom_violin() + geom_boxplot(width=0.2)

Creación de gráficos de densidad

• Los gráficos de densidad, creados con geom_density() calculan una estimación de la densidad, que es una versión suavizada de un histograma y se utiliza con datos continuos. ggplot2 también puede calcular versiones 2D de la densidad, incluyendo contornos y gráficos de densidad con estilo de polígono

1. En esta primera parte del ejercicio crearás un gráfico de densidad básico. Cargue el archivo StudyArea.csv y obtenga un subconjunto de columnas

dfWildfires <- read_csv("C:/Data/StudyArea.csv", 
                        col_types = list(UNIT = col_character()), col_names = TRUE)
df <- select(dfWildfires, ORGANIZATI, STATE, YEAR_, TOTALACRES, CAUSE)

2. Filtrar el marco de datos para que sólo se incluyan los incendios forestales de más de 1.000 acres

dfWildfires <- filter(dfWildfires, TOTALACRES >= 1000)

3. Cree un gráfico de densidad con la función geom_density()

ggplot(dfWildfires, aes(TOTALACRES)) + geom_density()

4. También puede crear el mismo gráfico de densidad con una versión registrada de los datos

ggplot(dfWildfires, aes(log(TOTALACRES))) + geom_density()

5. A continuación, crearás gráficos 2D de los datos empezando por los contornos. Añade el código que ves a continuación

ggplot(dfWildfires, aes(x=YEAR_, y=log(TOTALACRES))) + 
  geom_point() +
  geom_density_2d()

6. Por último, cree una superficie de densidad 2D utilizando stat_density_2d()

ggplot(dfWildfires, aes(x=YEAR_, y=log(TOTALACRES))) +
  geom_density_2d() + 
  stat_density_2d(geom="raster", aes(fill=..density..), contour=FALSE)

Conclusiones

• En este capítulo has aprendido varias técnicas de visualización de datos utilizando ggplot2. Comenzamos con gráficos de dispersión básicos, añadimos líneas de regresión, etiquetamos los gráficos de varias formas y creamos una leyenda. Además, se aprendió a crear gráficos de facetas y a trabajar con ggplot2 y las opciones de tematización. También aprendió a crear gráficos de barras, gráficos de violín y gráficos de densidad.

Ejercicio para entregar

• Considere los datos asociados a precios de casas house_prices.csv que aparece en el directorio RDataSets.zip en Github. Nótese que el archivo posee su respectiva descripción house_prices_description.txt que puede ser de gran utilidad. Realice un análisis exploratorio de los datos relacionados con los precios de casas, teniendo en cuenta cada ítem estudiado en esta sección. Debería agregar todas las visualizaciones extras, que sean necesarias y le permitan generar y refinar preguntas sobre los datos.

• Considere los datos asociados a precios de casas house_prices.csv que aparece en el directorio RDataSets.zip en Github. A manera de aplicación, utilice las técnicas básicas de visualización de datos estudiadas en cada ítem de esta sección, ahora aplicadas al conjunto de datos relacionado con la predicción del precio de casas, el cual puede ser influenciado por los distintos factores, presentados en las columnas del archivo house_prices.csv, cuyas descripciones puede encontrar en el archivo house_prices_description.txt. Puede agregar visualizaciones extras, que le permitan generar y refinar preguntas sobre los datos.

Creación de un mapa base

• Hay dos pasos básicos para crear un mapa con ggmap. En general sólo hay que descargar el mapa rasterizado (basemap) y luego trazar los datos operativos en el mapa base. El primer paso se consigue con la función get_map(), que puede utilizarse para crear un mapa base desde Google, Stamen, Open Street Map o CloudMade. Aprenderás cómo hacerlo en este paso. En un próximo paso aprenderás a añadir y estilizar los datos operativos de varias maneras

1. Cargue el paquete ggmap yendo al panel de paquetes en RStudio y haciendo clic en la casilla junto al nombre del paquete. También puede cargarlo desde la consola escribiendo

library(ggmap)

2. Crea una variable llamada myLocation y defínela como California, ubicación que usaremos más adelante

myLocation <- "California"

3. Llama a la función get_map() y pasa la variable de localización junto con un nivel de Zoom de 6. Necesitaras primero crear una API key en Google Cloud Platform. Si descargas la versión de desarrollo de ggmap desde GitHub, hay una función que guarda tu clave de la API de Google y la utiliza en las siguientes llamadas a la API a través de ggmap. Para instalarla desde la consola escribir

devtools::install_github("dkahle/ggmap")

2. Utilice las siguientes instrucciones para descargar el mapa. Cargue primero su API key y luego use la función get_googlemap para obtener el mapa de interés. Para que las siguientes instrucciones funcionen, debe realizar los siguientes pasos en Google Cloud Platform

• Debe registrarse en Google Cloud Platform. Desafortunadamente para el registro deberá usar una Tarjeta de credito. Obtendrá $300 gratis para usarlos por 90 días. Luego de esto se cargarán a su cuenta $5 mensuales, dependiendo del uso que haga de la API. Recomiendo usar una tarjeta crédito prepago como por ejemplo Daviplata o similiares

• Debe habilitar los siguientes servicios en la Biblioteca

  • Maps Static API
  • Geocoding API
  • Geolocation API
  • Maps Embed API

• Para finalizar genere una API restringida a cada unos de los servicios de localización habilitados anteriormente. Las siguientes imágenes muestran como quedaría la configuración requerida

• Para obtener mas información sobre el uso de la función get_googlemap() ver ggmap. Por favor, no use la API key mykey = "AIzaSyABAXSsDR-f9qby3LK2o6bh0BypFWlKrm4", pues está será cambiada por otra, en este caso es usada a manera de ejemplo. Cree su propia API Google Cloud Platform

library(ggmap)

mykey = "AIzaSyABAXSsDR-f9qby3LK2o6bh0BypFWlKrm4"
register_google(key = mykey)

get_googlemap("waco, texas") %>% ggmap()

3. Lo siguientes tipos de mapas también están disponibles para vista satelital en los mapas. Para ver más opciones visitar la documentación de get_googlemap

get_googlemap("Colombia", 
              maptype = "hybrid",
              scale = 2,
              zoom = 5) %>% ggmap()

get_googlemap("Universidad del Norte", 
              maptype = "roadmap",
              scale = 2,
              zoom = 12) %>% ggmap()

get_googlemap("Universidad del Norte", 
              maptype = "terrain",
              scale = 2,
              zoom = 16) %>% ggmap()

get_googlemap("Universidad del Norte", 
              maptype = "hybrid",
              scale = 2,
              zoom = 18) %>% ggmap()

Añadir capas de datos operativos

• ggmap() devuelve un objeto ggplot, lo que significa que actúa como capa base en el ggplot2. Esto permite la gama completa de capacidades de ggplot2 lo que significa que puede trazar puntos en el mapa, añadir contornos, mapas de calor 2D y mucho más. Examinaremos algunas de estas capacidades en esta sección.

1. El siguiente código es usado para producir un mapa de California que muestre los incendios forestales entre los años 1980-2016 que quemaron más de 1,000 hectáreas.

myLocation <- "California"
myMap <- get_map(location = myLocation, zoom = 6)

library(readr)
dfWildfires <- read_csv("C:/Data/StudyArea_SmallFile.csv", col_names = TRUE)

library(dplyr)
df <- select(dfWildfires, STATE, YEAR_, TOTALACRES, DLATITUDE, DLONGITUDE)
knitr::kable(head(df))

df <- filter(df, TOTALACRES >= 1000 & STATE == "California")

ggmap(myMap) + geom_point(data=df, aes(x = DLONGITUDE, y = DLATITUDE))

2. Ahora vamos a hacer algo un poco más interesante. En primer lugar, utilizaremos la función de dplyr mutate() para agrupar los incendios por década

df <- mutate(df, 
             DECADE = ifelse(YEAR_ %in% 1980:1989, "1980-1989", 
                      ifelse(YEAR_ %in% 1990:1999, "1990-1999", 
                      ifelse(YEAR_ %in% 2000:2009, "2000-2009", 
                      ifelse(YEAR_ %in% 2010:2016, "2010-2016", "-99")))))

3. A continuación, codifica por colores los incendios forestales por DECADE y crea un mapa de símbolos graduados en función del tamaño de cada incendio. La propiedad de color define la columna a utilizar para la agrupación, y la propiedad de tamaño define la columna a utilizar para el tamaño de cada símbolo.

• Nótese que además es posible hacer dinámico el mapa usando ggplotly. Para esto utilizamos la función ggplot_build para convertir el plot de ggmap a un objeto ggplot.

last_plot <- ggmap(myMap) + 
  geom_point(data=df, aes(x = DLONGITUDE, y = DLATITUDE, colour= DECADE, size = TOTALACRES))

library(plotly)

gg <- ggplot_build(last_plot)$plot
plotly_plot <- ggplotly(gg)
plotly_plot

4. Vamos a cambiar la vista del mapa para centrarnos más en el sur de California, y en en particular el área al norte de Los Ángeles.

myMap <- get_map(location = "Santa Clarita, California", zoom = 10)

ggmap(myMap) + 
  geom_point(data=df, 
             aes(x = DLONGITUDE, y = DLATITUDE, colour= DECADE, size = TOTALACRES))

5. A continuación, añadiremos capas de contorno y de calor. La función geom_density2d() se utiliza para crear los contornos mientras que la función stat_density2d() crea el mapa de calor. Puedes experimentar con los colores usando las propiedades scale_fill_gradient(low and high).

• Aquí los hemos establecido en verde y rojo respectivamente, pero puede cambiar el esquema de colores según su preferencia. En este ejemplo, la sintaxis fill = ..level.. se utiliza para rellenar los contornos en el gráfico basado en los niveles de densidad de los puntos de datos.

myMap <- get_map(location = "California", zoom = 6)
ggmap(myMap, extent = "device") + 
  geom_density2d(data = df, aes(x =DLONGITUDE, y = DLATITUDE), size = 0.3) +
  stat_density2d(data = df, 
                 aes(x = DLONGITUDE, y = DLATITUDE, fill = ..level.., 
                     alpha =..level..), 
                 size = 0.01, bins = 16, geom = "polygon") + 
  scale_fill_gradient(low = "green", high ="red") + 
  scale_alpha(range = c(0, 0.3), guide = FALSE)

6. Si prefiere ver el mapa de calor sin contornos, el código puede simplificarse como sigue

ggmap(myMap, extent = "device") + 
  stat_density2d(data = df, 
                 aes(x =DLONGITUDE, y = DLATITUDE, fill = ..level.., 
                     alpha = ..level..), 
                 size = 0.01,bins = 16, geom = "polygon") + 
  scale_fill_gradient(low = "green", high ="red") + 
  scale_alpha(range = c(0, 0.3), guide = FALSE)

7. Por último, vamos a crear un mapa de facetas que represente los puntos críticos de cada año de la década actual. El conjunto de datos contiene información hasta el año 2016. facet_wrap() es una función útil en ggplot2 para crear múltiples mapas mostrados juntos en un diseño de cuadrícula, cada uno representando un subconjunto de los datos.

df <- filter(df, YEAR_ %in% c(2010, 2011, 2012, 2013, 2014, 2015, 2016))

gg <- ggmap(myMap, extent = "device") + 
  stat_density2d(data = df, aes(x = DLONGITUDE, y = DLATITUDE, fill = ..level.., alpha = ..level..),
                 size = 0.01, bins = 16, geom = "polygon") +
  scale_fill_gradient(low = "green", high ="red") + 
  scale_alpha(range = c(0, 0.3), guide = FALSE) + 
  facet_wrap(~YEAR_, ncol = 4)
  theme(plot.title = element_text(size = 20)) +
  theme(plot.margin = margin(1, 1, 1, 1, "cm"))

## List of 2
##  $ plot.title :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : num 20
##   ..$ hjust        : NULL
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : NULL
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi FALSE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ plot.margin: 'margin' num [1:4] 1cm 1cm 1cm 1cm
##   ..- attr(*, "unit")= int 1
##  - attr(*, "class")= chr [1:2] "theme" "gg"
##  - attr(*, "complete")= logi FALSE
##  - attr(*, "validate")= logi TRUE

print(gg)

Añadir capas a partir de Shapefiles

• Aunque se trata de un formato de datos SIG algo más antiguo, los shapefiles siguen siendo utilizados comúnmente para representar características geográficas. Con un poco de manipulación, puede obtener datos de trazado de shapefiles en ggmap.

1. Para este ejercicio necesitarás instalar un paquete adicional llamado rgdal. Utilice el panel de paquetes para encontrar e instalar rgdal o introduzca el código que ve a continuación

install.packages("rgdal")

2. Cargue el paquete rgdal a través del panel de paquetes o introduzca el código que ve a continuación

library(rgdal)

3. La carpeta RDataSets que contiene los datos del ejercicio para este curso, contiene un archivo shape llamado S_USA.Wilderness.

• Estos archivos se combinan para crear lo que se llama un shapefile. Este archivo contiene los límites de las áreas silvestres designadas en los Estados Unidos. Utilice la función readOGR() de rgdal para cargar los datos en una variable.

wild <- readOGR("C:/Data", "S_USA.Wilderness")

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Data", layer: "S_USA.Wilderness"
## with 447 features
## It has 8 fields
## Integer64 fields read as strings:  WID

4. La función fortify(), que forma parte de ggplot2, convierte todos los puntos individuales que definen cada límite en un marco de datos que se puede utilizar para trazar los límites del polígono.

wild <- fortify(wild)

5. Utilice la función ggmap qmap() (qmap significa mapa rápido) para crear el mapa base que se utilizará como referencia para los límites de los espacios naturales. Centra el mapa en Montana.

montana <- qmap("Montana", zoom=6)

6. Antes de trazar los límites de los espacios naturales como polígonos en el mapa, observe el marco de datos creado por la función fortify() para entender mejor la estructura creada por esta función

knitr::kable(head(wild)) 

• Fíjese en la columna group. Esta columna identifica de forma exclusiva límites de los espacios naturales. Los límites de los espacios naturales son polígonos, y los polígonos están definidos por un conjunto de puntos que definen la estructura del polígono. Es algo así como como jugar a conectar los puntos, donde cada punto es un par de coordenadas de latitud/longitud definidas por las columnas long y lat en el marco de datos.

• Por ejemplo, mire el grupo 0.1. Observe que hay varias filas que contiene el valor 0.1, y que cada fila tiene valores únicos de longitud y lat. Estos son todos los puntos utilizados para definir los límites de ese polígono.

7. Ahora traza los límites de los espacios naturales en el mapa base. Observe el uso de la columna para agrupar los polígonos. El trazado de los límites en el mapa lleva algo de tiempo, así que sea paciente con este paso. Al final debería ver un mapa similar a la captura de pantalla de abajo.

montana + 
  geom_polygon(aes(x=long,y=lat, group=group, alpha=0.25), data=wild, fill="white") + 
  geom_polygon(aes(x=long, y=lat, group=group), data=wild, color="black", fill=NA)

Conclusión

• En este capítulo ha aprendido a utilizar el paquete ggmap para crear atractivas visualizaciones de datos en formato de mapa.

• Ha aprendido a *crear mapas base utilizando Google como fuente de datos, añadir capas de datos operativos, crear varios tipos de visualizaciones de mapas utilizando fuentes de datos externas, y cargar archivos de forma.

Ejercicio para entregar

• Considere el conjunto de datos violent_crimes.csv en RDataSets.zip el cual abarca una gran variedad de delitos, y realice visualizaciones geográficas para cada tipo de crimen.

• Considere cada ítem estudiado en esta sección para realizar las visualizaciones respectivas con este dataset.

• Verifique si existe algún patrón dentro de los delitos violentos que pueda explorarse visualmente. Los datos incluyen dentro de sus columnas la longitud y latitud, la categoría en la que se clasifica el delito, la fecha y hora.

Introducción

• Los datos faltantes es uno de los temas que se ignoran en la mayoría de los textos introductorios. Probablemente, parte de la razón por la que esto es así es que todavía abundan muchos mitos sobre el análisis con datos faltantes. Además, algunas de las investigaciones sobre técnicas de vanguardia es aún relativamente nueva. Una razón más legítima para su ausencia en los textos introductorios es que la mayoría de las metodologías son bastante complicadas, desde el punto de vista matemático. Sin embargo, la increíble ubicuidad de los problemas relacionados con los datos faltantes en el análisis de datos de la vida real requiere que se aborde el tema.

• Esta sección sirve como una suave introducción al tema y a una de las técnicas más eficaces para tratarla. Un refrán común sobre el tema es algo así como que, la mejor manera de tratar con los datos que faltan es no tenerlos. Es cierto que los datos que faltan son un tema complicado, y hay muchas maneras de hacerlo, de las cuales algunas pueden estar mal. Es importante no llevar este consejo al extremo, sin embargo, para eludir los problemas de los datos faltantes, algunos han impedido participantes en una encuesta, que por ejemplo, siguen sin responder a todas las preguntas de un formulario.

• Hay tratamientos para los datos que faltan, pero no hay tratamientos para los datos malos o incorrectos. El tratamiento estándar para el problema de los datos que faltan es reemplazar los datos que faltan por valores no ausentes. Este proceso se denomina imputación. En la mayoría de los casos, el objetivo de la imputación no es recrear el conjunto de datos faltantes completo, sino permitir que se realicen estimaciones o inferencias estadísticas válidas a partir de los datos faltantes.

• Por ello, la eficacia de las diferentes técnicas de imputación no puede evaluarse por su capacidad de recrear los datos con la mayor exactitud posible a partir de un conjunto de datos faltantes simulado, sino que deben juzgarse por su capacidad de apoyar las mismas inferencias estadísticas que se obtendrían del análisis de los datos completos que se extraen del análisis.

• De este modo, rellenar los datos que faltan es sólo un paso hacia el verdadero objetivo: el análisis. El conjunto de datos imputados rara vez se considera el objetivo final de la imputación. En la práctica, hay muchas formas diferentes de tratar los datos que faltan, algunas son buenas y otras no tanto. Algunas están bien en determinadas circunstancias, pero no en otras. Algunas implican la eliminación de datos faltantes, mientras que otras implican la imputación. Vamos a mencionar brevemente algunos de los métodos más comunes. Sin embargo, el objetivo final de esta sección, es iniciarle en lo que a menudo se describe como el estándar de oro de las técnicas de imputación: la imputación múltiple.

Visualización de los datos que faltan

• Para demostrar la visualización de patrones de datos ausentes, primero tenemos que crear algunos datos ausentes. Para mostrar cómo utilizar la imputación múltiple en un escenario semirealista, vamos a crear una versión del conjunto de datos mtcars con algunos valores faltantes. Configuramos seed (para la aleatoriedad determinista), y creamos una variable para mantener nuestro nuevo conjunto de datos.

• Los datos se extrajeron de la revista Motor Trend US de 1974 y comprenden el consumo de combustible y 10 aspectos del diseño y el rendimiento de 32 automóviles (modelos de 1973 a 1974) (ver MTCARS).

set.seed(2)
miss_mtcars <- mtcars

• En primer lugar, vamos a crear siete valores faltantes en drat (alrededor del 20 por ciento), cinco valores faltantes en la columna mpg (alrededor del 15 por ciento), cinco valores faltantes en la columna cyl, tres valores faltantes en wt (alrededor del 10 por ciento), y tres valores faltantes en vs

some_rows <- sample(1:nrow(miss_mtcars), 7)
miss_mtcars$drat[some_rows] <- NA

some_rows <- sample(1:nrow(miss_mtcars), 5)
miss_mtcars$mpg[some_rows] <- NA

some_rows <- sample(1:nrow(miss_mtcars), 5)
miss_mtcars$cyl[some_rows] <- NA

some_rows <- sample(1:nrow(miss_mtcars), 3)
miss_mtcars$wt[some_rows] <- NA

some_rows <- sample(1:nrow(miss_mtcars), 3)
miss_mtcars$vs[some_rows] <- NA

• Ahora, vamos a crear cuatro valores faltantes en qsec, pero sólo para los coches automáticos

only_automatic <- which(miss_mtcars$am==0)
some_rows <- sample(only_automatic, 4)
miss_mtcars$qsec[some_rows] <- NA

• Ahora, observemos las primeras filas del nuevo conjunto de datos:

head(miss_mtcars)

##                    mpg cyl disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4         21.0  NA  160 110 3.90 2.620 16.46 NA  1    4    4
## Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
## Datsun 710        22.8  NA  108  93 3.85 2.320 18.61  1  1    4    1
## Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
## Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
## Valiant           18.1   6  225 105   NA 3.460 20.22  1  0    3    1

str(miss_mtcars)

## 'data.frame':    32 obs. of  11 variables:
##  $ mpg : num  21 21 22.8 21.4 18.7 18.1 14.3 24.4 NA 19.2 ...
##  $ cyl : num  NA 6 NA 6 8 6 8 4 4 6 ...
##  $ disp: num  160 160 108 258 360 ...
##  $ hp  : num  110 110 93 110 175 105 245 62 95 123 ...
##  $ drat: num  3.9 3.9 3.85 3.08 3.15 NA 3.21 NA 3.92 3.92 ...
##  $ wt  : num  2.62 2.88 2.32 3.21 3.44 ...
##  $ qsec: num  16.5 17 18.6 19.4 17 ...
##  $ vs  : num  NA 0 1 1 0 1 NA 1 1 1 ...
##  $ am  : num  1 1 1 0 0 0 0 0 0 0 ...
##  $ gear: num  4 4 4 3 3 3 3 4 4 4 ...
##  $ carb: num  4 4 1 1 2 1 4 2 2 4 ...

• Ahora vamos a visualizar los datos faltantes. La primera forma en que vamos a visualizar el patrón de los datos faltantes es utilizando la función md.pattern del paquete mice (que es también el paquete que vamos a utilizar para imputar nuestros datos faltantes). Si no tiene el paquete instálelo antes

library(mice)
md.pattern(miss_mtcars, plot = TRUE, rotate.names = TRUE)

##    disp hp am gear carb wt vs qsec mpg cyl drat   
## 12    1  1  1    1    1  1  1    1   1   1    1  0
## 6     1  1  1    1    1  1  1    1   1   1    0  1
## 2     1  1  1    1    1  1  1    1   1   0    1  1
## 3     1  1  1    1    1  1  1    1   0   1    1  1
## 1     1  1  1    1    1  1  1    1   0   0    1  2
## 1     1  1  1    1    1  1  1    0   1   1    1  1
## 1     1  1  1    1    1  1  1    0   0   1    1  2
## 1     1  1  1    1    1  1  0    1   1   0    1  2
## 1     1  1  1    1    1  1  0    0   1   1    1  2
## 1     1  1  1    1    1  1  0    0   1   0    1  3
## 2     1  1  1    1    1  0  1    1   1   1    1  1
## 1     1  1  1    1    1  0  1    1   1   1    0  2
##       0  0  0    0    0  3  3    4   5   5    7 27

• Un patrón de datos faltantes por fila se refiere a las columnas que faltan en cada fila. Esta función agrega y contabiliza el número de filas con el mismo patrón de datos faltantes. Esta función produce una matriz binaria (0 y 1). Las celdas con un 1 representan datos no faltantes; los 0 representan datos que faltan. Como las filas están ordenadas en un orden creciente de ausencia, la primera fila siempre se refiere al patrón de datos ausentes que contiene la menor cantidad de datos faltantes.

• La columna de la izquierda es un conteo del número de filas que muestran el patrón de datos faltantes, y la columna más a la derecha es un recuento del número de puntos faltantes en ese patrón. La última fila contiene un recuento del número de puntos faltantes en cada columna. Como puede ver, 12 de las filas no contienen datos faltantes.

• Sólo hay seis filas que contienen más de un valor perdido. Sólo una de estas filas contiene más de dos valores faltantes (como se muestra en la penúltima fila). En cuanto a los conjuntos de datos con datos faltantes, este en particular no contiene mucho. No es raro que en algunos conjuntos de datos falte más del 30% de los datos. Este conjunto de datos no llega ni al 3%. Ahora vamos a visualizar el patrón de datos faltantes gráficamente utilizando el paquete VIM, el cual debe instalar en caso de no contar con este.

library(VIM)
aggr(miss_mtcars, numbers=TRUE)

• A simple vista, esta representación nos muestra, sin esfuerzo, que la columna drat representa la mayor proporción de datos faltantes, por columnas, seguida de mpg, cyl, qsec, vs y wt. El gráfico de la derecha nos muestra información similar a la de la salida de md.pattern. Esta representación, sin embargo, hace que sea más fácil saber si hay algún patrón sistemático de omisión. Las celdas azules representan los datos no ausentes, y las rojas representan los datos que faltan. Los números de la derecha del gráfico representan la proporción de filas que muestran ese patrón de datos faltantes. El 37.5% de las filas no contienen ningún tipo de dato faltante.

marginplot(miss_mtcars[c(1,2)])

Tipos de datos faltantes

• El paquete VIM nos permitió visualizar los patrones de datos faltantes. Un término relacionado, al mecanismo de datos faltantes, describe el proceso que determina la probabilidad de que cada punto de datos sea faltante. Hay tres categorías principales de mecanismos de datos faltantes: Missing Completely At Random (MCAR), Missing At Random (MAR), y Missing Not At Random (MNAR). La discriminación basada en el mecanismo de datos faltantes es crucial, ya que nos informa sobre las opciones para manejar los datos faltantes.

  • Faltantes completamente al azar (MCAR)

    • El primer mecanismo, MCAR, se produce cuando la falta de datos no está relacionada con la datos. Esto ocurriría, por ejemplo, si se borraran filas de una base de datos al azar, o si una ráfaga de viento se llevara una muestra aleatoria de los formularios de encuesta de un investigador. El mecanismo que rige la ausencia de drat, mpg, cyl, wt, y vs es MCAR, porque seleccionamos al azar los elementos que faltan. Este mecanismo, aunque es el más fácil de trabajar, rara vez es sostenible en la práctica.

  • Faltantes no al azar (MNAR)

    • El mecanismo MNAR se produce cuando la ausencia de una variable está relacionada con la variable en sí misma. Por ejemplo, supongamos que la báscula que pesa cada coche tiene una capacidad de sólo 3.700 libras y, por ello, los ocho coches que pesaban más se registraron como NA. Este es un ejemplo clásico del mecanismo MNAR, es el propio peso de la observación la causa de que falte.
    • Otro ejemplo sería si durante el transcurso del ensayo de un fármaco antidepresivo los participantes a los que no les ayudaba el fármaco se deprimían demasiado para continuar con el ensayo. Al final del ensayo, cuando se accede y se registra el nivel de depresión de todos los participantes, habrá valores que falten para los participantes cuyo motivo de ausencia está relacionado con su nivel de depresión.

  • Faltantes al azar (MAR)

    • El mecanismo, faltante al azar, tiene un nombre un tanto desafortunado, al contrario de lo que pueda parecer, significa que existe una relación sistemática entre la ausencia de una variable resultado y otras variables observadas, pero no la variable de resultado en sí misma.

    • Un ejemplo es que los hombres son menos propensos a completar una encuesta sobre depresión, pero esto no tiene nada que ver con su nivel de depresión, después de tener en cuenta la masculinidad.

    • Por ejemplo, puede que todas las observaciones de nuestro conjunto de datos con un valor ausente de fiebre no se hayan registrado porque se asumió que todos los pacientes con escalofríos y dolores tenían fiebre, por lo que nunca se les tomó la temperatura. Si es cierto, podríamos predecir fácilmente que cada observación que falta con escalofríos y dolores también tiene fiebre y utilizar esta información para imputar nuestros datos que faltan.

• Observación

  • Quizá haya observado que el lugar que ocupa un conjunto de datos concreto en la taxonomía del mecanismo de datos faltantes depende de las variables que incluye. Por ejemplo, sabemos que el mecanismo detrás de qsec es MAR, pero si el conjunto de datos no incluyera am, sería MNAR. Como somos nosotros los que creamos los datos, sabemos el procedimiento que dio lugar a los valores faltantes de qsec. Si no fuéramos nosotros los que creamos los datos, como ocurre en el mundo real, y el conjunto de datos no contuviera la columna am simplemente veríamos una cantidad de valores qsec que faltan arbitrariamente.
  • Esto podría llevarnos a creer que los datos son MCAR. Sin embargo, no lo es; sólo porque la variable con la cual otra variable faltante está sistemáticamente relacionada no se observa, no quiere decir que esta no exista. Esto plantea una cuestión crítica: ¿Podemos estar seguros de que nuestros datos no son MCAR? La respuesta desafortunada es no. Dado que los datos que necesitamos para demostrar o refutar el MNAR faltan ipso facto, la suposición de MNAR nunca puede ser desconfirmada de forma concluyente. Es nuestro trabajo, como analistas de datos con pensamiento crítico, preguntar si es probable que haya un mecanismo MNAR o no.

data <- airquality
knitr::kable(head(data))

# data <- airquality
data[4:10, 3] <- rep(NA, 7)
data[1:5, 4] <- NA

data <- data[-c(5,6)]
summary(data)

##      Ozone           Solar.R           Wind             Temp      
##  Min.   :  1.00   Min.   :  7.0   Min.   : 1.700   Min.   :57.00  
##  1st Qu.: 18.00   1st Qu.:115.8   1st Qu.: 7.400   1st Qu.:73.00  
##  Median : 31.50   Median :205.0   Median : 9.700   Median :79.00  
##  Mean   : 42.13   Mean   :185.9   Mean   : 9.806   Mean   :78.28  
##  3rd Qu.: 63.25   3rd Qu.:258.8   3rd Qu.:11.500   3rd Qu.:85.00  
##  Max.   :168.00   Max.   :334.0   Max.   :20.700   Max.   :97.00  
##  NA's   :37       NA's   :7       NA's   :7        NA's   :5

knitr::kable(head(data))

pMiss <- function(x){sum(is.na(x))/length(x)*100}
apply(data,2,pMiss)

##     Ozone   Solar.R      Wind      Temp 
## 24.183007  4.575163  4.575163  3.267974

library(mice)
md.pattern(data, rotate.names = TRUE)

##     Temp Solar.R Wind Ozone   
## 104    1       1    1     1  0
## 34     1       1    1     0  1
## 3      1       1    0     1  1
## 1      1       1    0     0  2
## 4      1       0    1     1  1
## 1      1       0    1     0  2
## 1      1       0    0     1  2
## 3      0       1    1     1  1
## 1      0       1    0     1  2
## 1      0       0    0     0  4
##        5       7    7    37 56

library(VIM)

aggr_plot <- aggr(data, col=c('navyblue','red'), numbers=TRUE, sortVars=TRUE, labels=names(data),
                  cex.axis=.7, gap=3, ylab=c("Histogram of missing data","Pattern"))

## 
##  Variables sorted by number of missings: 
##  Variable      Count
##     Ozone 0.24183007
##   Solar.R 0.04575163
##      Wind 0.04575163
##      Temp 0.03267974

marginplot(data[c(1,2)])

airquality_omit <- na.omit(data)
head(airquality_omit)

##    Ozone Solar.R Wind Temp
## 12    16     256  9.7   69
## 13    11     290  9.2   66
## 14    14     274 10.9   68
## 15    18      65 13.2   58
## 16    14     334 11.5   64
## 17    34     307 12.0   66

library(ggplot2)
library(tidyverse)
library(hrbrthemes)
library(gridExtra)

ggp1 <- ggplot(data.frame(value=data$Ozone), aes(x=value)) +
  geom_histogram(fill="#FBD000", color="#E52521", alpha=0.9) +
  ggtitle("Original data") +
  xlab('Ozone') + ylab('Frequency') +
  theme_ipsum() +
  theme(plot.title = element_text(size=15))

ggp2 <- ggplot(data.frame(value=airquality_omit$Ozone), aes(x=value)) +
  geom_histogram(fill="#43B047", color="#049CD8", alpha=0.9) +
  ggtitle("Listwise Deletion") +
  xlab('Ozone') + ylab('Frequency') +
  theme_ipsum() +
  theme(plot.title = element_text(size=15))

grid.arrange(ggp1, ggp2, ncol = 2)

Imputación por emparejamiento predictivo medio

• Algunos científicos de datos o estadísticos pueden buscar una solución rápida sustituyendo los datos que faltan por la media. La media se utiliza a menudo para imputar datos no categóricos. Por ejemplo, en el conjunto de datos de calidad del aire, supongamos que queremos imputar la media de sus valores faltantes. En este caso, utilizamos el paquete R mice. Cambiando el argumento method = mean, se especifica la imputación de la media, el argumento m = 1 cambia las iteraciones a 1, lo que significa no (iteración).

• El fundamento teórico de utilizar la media para imputar los datos faltantes es que la media es una buena estimación para seleccionar aleatoriamente una observación de una distribución normal. Ahora, intentaremos imputar la media para las variables Ozono y Solar.R del conjunto de datos airquality. En primer lugar, vamos a cargar los paquetes mice y mipackages (instale el paquete de CRAN primero usando install.packages("mice") y install.packages("mi")). Para mas información sobre el paquete utilizado (ver Mice). Podemos recuperar el conjunto de datos completo utilizando la función complete().

library(mice)
library(foreign)

• Los datos faltantes para Ozono y Solar.R pueden ser imputados por la función mice().

imp <- mice(data, m=5, maxit=50, method ='pmm', seed=500, printFlag = FALSE)

• El parámetro m se refiere al número de conjuntos de datos imputados que se generarán. Es decir, m representa la cantidad de veces que se imputarán los datos para crear múltiples conjuntos de datos imputados. Cada conjunto de datos imputados contendrá valores imputados diferentes para los datos faltantes, lo que permite capturar la incertidumbre en el proceso de imputación o falta de certeza sobre los valores reales de los datos faltantes.

• El parámetro maxit en el paquete mice de R controla cuántas veces se repetirá el proceso de imputación en cada conjunto de datos para alcanzar una convergencia o un límite máximo de iteraciones. Durante cada iteración, el algoritmo intenta mejorar la imputación de los valores faltantes. Si el número de iteraciones alcanza el valor especificado en maxit y el proceso no converge, la imputación se detendrá y se generará un mensaje de advertencia. Ajustar adecuadamente el valor de maxit puede ayudar a controlar el tiempo de ejecución y garantizar una convergencia adecuada del algoritmo de imputación.

• El valor por defecto es m=5. method='pmm' (ver mice.impute.pmm) se refiere al método de imputación. En este caso, utilizamos el método de imputación por emparejamiento predictivo medio (PMM). Se pueden utilizar otros métodos de imputación, escriba methods(mice) para obtener una lista de los métodos de imputación disponibles.

methods(mice)

##  [1] mice.impute.2l.bin              mice.impute.2l.lmer            
##  [3] mice.impute.2l.norm             mice.impute.2l.pan             
##  [5] mice.impute.2lonly.mean         mice.impute.2lonly.norm        
##  [7] mice.impute.2lonly.pmm          mice.impute.cart               
##  [9] mice.impute.jomoImpute          mice.impute.lasso.logreg       
## [11] mice.impute.lasso.norm          mice.impute.lasso.select.logreg
## [13] mice.impute.lasso.select.norm   mice.impute.lda                
## [15] mice.impute.logreg              mice.impute.logreg.boot        
## [17] mice.impute.mean                mice.impute.midastouch         
## [19] mice.impute.mnar.logreg         mice.impute.mnar.norm          
## [21] mice.impute.mpmm                mice.impute.norm               
## [23] mice.impute.norm.boot           mice.impute.norm.nob           
## [25] mice.impute.norm.predict        mice.impute.panImpute          
## [27] mice.impute.passive             mice.impute.pmm                
## [29] mice.impute.polr                mice.impute.polyreg            
## [31] mice.impute.quadratic           mice.impute.rf                 
## [33] mice.impute.ri                  mice.impute.sample             
## [35] mice.mids                       mice.theme                     
## see '?methods' for accessing help and source code

imp_df <- complete(imp)
head(imp_df)

##   Ozone Solar.R Wind Temp
## 1    41     190  7.4   87
## 2    36     118  8.0   80
## 3    12     149 12.6   74
## 4    18     313 10.9   66
## 5    13      81 16.6   57
## 6    28      78  7.4   66

• La siguiente imagen ilustra el funcionamiento del método de imputación múltiple por ecuaciones encadenadas

• Además, la siguiente figura ilustra el funcionamiento del método de imputación por emparejamiento predictivo medio usado en el ejemplo para airquality

• Observación: El método PMM utiliza aleatoriedad en la imputación para introducir variabilidad en las estimaciones de los valores faltantes. Esto ayuda a evitar la sobreestimación o subestimación sistemática de los valores imputados y a capturar la incertidumbre asociada con la imputación de datos faltantes. La introducción de aleatoriedad también puede ayudar a mejorar la robustez y generalización del modelo de imputación, especialmente en situaciones donde los datos faltantes siguen patrones complejos o no aleatorios.

• Tracemos histograma y diagrama de dispersión del aspecto del conjunto de datos de la calidad del aire tras la imputación mediante método PMM. Los siguientes son los histogramas asociados

ggp1 <- ggplot(data.frame(value=data$Ozone), aes(x=value)) +
  geom_histogram(fill="#FBD000", color="#E52521", alpha=0.9) +
  ggtitle("Original data") +
  xlab('Ozone') + ylab('Frequency') +
  theme_ipsum() +
  theme(plot.title = element_text(size=15))

ggp2 <- ggplot(data.frame(value=imp_df$Ozone), aes(x=value)) +
  geom_histogram(fill="#43B047", color="#049CD8", alpha=0.9) +
  ggtitle("PMM imputation") +
  xlab('Ozone') + ylab('Frequency') +
  theme_ipsum() +
  theme(plot.title = element_text(size=15))

grid.arrange(ggp1, ggp2, ncol = 2)

• Vamos a comparar las distribuciones de los datos originales e imputados utilizando algunos gráficos útiles. En primer lugar, podemos utilizar un gráfico de dispersión y comparar Ozone con todas las demás variables.

library(lattice)
xyplot(imp, Ozone ~ Wind + Temp + Solar.R, pch=18, cex=1)

• Los valores de pch están asociados con los símbolos gráficos desplegados. cex controla el tamaño de los símbolos, por defecto el tamaño es 1, podemos poner valores menores o superiores.

• Lo que queremos ver es que la forma de los puntos rojos (imputados) coincida con la de los azules (observados). La coincidencia de la forma nos indica que los valores imputados son realmente valores plausibles. Otro gráfico útil es el de la densidad:

densityplot(imp)

• La densidad de los datos imputados para cada conjunto de datos imputados se muestra en color rojo, mientras que la densidad de los datos observados se muestra en azul. De nuevo, según nuestros supuestos anteriores, esperamos que las distribuciones sean similares.

• Es muy recomendable comprobar visualmente la convergencia. Lo comprobamos cuando llamamos a la función de trazado a la que asignamos la salida de mice, para mostrar gráficos de seguimiento de la media y desviación estándar de todas las variables implicadas en las imputaciones.

plot(imp);

• Cada línea es una de las m conjuntos de datos imputados, a saber m=5 por defecto. Como se puede ver en el gráfico de trazado anterior sobre imp, no hay tendencias claras y las variables se superponen de una iteración a la siguiente. Dicho de otro modo, la varianza dentro de una cadena (hay m cadenas) debería ser aproximadamente igual a la varianza entre las cadenas. Esto indica que se ha logrado la convergencia. Si no se logra la convergencia, se puede aumentar el número de iteraciones que mice emplea especificando explícitamente el parámetro maxit a la función mice.

data <- data[, c("Solar.R", "Wind")]
imp.regress <- mice(data, method="norm.predict", m=1, maxit=1)

## 
##  iter imp variable
##   1   1  Solar.R  Wind

imp.regress$imp$Wind

##           1
## 4  9.736111
## 5  9.806622
## 6  9.805925
## 7  9.743832
## 8  9.854140
## 9  9.898263
## 10 9.801744

Detección de valores atípicos

• Un valor atípico es un valor o una observación que se aleja de otras observaciones, es decir, un punto de datos que difiere significativamente de otros puntos de datos. Enderlein (1987) va incluso más allá, ya que el autor considera que los valores atípicos son aquellos que se desvían tanto de otras observaciones que se podría suponer un mecanismo de muestreo subyacente diferente.

• Una observación debe compararse siempre con otras observaciones realizadas sobre el mismo fenómeno antes de calificarla realmente de atípica. En efecto, una persona de 200 cm de altura (1,90 m en EE.UU.) se considerará probablemente un valor atípico en comparación con la población general, pero esa misma persona podría no considerarse un valor atípico si midiéramos la altura de los jugadores de baloncesto.

• En esta sección, presentamos varios enfoques para detectar valores atípicos en R, desde técnicas sencillas como la estadística descriptiva (que incluye el mínimo, el máximo, el histograma, el boxplot y los percentiles) hasta técnicas más formales como el filtro de Hampel, el Grubbs, el Dixon y las pruebas de Rosner para detectar valores atípicos. Algunas pruebas estadísticas exigen la ausencia de valores atípicos para sacar conclusiones sólidas, pero la eliminación de valores atípicos no se recomienda en todos los casos y debe hacerse con precaución.

dat <- ggplot2::mpg
head(dat)

## # A tibble: 6 × 11
##   manufacturer model displ  year   cyl trans      drv     cty   hwy fl    class 
##   <chr>        <chr> <dbl> <int> <int> <chr>      <chr> <int> <int> <chr> <chr> 
## 1 audi         a4      1.8  1999     4 auto(l5)   f        18    29 p     compa…
## 2 audi         a4      1.8  1999     4 manual(m5) f        21    29 p     compa…
## 3 audi         a4      2    2008     4 manual(m6) f        20    31 p     compa…
## 4 audi         a4      2    2008     4 auto(av)   f        21    30 p     compa…
## 5 audi         a4      2.8  1999     6 auto(l5)   f        16    26 p     compa…
## 6 audi         a4      2.8  1999     6 manual(m5) f        18    26 p     compa…

summary(dat$hwy)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   12.00   18.00   24.00   23.44   27.00   44.00

min(dat$hwy)

## [1] 12

max(dat$hwy)

## [1] 44