Introducción a R
Ángel Borrego
10/01/2022
Este documento presenta los materiales del curso de introducción a R para miembros del Centre de Recerca en Informació, Comunicació i Cultural (CRICC) y estudiantes de máster y doctorado de la Facultat d’Informació i Mitjans Audiovisuals de la Universitat de Barcelona.
1 Introducción
1.1 ¿Qué es R i RStudio?
R es un entorno y lenguaje de programación con un enfoque al análisis estadístico muy popular para el análisis de datos masivos. Cualquier usuario puede publicar paquetes que extienden la configuración básica.
El objetivo de este curso es familiarizarnos con el entorno de trabajo de R y el integrated development environment (IDE) RStudio. Podemos descargar ambos programas desde las siguientes direcciones. Es importante hacer la instalación en este orden, primero instalar R y, después, RStudio.
RStudio: https://www.rstudio.com
RStudio está dividido en cuatro paneles:
Interficie de RStudio
Consola (1): aquí es donde R espera que le demos instrucciones. Para ejecutarlas y obtener el resultado pulsamos «Intro». La primera vez que abramos RStudio esta sección ocupará toda la parte izquierda de la pantalla.
Scripts (2): trabajar en la consola es muy limitado ya que las instrucciones se han de introducir una a una. Lo habitual es trabajar con scripts o ficheros de instrucciones que podemos guardar y compartir. Para ver esta sección en la pantalla hay que abrir un script nuevo («File > New File > R Script»). Otra opción, aún mejor, es crear ficheros RMarkdown con código y narrativa.
Entorno, etc. (3) que incluye diversas pestañas: Environment (donde se irán registrando los objetos creados en la sesión de trabajo); History (se registran las instrucciones ejecutadas); Connections y Tutorial.
Ficheros, etc. (4) que incluye las pestañas Files, Plots, Packages, Help y Viewer. De momento, destacaremos la pestaña Packages, que proporciona un listado de los paquetes disponibles en R y los que han sido cargados en la sesión. A través de las opciones de esta pestaña podemos instalar nuevos paquetes o actualizar los existentes.
Lo más habitual será trabajar con scripts en el panel 2 de la figura. Allí escribiremos las funciones que queremos ejecutar. Asimismo, podemos añadir comentarios que, para diferenciarlos de las funciones, irán precedidos del signo # (para convertir un texto en comentario también es posible seleccionarlo y pulsar «Ctrl + Shift + C»). Aunque se pueden ejecutar varias funciones de una sola vez, normalmente ejecutaremos los script línea a línea. De esta manera, si se produce un error, sabremos inmediatamente en qué línea se encuentra. Para ejecutar una línea de código, situamos en ella el cursor y pulsamos «Run» o «Ctrl + Intro».
A menudo utilizaremos paquetes (libraries) que extienden la configuración básica de R. Para utilizar un paquete primero tendremos que instalarlo (esto lo haremos una única vez) y llamarlo (esto habrá que hacerlo cada vez que abramos una sesión en R). Es recomendable hacer la llamada al paquete al inicio del script. Por ejemplo, estos son los dos comandos para instalar y ejecutar el paquete rworldmap:
install.packages("rworldmap")
library(rworldmap)A medida que trabajemos con R iremos acumulando ficheros (datos, scripts, resultados de los análisis, figuras, etc.). Para mantener organizada la información, es recomendable trabajar con «proyectos» (projects). Cada proyecto tiene su propio directorio, de manera que, al indicar un fichero que queremos utilizar, no es necesario detallar la ruta para acceder, sino que R lo buscará en el directorio del proyecto. En el caso de proyectos de envergadura, es conveniente crear subdirectorios de datos, scripts y outputs.
1.2 Elementos básicos de R
En su versión más simple, es posible utilizar R como una calculadora. Si introducimos en la consola 2+2, R nos devuelve el resultado.
2+2## [1] 4Al ejecutar el código, el resultado se muestra en la consola. Sin embargo, muchas veces será útil guardarlo en un objeto para poder reutilizarlo. En R el concepto de objeto es esencial. En R cualquier cosa puede ser un objeto: una cifra, una palabra, una tabla de datos, una figura, etc. Por tanto, un objeto puede ser algo muy simple o algo muy complejo.
Para crear un objeto, simplemente le damos un nombre y le asignamos un valor usando el operador <-
objeto_1 <- 48Acabamos de crear un objeto llamado objeto_1 y le hemos asignado el valor 48. Es importante evitar los espacios en blanco en los nombres de objetos. Generalmente se sustituyen por «_» o por «.» Tampoco podemos iniciar el nombre de un objeto con un dígito.
Ahora creamos un segundo objeto:
objeto_2 <- "buenos días"En el «Entorno» de RStudio podemos ver los objetos creados, sus valores y el tipo de objeto que RStudio les ha asignado. Si asignamos un nuevo valor a un objeto, R lo sustituye sin avisarnos por lo que hay que ser precavido para evitar borrar información involuntariamente:
objeto_2 <- 34Ahora tenemos dos objetos del mismo tipo y podemos combinarlos:
objeto_3 <- objeto_1 + objeto_2
print(objeto_3)## [1] 82Otro tipo de objeto son los vectores. Vamos a crear un vector concatenando 8 cifras:
vector_1 <- c(2, 3, 1, 6, 4, 3, 3, 7)Pedimos la media del vector y la guardamos en un nuevo objeto llamado media_vector_1:
media_vector_1 <- mean(vector_1)
print(media_vector_1)## [1] 3.625En los ejemplos anteriores hemos visto objetos que contenían dos tipos de valores: numéricos (numeric) y caracteres (character). R considera dos tipos más de información: enteros (integers, números enteros, sin decimales) y lógicos (logical, que adopta dos posibles valores: TRUE o FALSE).
Por lo que hace a las estructuras de datos, ya hemos visto los vectores (vectors). Además, R incluye datos en matrices (matrices o arrays), listas (lists) y dataframes. Este último tipo de estructura será la que utilizaremos con mayor asiduidad durante este curso. Un dataframe es un objeto bidimensional con filas y columnas. Es similar a una matriz, con la diferencia de que mientras la matriz sólo contiene datos de un tipo, los dataframes pueden combinar diferentes tipos de datos. Normalmente, en un dataframe cada fila corresponde a una observación individual y cada columna a una variable. En definitiva, un dataframe es similar a una hoja de cálculo en Microsoft Excel o en LibreOffice Calc.
A continuación se muestra un ejemplo de dataframe.
## Nombre Color Edad Peso
## 1 Amelia Amarillo 66 74.5
## 2 Daniel Rojo 25 69.1
## 3 Emma Verde 22 77.8
## 4 Lucas Rojo 30 71.1
## 5 Matilde Rojo 38 72.3
## 6 Pablo Verde 35 67.2A pesar de que es posible introducir datos directamente en R (por ejemplo, creando vectores), lo más habitual es importarlos desde una hoja de cálculo. Hay múltiples funciones para hacerlo. Una forma sencilla es guardar el fichero como un texto separado por tabulaciones e importarlo desde RStudio con la función read.table. En el ejemplo que viene a continuación, se ha importado desde R un fichero llamado «muestra.txt» con tres parámetros que indican que los datos tienen encabezamientos, están separados por tabulaciones y el símbolo que utilizamos para separar los decimales es la coma. Los datos importados se guardan en un objeto llamado muestra.
muestra <- read.table (file = "muestra.txt", header = TRUE, sep = "\t", dec = ",")Ahora podemos trabajar con el objeto muestra, obteniendo un resumen estadístico de sus datos:
summary(muestra)## Nombre Color Edad Peso
## Length:6 Length:6 Min. :22.00 Min. :67.20
## Class :character Class :character 1st Qu.:26.25 1st Qu.:69.60
## Mode :character Mode :character Median :32.50 Median :71.70
## Mean :36.00 Mean :72.00
## 3rd Qu.:37.25 3rd Qu.:73.95
## Max. :66.00 Max. :77.80También podemos obtener datos de un campo concreto (que se indica con el símbolo «$»), como la media de las edades:
mean(muestra$Edad)## [1] 36Como ya hemos visto, esta media puede guardarse en otro objeto para seguir trabajando con ella:
media_edad <- mean(muestra$Edad)
media_edad## [1] 36Para hacer un recuento de datos:
recuento <- table(muestra$Color)
recuento##
## Amarillo Rojo Verde
## 1 3 2Podríamos filtrar los registros que cumplen una determinada condición, ordenar los datos siguiendo algún criterio, etc.
2 Gráficos
2.1 Gráficos de base en R
R incorpora algunas funciones para crear gráficos como, por ejemplo, plot.
data(cars)
head(cars)## speed dist
## 1 4 2
## 2 4 10
## 3 7 4
## 4 7 22
## 5 8 16
## 6 9 10plot(cars)
La primera variable aparece representada en el eje horizontal (x) y la segunda en el vertical (y), pero, si lo deseamos, lo podemos modificar.
plot(cars$dist, cars$speed)
Es posible hacer gráficos de puntos (p), líneas (l) o combinar ambos (b). También podemos añadir color.
plot(cars, type = "p", col = "red")
Veamos ahora como etiquetar los ejes. El argumento “pch” permite modificar el tipo de marcador.
plot(cars, type = "p", col = "steelblue",
main = "Relación entre velocidad y distancia de frenado",
xlab = "Velocidad (km/h)",
ylab = "Distancia (metros)",
pch = 19)
Veamos ahora otros tipos de gráficos.
boxplot(cars, main = "Diagrama de caja")
hist(cars$speed) # Histograma
pie(recuento) # Gráfico de "tarta"
2.2 Gráficos con ggplot2
El paquete para gráficos ggplot2 se puede instalar por separado o conjuntamente con otros dentro de tidyverse.
library(tidyverse)Los gráficos ggplot tienen tres argumentos:
data: un dataframe
aesthetics (aes): variables (x, y), color, tamaño, forma, etc.
geometry (geom): tipo de gráfico (dispersión, barras, lineas, etc.)
Creamos un gráfico de dispersión al que denominamos «coches» y lo visualizamos.
coches <- ggplot(cars, aes(speed, dist)) +
geom_point() +
ggtitle("Diagrama de dispersión")
coches
Trabajaremos ahora con otro dataset también disponible en R.
data(mpg) # Cargamos los datos
str(mpg) # Observamos la estructura del dataset## tibble [234 x 11] (S3: tbl_df/tbl/data.frame)
## $ manufacturer: chr [1:234] "audi" "audi" "audi" "audi" ...
## $ model : chr [1:234] "a4" "a4" "a4" "a4" ...
## $ displ : num [1:234] 1.8 1.8 2 2 2.8 2.8 3.1 1.8 1.8 2 ...
## $ year : int [1:234] 1999 1999 2008 2008 1999 1999 2008 1999 1999 2008 ...
## $ cyl : int [1:234] 4 4 4 4 6 6 6 4 4 4 ...
## $ trans : chr [1:234] "auto(l5)" "manual(m5)" "manual(m6)" "auto(av)" ...
## $ drv : chr [1:234] "f" "f" "f" "f" ...
## $ cty : int [1:234] 18 21 20 21 16 18 18 18 16 20 ...
## $ hwy : int [1:234] 29 29 31 30 26 26 27 26 25 28 ...
## $ fl : chr [1:234] "p" "p" "p" "p" ...
## $ class : chr [1:234] "compact" "compact" "compact" "compact" ...Creamos un gráfico de dispersión con tres variables. Representamos las dos primeras en los dos ejes y la tercera mediante el color (colour) o la forma (shape).
fig1 <- ggplot(mpg, aes(displ, hwy, colour = class)) +
geom_point()
fig1
Realizamos ahora un gráfico de barras.
fig2 <- ggplot(mpg, aes(x=manufacturer)) +
geom_bar(stat = "count")
fig2
Introducimos algunos cambios en el formato.
fig3 <- fig2 + geom_bar(fill="steelblue") + theme_minimal()
fig3
Convertimos las columnas en horizontales.
fig4 <- fig3 + coord_flip()
fig4
Acabaremos este apartado con un mapa de color (heatmap) que muestra el retardo mensual de algunas compañía aéreas. Estos son los datos que hemos recopilado:
delays <- read.table("delays.txt", header = TRUE, sep = "\t", dec = ",")
delays| Companyia | Gener | Febrer | Març | Abril | Maig | Juny | Juliol | Agost | Setembre | Octubre | Novembre | Desembre |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Alitalia | 28 | 30 | 26 | 20 | 20 | 32 | 15 | 12 | 6 | 2 | 3 | 16 |
| Easyjet | 118 | 176 | 121 | 75 | 143 | 5 | 8 | 163 | 145 | 173 | 148 | 145 |
| Emirates | 20 | 219 | 211 | 23 | 123 | 50 | 60 | 90 | 150 | 200 | 240 | 308 |
| Iberia | 200 | 490 | 400 | 390 | 440 | 120 | 360 | 370 | 500 | 310 | 250 | 475 |
| KLM | 9 | 10 | 7 | 5 | 8 | 14 | 7 | 5 | 1 | 2 | 4 | 7 |
| Lufthansa | 365 | 247 | 180 | 113 | 117 | 177 | 330 | 382 | 348 | 345 | 100 | 184 |
| Ryanair | 432 | 454 | 431 | 381 | 367 | 367 | 318 | 467 | 432 | 455 | 424 | 415 |
| Vueling | 117 | 177 | 330 | 382 | 348 | 345 | 100 | 184 | 365 | 247 | 180 | 345 |
Cambiamos el formato de los datos para elaborar el mapa de calor:
library(reshape)
delays2 <- reshape::melt(delays)
head(delays2)## Companyia variable value
## 1 Alitalia Gener 28
## 2 Easyjet Gener 118
## 3 Emirates Gener 20
## 4 Iberia Gener 200
## 5 KLM Gener 9
## 6 Lufthansa Gener 365Y elaboramos el gráfico:
p <- ggplot(delays2, aes(variable, Companyia, fill=value)) +
geom_tile()
print(p)
Ordenamos alfabéticamente las compañías y mejoramos el formato.
# Primero la ordenación
delays2$Companyia <- with(delays2,factor(Companyia,levels = rev(sort(unique(Companyia)))))
# A continuación, repetimos el gráfico, ahora ya ordenado
p2 <- ggplot(delays2, aes(variable, Companyia, fill= value)) +
geom_tile()
# Añadimos algunos títulos
p3 <- p2+xlab("Mes del año") +
ylab("Compañía aérea") +
ggtitle("Retraso mensual") +
labs(fill = "Media en minutos")
# Por último, cambiamos los colores
p3+scale_fill_gradient(low="white",high="steelblue")
2.3 Gráficos interactivos
Para dotar a los gráficos de interactividad, utilizaremos el paquete plotly.
library(plotly)
fig1_interactiva <- plotly::ggplotly(fig1)
fig1_interactivaVeamos otro ejemplo.
fig3_int <- plotly::ggplotly(fig3)
fig3_intAhora crearemos un gráfico de líneas desde cero.
bitcoin <- read.table("bitcoin.txt", header = TRUE) # Importamos el dataset
bitcoin$date <- as.Date(bitcoin$date) # Cambiamos el formato del campo fecha
head(bitcoin) # Visualizamos los primeros registros## date value
## 1 2013-04-28 135.98
## 2 2013-04-29 147.49
## 3 2013-04-30 146.93
## 4 2013-05-01 139.89
## 5 2013-05-02 125.60
## 6 2013-05-03 108.13fig5 <- ggplot(bitcoin, aes(x=date, y=value)) + # Creamos el gráfico
geom_line() +
xlab("Evolución del precio del bitcoin") +
ylab("Dólares") +
theme_minimal()
fig5_int <- plotly::ggplotly(fig5) # Hacemos el gráfico interactivo
fig5_int # Visualizamos el gráfico2.4 Gráficos animados
Utilizaremos el paquete gganimate y diversos datasets disponibles en R.
library(gganimate) # Llamamos al paquete
library(babynames) # Cargamos los datos
library(tidyverse)
data(babynames)
# Filtramos las frecuencias relativas a tres nombres y vemos el formato
nombres <- babynames %>%
filter(name %in% c("Ashley", "Patricia", "Helen")) %>%
filter(sex=="F")
head(nombres)## # A tibble: 6 x 5
## year sex name n prop
## <dbl> <chr> <chr> <int> <dbl>
## 1 1880 F Helen 636 0.00652
## 2 1881 F Helen 612 0.00619
## 3 1882 F Helen 838 0.00724
## 4 1883 F Helen 862 0.00718
## 5 1884 F Helen 986 0.00717
## 6 1884 F Patricia 6 0.0000436# Los representamos en un gráfico de lineas
fig_nombres <- ggplot(nombres, aes(x=year, y=n, group=name, color=name)) +
geom_line() +
geom_point() +
ggtitle("Popularidad de tres nombres en Estados Unidos") +
xlab("Año") +
ylab("Número de bebés")
fig_nombres + gganimate::transition_reveal(year) # Animamos el gráfico
Podemos guardar el gráfico animado en formato GIF en nuestro directorio de trabajo:
anim_save(nombre.gif)Veamos otro ejemplo en el que utilizaremos el paquete gifski.
library(gifski)
library(gapminder)
data(gapminder)
str(gapminder)## tibble [1,704 x 6] (S3: tbl_df/tbl/data.frame)
## $ country : Factor w/ 142 levels "Afghanistan",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ continent: Factor w/ 5 levels "Africa","Americas",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ year : int [1:1704] 1952 1957 1962 1967 1972 1977 1982 1987 1992 1997 ...
## $ lifeExp : num [1:1704] 28.8 30.3 32 34 36.1 ...
## $ pop : int [1:1704] 8425333 9240934 10267083 11537966 13079460 14880372 12881816 13867957 16317921 22227415 ...
## $ gdpPercap: num [1:1704] 779 821 853 836 740 ...# Hacemos un gráfico con ggplot2
figura10 <- ggplot(gapminder, aes(gdpPercap, lifeExp, size = pop, colour = country)) +
geom_point(alpha = 0.7, show.legend = FALSE) +
scale_colour_manual(values = country_colors) +
scale_size(range = c(2, 12)) +
scale_x_log10() +
# Esta parte es la de la animación:
labs(title = 'Year: {frame_time}', x = 'GDP per capita', y = 'Esperanza de vida') +
transition_time(year) +
ease_aes('linear')
figura10 # Visualizamos el gráfico
3 Mapas
3.1 Mapas con rworldmap
Comenzaremos haciendo un mapa en rworldmap. Utilizaremos un dataset con datos mundiales disponible en R.
library(rworldmap) # Este es el paquete que utilizaremos
data(countryExData) # Y estos los datos
str(countryExData)## 'data.frame': 149 obs. of 80 variables:
## $ ISO3V10 : chr "AGO" "ALB" "ARE" "ARG" ...
## $ Country : chr "Angola" "Albania" "United Arab Emirates " "Argentina" ...
## $ EPI_regions : chr "Sub-Saharan Africa" "Central and Eastern Europ" "Middle East and North Africa" "Latin America and Caribbe" ...
## $ GEO_subregion : chr "Southern Africa" "Central Europe" "Arabian Peninsula" "South America" ...
## $ Population2005 : num 15941 3130 4496 38747 3016 ...
## $ GDP_capita.MRYA : num 2314 4955 22698 13652 5011 ...
## $ landlock : int 0 0 0 0 1 0 1 1 1 0 ...
## $ landarea : num 1251896 28346 74777 2736296 28273 ...
## $ density : num 0.2 34.3 8.7 1.3 30.3 0.3 16.3 14.6 91 71.2 ...
## $ EPI : num 39.5 84 64 81.8 77.8 79.8 89.4 72.2 54.7 78.4 ...
## $ ENVHEALTH : num 8.9 89.3 89.8 91.1 88 99.3 98.1 76.4 37.6 98.8 ...
## $ ECOSYSTEM : num 70.1 78.6 38.2 72.5 67.5 60.4 80.7 67.9 71.7 58 ...
## $ ENVHEALTH.1 : num 8.9 89.3 89.8 91.1 88 99.3 98.1 76.4 37.6 98.8 ...
## $ AIR_E : num 49.2 99.1 85.1 87.3 99.4 84.9 97 97.7 99.5 50.2 ...
## $ WATER_E : num 61.6 96.5 27.1 74.9 28 62.5 79.9 48.5 62.8 52.3 ...
## $ BIODIVERSITY : num 58.9 4 36.6 33.6 16 78.1 71.6 29 62.5 10 ...
## $ PRODUCTIVE_NATURAL_RESOURCES: num 81.3 79.4 74.1 71.5 82.1 91.8 88.2 85.7 48 76.1 ...
## $ CLIMATE : num 74.6 93.4 26.6 82.3 87.2 42.5 79.9 77.1 81.5 69.5 ...
## $ DALY_SC : num 0 99.5 98.9 98 98.2 99.6 99.8 93 26.1 99.6 ...
## $ WATER_H : num 19.8 91.3 98.8 91.3 83.3 100 100 53.6 44.7 100 ...
## $ AIR_H : num 16 66.8 62.4 76.9 72.5 97.9 92.8 66.2 53.5 96 ...
## $ AIR_E.1 : num 49.2 99.1 85.1 87.3 99.4 84.9 97 97.7 99.5 50.2 ...
## $ WATER_E.1 : num 61.6 96.5 27.1 74.9 28 62.5 79.9 48.5 62.8 52.3 ...
## $ BIODIVERSITY.1 : num 58.9 4 36.6 33.6 16 78.1 71.6 29 62.5 10 ...
## $ FOREST : num 95.4 100 100 75.9 70.1 100 100 100 0 100 ...
## $ FISH : num 87.3 62.5 50 58.8 NA 96.7 NA NA NA 47.4 ...
## $ AGRICULTURE : num 61.3 75.6 72.3 79.9 94.2 78.7 76.4 71.4 95.9 80.8 ...
## $ CLIMATE.1 : num 74.6 93.4 26.6 82.3 87.2 42.5 79.9 77.1 81.5 69.5 ...
## $ ACSAT_pt : num 19.3 89.5 97.7 89.5 80.1 100 100 46.2 25.1 100 ...
## $ WATSUP_pt : num 20.2 93.2 100 93.2 86.4 100 100 61 64.3 100 ...
## $ DALY_pt : num 0 99.5 98.9 98 98.2 99.6 99.8 93 26.1 99.6 ...
## $ INDOOR_pt : num 0 47.4 94.7 94.7 72.2 94.7 94.7 48.4 0 94.7 ...
## $ PM10_pt : num 40 70.1 11.2 51.3 59 100 87.8 67 84.1 95.4 ...
## $ OZONE_H_pt : num 0 99.1 100 92.4 100 100 99.2 100 99.4 99.7 ...
## $ SO2_pt : num 98.4 98.5 70.2 98.8 98.8 69.9 94.4 95.4 99.3 0.6 ...
## $ OZONE_E_pt : num 0 99.8 100 75.7 100 100 99.6 100 99.6 99.8 ...
## $ WATQI_pt : num 29.4 93 0 76.4 31.7 75.3 59.8 31.7 25.6 59.6 ...
## $ WATSTR_pt : num 98.3 90.3 100 100 100 73.4 100 100 63 100 ...
## $ WATQI_GEMS.station.data : num NA 93 NA 76.4 NA 75.3 59.8 NA NA 59.6 ...
## $ FORGRO_pt : num 95.4 100 100 75.9 70.1 100 100 100 0 100 ...
## $ CRI_pt : num 99.7 5.5 100 39.8 37.7 86.1 80.1 46.2 84.1 9.6 ...
## $ EFFCON_pt : num 95.7 1.6 2.3 33.9 10.4 79 63 11.9 40.9 11.5 ...
## $ AZE_pt : num 0 NA NA 40 0 69.4 NA NA NA NA ...
## $ MPAEEZ_pt : num 14 6 1 2 100 78 100 100 100 0 ...
## $ EEZTD_pt : num 74.5 25.1 0 17.5 NA 93.5 NA NA NA 0 ...
## $ MTI_pt : num 100 100 100 100 NA 100 NA NA NA 94.9 ...
## $ IRRSTR_pt : num 97.5 100 51.8 74.6 97 50.7 100 82.9 100 100 ...
## $ AGINT_pt : num 100 90.2 100 78.4 94.5 79.6 63.2 91.1 92 87.1 ...
## $ AGSUB_pt : num 100 100 100 100 100 99.9 22.8 100 100 22.8 ...
## $ BURNED_pt : num 0 78.9 96.1 55.7 79.5 63.3 96 78.4 87.7 98.6 ...
## $ PEST_pt : num 9.1 9.1 13.6 90.9 100 100 100 4.5 100 95.5 ...
## $ GHGCAP_pt : num 65.8 98.8 38.6 87.1 98 45.4 81.6 88.7 94 77.7 ...
## $ CO2IND_pt : num 95 85 32.1 92.7 78.3 76.2 82.3 97.1 100 59.7 ...
## $ CO2KWH_pt : num 63 96.3 9 67 85.1 5.9 75.7 45.6 50.5 71.1 ...
## $ ACSAT : num 31 91 98 91 83 100 100 54 36 100 ...
## $ WATSUP : num 53 96 100 96 92 100 100 77 79 100 ...
## $ DALY : num 109 0.3 0.6 1.1 1 0.2 0.1 3.9 41 0.2 ...
## $ INDOOR : num 95 50 5 5 26.4 5 5 49 95 5 ...
## $ PM10 : num 91.4 55.5 125.6 77.9 68.7 ...
## $ OZONE_H : num 4948.8 15.8 0 140.4 0 ...
## $ SO2 : num 0.7 0.6 12.6 0.5 0.5 12.7 2.4 1.9 0.3 41.9 ...
## $ OZONE_E : num 1.36e+09 6.81e+05 2.63e+01 9.96e+07 0.00 ...
## $ WATQI : num 57.5 95.8 39.9 85.8 58.9 85.2 75.9 58.9 55.3 75.7 ...
## $ WATQI_GEMS.station.data.1 : num NA 95.8 NA 85.8 NA 85.2 75.9 NA NA 75.7 ...
## $ WATSTR : num 5.5 0 41.6 24.1 68.6 45.7 0 31.4 0 49.8 ...
## $ FORGRO : num 1 1 1 0.9 0.9 1 1.1 1 0.6 1.1 ...
## $ CRI : num 0.5 0 0.5 0.2 0.2 0.4 0.4 0.2 0.4 0 ...
## $ EFFCON : num 9.6 0.2 0.2 3.4 1 7.9 6.3 1.2 4.1 1.2 ...
## $ AZE : num 0 NA NA 40 0 69.4 NA NA NA NA ...
## $ MPAEEZ : num 1.4 0.6 0.1 0.2 10 7.8 10 10 10 0 ...
## $ EEZTD : num 0.255 0.749 1 0.825 NA ...
## $ MTI : num 0.0016 0 0.0034 0.0044 NA 0.0014 NA NA NA -0.001 ...
## $ IRRSTR : num 2.2 0 41 21.6 2.5 41.9 0 14.6 0 0 ...
## $ AGINT : num 0 6.2 0 13.7 3.5 12.9 23.3 5.6 5.1 8.2 ...
## $ AGSUB : num 0 0 0 0 0 0 36 0 0 36 ...
## $ BURNED : num 15.3 2.9 0.5 6 2.8 5 0.5 2.9 1.7 0.2 ...
## $ PEST : num 2 2 3 20 22 22 22 1 22 21 ...
## $ GHGCAP : num 20 2.9 34.1 8.9 3.3 30.5 11.8 8.1 5.3 13.8 ...
## $ CO2IND : num 1.2 1.9 5.5 1.4 2.3 2.5 2.1 1.1 0.8 3.6 ...
## $ CO2KWH : num 343 34 844 306 138 873 225 505 459 268 ...Creamos un objeto combinando nuestro dataset con el mapa. Haremos el cruce por dos campos diferentes para comprobar cual funciona mejor. En primer lugar, por el nombre del país. El argumento verbose nos devuelve los errores detectados.
p1 <- rworldmap::joinCountryData2Map(countryExData,
joinCode = "NAME",
nameJoinColumn = "Country",
verbose = TRUE)## 145 codes from your data successfully matched countries in the map
## 4 codes from your data failed to match with a country code in the map
## failedCodes failedCountries
## [1,] NA "Dem. Rep. Congo"
## [2,] NA "Czech Rep."
## [3,] NA "Dominican Rep."
## [4,] NA "Trinidad & Tobago"
## 98 codes from the map weren't represented in your dataEn segundo lugar, haremos el cruce por el campo ISO3.
p2 <- rworldmap::joinCountryData2Map(countryExData,
joinCode = "ISO3",
nameJoinColumn = "ISO3V10",
verbose = TRUE)## 149 codes from your data successfully matched countries in the map
## 0 codes from your data failed to match with a country code in the map
## failedCodes failedCountries
## 94 codes from the map weren't represented in your dataDado que la segunda estrategia ha funcionado mejor, será éste el objeto que utilizaremos para crear un mapamundi con información de la densidad de población por países.
mapa <- rworldmap::mapCountryData(p2, nameColumnToPlot="density")
Veamos algunos argumentos adicionales.
mapa2 <- rworldmap::mapCountryData(p2, nameColumnToPlot="density",
catMethod=c(0,20,40,60,80,100), # Igualar la franjas de la leyenda
mapTitle = "Densidad de población", # Título del gráfico
oceanCol='lightblue', # Océanos en color azul
missingCountryCol='white') # Países sin datos en blanco
3.2 Mapas con ggplot
En el apartado sobre gráficos utilizamos el paquete ggplot. Veremos ahora como emplearlo para crear gráficos.
library(tidyverse) # tidyverse contiene el paquete ggplot
ggplot() + borders("world")
Ahora ubicaremos en el mapa algunas ciudades mediante sus coordenadas geográficas
ciudades <- read.table("ciudades.txt", header = TRUE, sep = "\t")
head(ciudades)## ciudad lat lon
## 1 Madrid 40.41650 -3.70256
## 2 Barcelona 41.38879 2.15899
## 3 Oporto 41.16110 -8.62916
## 4 Lisboa 38.72960 -9.14235
## 5 Vitoria 42.85020 -2.67290ggplot() + borders("world", c("spain", "portugal"), fill = "white") +
geom_point(data = ciudades, aes(x = lon, y = lat)) +
geom_text(data = ciudades, aes(x = lon, y = lat, label = ciudad))
3.3 Mapa de comunidades autónomas
Necesitaremos los paquetes rgdal y broom, así como los datos geográficos de las comunidades autónomas, que se pueden descargar desde: https://www.arcgis.com/home/item.html?id=5f689357238847bc823a2fb164544a77
library(rgdal)
library(broom)A continuación trabajamos con los datos del fichero SHP
shapefile_ccaa <- rgdal::readOGR("Comunidades_Autonomas_ETRS89_30N.shp") # Leemos los datos## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\borrego\OneDrive - Universitat de Barcelona\CURSOS\CURS_R_2022\Curs R\Comunidades_Autonomas_ETRS89_30N.shp", layer: "Comunidades_Autonomas_ETRS89_30N"
## with 19 features
## It has 3 fieldsdata_ccaa <- broom::tidy(shapefile_ccaa) # Los convertimos en un dataframe
ggplot(data_ccaa, aes(x= long, y = lat, group = group)) + # Hacemos el mapa
geom_polygon(fill = "violetred4", color = "white") +
theme_minimal() +
theme(axis.line = element_blank(),
axis.text = element_blank(),
axis.title = element_blank(),
axis.ticks = element_blank())
En la página web del Instituto Nacional de Estadística (https://www.ine.es/prensa/eb_2018.pdf) hemos obtenido la cifra de bibliotecas por Comunidad Autónoma. Importamos los datos y los representamos gráficamente.
biblios <- read.table("bibliotecas.txt", header = TRUE, sep = "\t")
biblios$id <- as.character(biblios$id)
bibliotecas_grafico <- data_ccaa %>%
left_join(biblios, by= "id")
graf1 <- ggplot(bibliotecas_grafico, aes(x = long, y = lat, group = group)) +
geom_polygon(aes(fill=Bibliotecas), color= "white", size = 0.2) +
scale_fill_distiller(palette = "Blues", direction = 1) +
labs( title = "Bibliotecas públicas por Comunidades Autónomas",
subtitle = "Unidades: Número de bibliotecas",
caption = "Fuente: INE",
fill = "Bibliotecas públicas") +
theme_minimal() +
theme(axis.line = element_blank(),
axis.text = element_blank(),
axis.title = element_blank(),
axis.ticks = element_blank())
graf1
Repetiremos el gráfico sin los datos de las Islas Canarias para que quede centrado en la península y lo haremos interactivo con plotly.
bibliotecas_grafico_sc <- bibliotecas_grafico[bibliotecas_grafico$id != 4, ]
graf_sc <- ggplot(bibliotecas_grafico_sc, aes(x = long, y = lat, group = group)) +
geom_polygon(aes(fill=Bibliotecas), color= "white", size = 0.2) +
scale_fill_distiller(palette = "Blues", direction = 1) +
labs( title = "Bibliotecas públicas por Comunidades Autónomas",
subtitle = "Unidades: Número de bibliotecas",
caption = "Fuente: INE",
fill = "Bibliotecas públicas") +
theme_minimal() +
theme(axis.line = element_blank(),
axis.text = element_blank(),
axis.title = element_blank(),
axis.ticks = element_blank())
graf_sc_int <- plotly::ggplotly(graf_sc)
graf_sc_int3.4 Mapa con OpenStreetMap
Acabaremos creando un mapa en el que ubicamos la Facultat.
library(leaflet)
library(dplyr)
leaflet() %>%
setView(lng=2.139076193999787, lat=41.381288320297756, zoom = 16) %>%
addTiles() %>%
addMarkers(lng=2.139076193999787, lat=41.381288320297756, popup="Estamos en la Facultat FIMA de la UB haciendo un curso de introducción a R")4 Extracción de datos de la web
En apartados anteriores hemos visto cómo importar ficheros para trabajar con los datos. Veremos ahora cómo obtenerlos desde la web.
Vamos a importar los datos de los 100 mejores libros de todos los tiempos según la web https://www.goodreads.com/list/show/1.Best_Books_Ever. Para ello, utilizaremos el paquete rvest.
library(rvest)Almacenamos el código HTML de la página de la que obtener los datos.
webpage <- rvest::read_html("https://www.goodreads.com/list/show/1.Best_Books_Ever")A continuación, extraemos los títulos de los libros. Para conocer las etiquetas, podéis utilizar una extensión del navegador como SelectorGadget (o la opción “Inspecciona” con el botón derecho del ratón).
titulos_html <- rvest::html_nodes(webpage,'.bookTitle')Convertimos los títulos de HTML a texto y los visualizamos.
titulos <- rvest::html_text(titulos_html)
head(titulos, 5)## [1] "\n The Hunger Games (The Hunger Games, #1)\n"
## [2] "\n Harry Potter and the Order of the Phoenix (Harry Potter, #5)\n"
## [3] "\n To Kill a Mockingbird\n"
## [4] "\n Pride and Prejudice\n"
## [5] "\n Twilight (The Twilight Saga, #1)\n"Vamos a eliminar el símbolo “|n” y los espacios iniciales.
titulos2 <- gsub("\n ","",titulos)
head(titulos2, 5)## [1] "The Hunger Games (The Hunger Games, #1)\n"
## [2] "Harry Potter and the Order of the Phoenix (Harry Potter, #5)\n"
## [3] "To Kill a Mockingbird\n"
## [4] "Pride and Prejudice\n"
## [5] "Twilight (The Twilight Saga, #1)\n"También eliminados el símbolo “|n” al final. Comprobamos que es correcto.
titulos3 <- gsub("\n","",titulos2)
head(titulos3, 5)## [1] "The Hunger Games (The Hunger Games, #1)"
## [2] "Harry Potter and the Order of the Phoenix (Harry Potter, #5)"
## [3] "To Kill a Mockingbird"
## [4] "Pride and Prejudice"
## [5] "Twilight (The Twilight Saga, #1)"Hacemos lo mismo con los autores.
autores_html <- rvest::html_nodes(webpage,'.authorName')
autores <- rvest::html_text(autores_html)
head(autores)## [1] "Suzanne Collins" "J.K. Rowling" "Harper Lee" "Jane Austen"
## [5] "Stephenie Meyer" "Markus Zusak"Y lo mismo con los ratings.
rating_html <- rvest::html_nodes(webpage,'.minirating')
rating <- rvest::html_text(rating_html)
head(rating)## [1] " 4.32 avg rating — 7,178,725 ratings"
## [2] " 4.50 avg rating — 2,852,308 ratings"
## [3] " 4.27 avg rating — 5,115,277 ratings"
## [4] " 4.27 avg rating — 3,475,155 ratings"
## [5] " 3.62 avg rating — 5,586,909 ratings"
## [6] " 4.38 avg rating — 2,107,898 ratings"Combinamos en un dataframe títulos, autores y ratings. Visualizamos los primeros registros.
libros <- data.frame(Titulo = titulos3, Autor = autores, Rating = rating)
head(libros)## Titulo Autor
## 1 The Hunger Games (The Hunger Games, #1) Suzanne Collins
## 2 Harry Potter and the Order of the Phoenix (Harry Potter, #5) J.K. Rowling
## 3 To Kill a Mockingbird Harper Lee
## 4 Pride and Prejudice Jane Austen
## 5 Twilight (The Twilight Saga, #1) Stephenie Meyer
## 6 The Book Thief Markus Zusak
## Rating
## 1 4.32 avg rating — 7,178,725 ratings
## 2 4.50 avg rating — 2,852,308 ratings
## 3 4.27 avg rating — 5,115,277 ratings
## 4 4.27 avg rating — 3,475,155 ratings
## 5 3.62 avg rating — 5,586,909 ratings
## 6 4.38 avg rating — 2,107,898 ratingsPara exportar el dataframe a un fichero, podemos usar el siguiente código. R guardará el fichero «Libros.txt» en nuestro directorio de trabajo.
write.table(libros, file = "Libros", sep = "\t")5 Extracción de datos de Twitter
Actualmente es muy común extraer datos de aplicaciones web como Twitter a través de APIs (Aplication Programming Interface), ya que la mayoría de organismos y empresas tienen una o varias.
En este apartado nos centraremos en la obtención de datos de Twitter a través de su API pública. Para ello utilizaremos el paquete rtweet.
library(rtweet)Lo primero que necesitarás para consultar la API será una cuenta en Twitter. Si no la tienes, tendrás que crear una.
A continuación, deberás solicitar una cuenta de desarrollador. Puedes hacerlo desde la dirección: https://developer.twitter.com/en/apply-for-access. Cuando tengas acceso, debes crear una nueva app. Una vez hayas llegado al final del proceso, tendrás las claves que necesitarás en los siguientes pasos para extraer los datos.
Comenzamos autenticándonos. Debemos introducimos tres datos: el nombre de la app, la API key y la API secret key.
appname <- "XXXXXXXXXX"
key <- "XXXXXXXXXX"
secret <- "XXXXXXXXXX"Prueba a escribir un tweet desde R. Al ejecutar la función, se abrirá una ventana del navegador para autenticarte. Si funciona correctamente recibirás el mensaje: Authentication complete. Please close this page and return to R.
post_tweet("Escribo este tweet desde R!")5.1 Usuarios a los que sigue una cuenta
Veamos ahora como obtener los identificadores de los usuarios a los que siguien una o varias cuentas, por ejemplo los que siguen la cuentas de la Facultat d’Informació i Mitjans Audiovisuals de la UB @FIMA_UB i la del Col·legi Oficial de Bibliotecaris i Documentalistes de Catalunya @cobdc.
following <- rtweet::get_friends(c("fima_ub", "cobdc"))5.2 Seguidores de una cuenta
De manera similar, se puede obtener el listado de seguidores (followers) de una cuenta:
followers <- rtweet::get_followers("fima_ub")5.3 Recuperación de tweets
Buscaremos ahora los tweets que mencionan diversas ciudades combinando sus nombres mediante operadores booleanos:
barcelona <- rtweet::search_tweets(q = "barcelona")
barcelona_y_madrid <- rtweet::search_tweets(q = "barcelona AND madrid")
barcelona_o_madrid <- rtweet::search_tweets(q = "barcelona OR madrid")Por defecto, la API pública de Twitter recupera los últimos 100 tweets. No obstante, este límite se puede ampliar hasta 18.000 tweets.
barcelona_2 <- rtweet::search_tweets(q = "barcelona", n = 500)Vamos a limitar la búsqueda a tweets en catalán, excluyendo los retweets:
rtweet::search_tweets("barcelona", lang = "ca", include_rts = FALSE)5.4 Otras búsquedas
Veamos ahora, de forma gráfica, la evolución cronológica por horas de los últimos 150 tweets publicados por la cuenta de “La Vanguardia”.
lavanguardia <- rtweet::get_timeline("lavanguardia", n = 150)
ts_plot(lavanguardia, "hours")También podemos observar la actividad en torno a tres tweets recientes publicados en la cuenta de la Facultad:
status_ids <- c("1375086465626091525", "1371816665550356480", "1370323578365829122")
rtweet::lookup_tweets(status_ids)De manera similar, podemos consultar la actividad reciente de algunos usuarios.
users <- c("lavanguardia", "fima_ub", "unibarcelona")
usr <- rtweet::lookup_users(users)Recuperamos ahora los últimos 500 tweets que mencionan la cuenta de la UB.
ub <- rtweet::search_tweets(q = "unibarcelona", n = 500)¿Quienes son los autores de estos tweets?
table(ub$screen_name)6 Consulta de fuentes bibliográficas
En este apartado vamos a consultar dos fuentes de datos bibliográficas que disponen de APIs que se pueden interpelar desde R: Crossref y Unpaywall.
6.1 Crossref
Crossref es una organización dedicada a enlazar metadatos de publicaciones académicas. Para interactuar con su API vamos a utilizar el paquete rcrossref. Lo primero, como siempre, es instalar el paquete y llamarlo.
library(rcrossref)La función cr_works() permite buscar documentos por DOI. Podemos buscar directamente algunos DOIs o guardarlos previamente en un objeto para buscarlos todos conjuntamente como haremos a continuación.
dois <- c("10.1038/s41598-019-57002-9",
"10.1155/2019/4146362",
"10.1186/s13100-019-0197-9",
"10.1038/s41598-019-56941-7",
"10.1177/1087054719894378")
referencias <- rcrossref::cr_works(dois)
head(referencias)## $meta
## NULL
##
## $data
## # A tibble: 5 x 35
## alternative.id container.title created deposited published.print
## <chr> <chr> <chr> <chr> <chr>
## 1 57002 Scientific Reports 2019-1~ 2021-12-~ 2019-12
## 2 4146362,4146362 Journal of Advance~ 2019-1~ 2019-12-~ 2019-12-31
## 3 197 Mobile DNA 2019-1~ 2020-12-~ 2019-12
## 4 56941 Scientific Reports 2019-1~ 2021-12-~ 2019-12
## 5 10.1177/1087054719894378 Journal of Attenti~ 2019-1~ 2021-06-~ 2021-07
## # ... with 30 more variables: published.online <chr>, doi <chr>, indexed <chr>,
## # issn <chr>, issue <chr>, issued <chr>, member <chr>, prefix <chr>,
## # publisher <chr>, score <chr>, source <chr>, reference.count <chr>,
## # references.count <chr>, is.referenced.by.count <chr>, subject <chr>,
## # title <chr>, type <chr>, update.policy <chr>, url <chr>, volume <chr>,
## # abstract <chr>, language <chr>, short.container.title <chr>,
## # assertion <list>, author <list>, link <list>, license <list>, ...
##
## $facets
## NULLLa siguiente es otra manera de hacer lo mismo con una única función.
referencias <- rcrossref::cr_works(dois=c("10.1002/leap.1347",
"10.1108/OIR-09-2019-0291",
"10.1353/pla.2020.0031",
"10.1016/j.acalib.2018.07.012"))La información contenida en el objeto referencias está organizada en tres partes (meta, data y facets), pero la más interesante es la denominada data. Usaremos la función pluck() del paquete purrr para extraer esta información y guardarla en un nuevo objeto. Para concatenar las funciones usaremos un pipe %>%. Se trata de un operador que toma el resultado de una acción para convertirlo inmediatamente en el input de la siguiente.
library(purrr)
library(tidyverse)
referencias2 <- cr_works(dois) %>%
purrr::pluck("data")De entre la información disponible, vamos a extraer tres campos: título, doi e issn.
referencias3 <- referencias2 %>%
dplyr::select(title, doi, issn)
head(referencias3)## # A tibble: 5 x 3
## title doi issn
## <chr> <chr> <chr>
## 1 Short- and medium-term impact of bariatric sur~ 10.1038/s41598-0~ 2045-2322
## 2 Operations Research for Transportation and Sus~ 10.1155/2019/414~ 0197-6729,2~
## 3 A benchmark of transposon insertion detection ~ 10.1186/s13100-0~ 1759-8753
## 4 Duplication and parallel evolution of the panc~ 10.1038/s41598-0~ 2045-2322
## 5 The Role of ADHD Symptomatology and Emotion Dy~ 10.1177/10870547~ 1087-0547,1~Veamos algunos ejemplos de búsquedas por campos guardando los resultados en un nuevo objeto.
ecologia <- cr_works(query = "ecologia")
urbano <- cr_works(flq = c(query.author = "urbano"))Y un par de ejemplos de consultas que nos devuelven el resultado en la consola.
cr_citation_count("10.1002/leap.1099")## doi count
## 1 10.1002/leap.1099 34cr_journals(issn = "1575-5886")## $meta
## NULL
##
## $data
## title
## 1 BiD textos universitaris de biblioteconomia i documentaci<U+FFFD>
## publisher issn last_status_check_time
## 1 Edicions de la Universitat de Barcelona 1575-5886 2022-01-10
## deposits_abstracts_current deposits_orcids_current deposits
## 1 TRUE FALSE TRUE
## deposits_affiliations_backfile deposits_update_policies_backfile
## 1 FALSE FALSE
## deposits_similarity_checking_backfile deposits_award_numbers_current
## 1 FALSE FALSE
## deposits_resource_links_current deposits_ror_ids_current deposits_articles
## 1 FALSE FALSE TRUE
## deposits_affiliations_current deposits_funders_current
## 1 FALSE FALSE
## deposits_references_backfile deposits_ror_ids_backfile
## 1 TRUE FALSE
## deposits_abstracts_backfile deposits_licenses_backfile
## 1 TRUE FALSE
## deposits_award_numbers_backfile deposits_open_references_backfile
## 1 FALSE TRUE
## deposits_open_references_current deposits_descriptions_current
## 1 TRUE FALSE
## deposits_references_current deposits_resource_links_backfile
## 1 FALSE FALSE
## deposits_descriptions_backfile deposits_orcids_backfile
## 1 FALSE FALSE
## deposits_funders_backfile deposits_update_policies_current
## 1 FALSE FALSE
## deposits_similarity_checking_current deposits_licenses_current
## 1 FALSE FALSE
## affiliations_current similarity_checking_current descriptions_current
## 1 0 0 0
## ror_ids_current funders_backfile licenses_backfile funders_current
## 1 0 0 0 0
## affiliations_backfile resource_links_backfile orcids_backfile
## 1 0 0 0
## update_policies_current open_references_backfile ror_ids_backfile
## 1 0 1 0
## orcids_current similarity_checking_backfile references_backfile
## 1 0 0 0.11875
## descriptions_backfile award_numbers_backfile update_policies_backfile
## 1 0 0 0
## licenses_current award_numbers_current abstracts_backfile
## 1 0 0 0.009375
## resource_links_current abstracts_current open_references_current
## 1 0 0.2941176 1
## references_current current_dois backfile_dois total_dois
## 1 0 51 320 371
##
## $facets
## NULL6.2 Unpaywall
Unpaywall es una extensión de navegador para localizar versiones en acceso abierto de artículos científicos. Para consultarlo desde R utilizaremos el paquete roadoi.
library(roadoi)Haremos una consulta sobre la disponibilidad de los artículos cuyos DOIs teníamos guardados en el objeto mis_dois. Aunque la consulta es libre, nos pide que nos identifiquemos con nuestro correo electrónico.
dois_oa <- roadoi::oadoi_fetch(dois = dois, email = "borrego@ub.edu")Consultamos los campos en los que está organizada la información que hemos recuperado. Para conocer el significado puedes consultar la página web de Unpaywall.
names(dois_oa)## [1] "doi" "best_oa_location" "oa_locations"
## [4] "oa_locations_embargoed" "data_standard" "is_oa"
## [7] "is_paratext" "genre" "oa_status"
## [10] "has_repository_copy" "journal_is_oa" "journal_is_in_doaj"
## [13] "journal_issns" "journal_issn_l" "journal_name"
## [16] "publisher" "published_date" "year"
## [19] "title" "updated_resource" "authors"Del objeto con la información podemos eliminar los documentos para los que no existe una versión en acceso abierto.
dois_oa2 <- dois_oa %>% dplyr::filter(is_oa == TRUE)
head(dois_oa2)## # A tibble: 4 x 21
## doi best_oa_location oa_locations oa_locations_emb~ data_standard is_oa
## <chr> <list> <list> <list> <int> <lgl>
## 1 10.1038/~ <tibble [1 x 10]> <tibble [6 ~ <tibble [0 x 0]> 2 TRUE
## 2 10.1155/~ <tibble [1 x 9]> <tibble [4 ~ <tibble [0 x 0]> 2 TRUE
## 3 10.1186/~ <tibble [1 x 10]> <tibble [5 ~ <tibble [0 x 0]> 2 TRUE
## 4 10.1038/~ <tibble [1 x 10]> <tibble [4 ~ <tibble [0 x 0]> 2 TRUE
## # ... with 15 more variables: is_paratext <lgl>, genre <chr>, oa_status <chr>,
## # has_repository_copy <lgl>, journal_is_oa <lgl>, journal_is_in_doaj <lgl>,
## # journal_issns <chr>, journal_issn_l <chr>, journal_name <chr>,
## # publisher <chr>, published_date <chr>, year <chr>, title <chr>,
## # updated_resource <chr>, authors <list>7 Clasificadores bayesianos
7.1 Estadística bayesiana
Muchos algoritmos de clasificación se basan en «entrenar» un modelo para que, a continuación, haga predicciones usando los ejemplos que ya han sido clasificados. Un modelo de este tipo puede servir, por ejemplo, para decidir si un correo electrónico es spam o no. Muchos de estos modelos se basan en la estadística bayesiana, un tipo de inferencia estadística en la que las evidencias u observaciones se emplean para actualizar o inferir la probabilidad de que una hipótesis sea cierta.
Comencemos con un ejemplo relativo a dos sucesos independientes. Si tiro una moneda al aire dos veces consecutivas, ¿cuál es la probabilidad de obtener dos caras? Dado que el resultado del primer lanzamiento (A) no influye sobre el resultado del segundo (B):
Posibles resultados de lanzar una moneda en dos ocasiones consecutivas
\[p (AB) = p (A) × p (B) = 1/2 × 1/2 = 1/4\]
La probabilidad conjunta de dos sucesos independientes (A y B) es la probabilidad de que se dé A multiplicada por la probabilidad de que ocurra B dado A y se representa con la siguiente fórmula:
\[p (AB) = p (A) × p (B | A)\]
En el caso del lanzamiento de la moneda, la probabilidad de que salgan dos caras es igual a la probabilidad de que salga cara \(p (A) = 0,5\), multiplicada por la probabilidad de que vuelva a salir cara si en el primer lanzamiento ha salido una cara \(p (B | A) = 0,5\). Por tanto, \(0,5 × 0,5 = 0,25\).
Dado que A y B son sucesos independientes, siguiendo la misma lógica podemos afirmar que la probabilidad conjunta \(p (AB)\) es la probabilidad de que se dé B multiplicada por la probabilidad de que ocurra A dado B:
\[p (AB) = p (B) × p (A | B)\]
Por tanto:
\[p (A) × p (B | A) = p (AB) = p (B) × p (A | B)\]
De lo anterior, se deriva el teorema de Bayes:
\[p (A | B) = \frac{p (B | A) × p (A)} {p (B)}\]
Veamos un ejemplo. Asumamos que la probabilidad de padecer cáncer es de 0,05. Es decir, el 5% de la población padece cáncer \(p (A) = 0,05\). Ahora asumamos que la probabilidad de ser fumador es de 0,10. Es decir, el 10% de la población fuma \(p (B) = 0,10\). Sabemos que el 20% de los pacientes con cáncer son fumadores, de manera que \(p (fumador | cáncer) = 0,20\). O, lo que es lo mismo, \(p (B | A) = 0,20\). ¿Cuál es la probabilidad de padecer cáncer si se es fumador? Aplicando el teorema de Bayes:
\[p(\text{cáncer|fumador}) = \frac{p (\text{fumador│cáncer}) × p (\text{cáncer})} {p (\text{fumador})} = \frac {0,20 × 0,05} {0,10} = 0,10\]
Es decir, puesto que la ratio de tabaquismo es el doble entre las personas que padecen cáncer que en el conjunto de la población, se duplica la probabilidad de padecer cáncer entre las personas que fuman.
Veamos un segundo ejemplo. Imaginemos un test de dopaje con una fiabilidad del 95%. Es decir, el test clasifica correctamente al 95% de los atletas que se dopan y al 95% de los que no. Supongamos que 1 de cada 50 atletas (el 2%) recurre al dopaje. Los resultados de aplicar el test a 1000 atletas se muestran en la figura.
Resultados de aplicar a 1000 atletas el test de dopaje
El test nos ofrecerá 19+49 = 68 positivos. Sin embargo, solo 19 corresponden realmente a atletas que se dopan (19/68 = 28%). El 72% restante son falsos positivos. A pesar del 95% de fiabilidad del test, la mayoría de los positivos son falsos.
Es importante no confundir la probabilidad de que un test arroje un resultado positivo dado que el atleta se dopa \(p (A | B) = 95\%\) con la probabilidad de que un atleta se dope dado un test positivo \(p (B | A) = 28\%\). Veamos el mismo resultado aplicando el teorema de Bayes:
\[p (B | A)= \frac{p (A | B) × p (B)} {p (A)}\]
\[p(\text{dopaje|positivo}) = \frac{p (\text{positivo│dopaje}) × p(\text{dopaje})} {p (\text{positivo})} = \frac{0,95 × 0,02} {0,068} = 0,28\]
El teorema de Bayes tiene múltiples aplicaciones en ciencia de datos, entre otras con finalidades de clasificación. Una de estas utilidades es la construcción de detectores de correo basura (spam). El sistema consiste en analizar las características de los correos spam: quién es el emisor, número de destinatarios, existencia de archivos adjuntos, errores ortográficos, palabras que aparecen en el mensaje, etc.
Nos centraremos ahora en la aparición de palabras, por ejemplo, los términos «barato» y «comprar». En una muestra de 100 mensajes, observamos que 25 son spam y 75 no lo son. De los 25 mensajes que sí son spam, la palabra «barato» aparece en 20 de ellos. Entre los 75 que no son spam, la palabra «barato» aparece en 5 de ellos.
\[p (\text{barato | spam}) = 20/25 = 0,80\]
\[p (\text{barato | no spam}) = 5/75 = 0,066\]
En el caso de la palabra «comprar», su probabilidad de aparición en un mensaje de spam es 15/25 y en un mensaje que no sea spam de 10/75.
\[p (\text{comprar | spam}) = 15/25 = 0,60\]
\[p (\text{comprar | no spam}) = 10/75 = 0,133\]
Si suponemos que la aparición de las dos palabras son sucesos independientes, la probabilidad de que aparezcan conjuntamente será el producto de las probabilidades individuales. Así, en el caso de los mensajes de spam, la probabilidad de que aparezcan juntas ambas palabras es la siguiente:
\[p (\text{barato y comprar | spam}) = 0,8 × 0,6 = 0,48\]
La probabilidad es de 0,48, es decir, que 12 de cada 25 mensajes de spam contendrán las palabras «barato» y «comprar» simultáneamente.
Siguiendo el mismo procedimiento, podemos calcular la frecuencia conjunta de ambas palabras en mensajes que no sean spam. Vemos que la probabilidad es muy inferior, no llega ni a un correo electrónico de los 75:
\[p (\text{barato y comprar | no spam}) = 0,066 × 0,133 = 0,009\]
Aplicando el teorema de Bayes ya podemos calcular la probabilidad de que un mensaje sea spam —o no lo sea— según la aparición de la palabra «barato», «comprar» o ambas.
\[p (\text{spam | barato}) = \frac{p (\text{barato │ spam}) × p(\text{spam})} {p (\text{barato})} = \frac{0,8 × 0,25} {0,25} = 0,80\]
\[p(\text{spam | comprar}) = \frac{p (\text{comprar │ spam}) × p(\text{spam})} {p (\text{comprar})} = \frac{0,6 × 0,25} {0,25} = 0,60\]
Finalmente, podemos calcular la posibilidad de que aparezcan las dos palabras conjuntamente:
\(p(\text{spam│barato y comprar})\)
\(= \frac{p (\text{barato | spam)} × p (comprar | spam) × p (spam)} {p(barato│spam) × p(comprar|spam) × p(spam) + p(barato|no spam) × p(comprar|no spam) × p(no spam)}\)
\(= \frac{(0,8 × 0,6 × 0,25)} {(0,8 × 0,6 × 0,25)+(0,066 ×0,133 ×0,75)} = 0,94\)
Es decir, según nuestro algoritmo, si un mensaje contiene las palabras «barato» y «comprar» hay un 94% de posibilidades de que sea spam (y 6% de que no lo sea).
Podemos ampliar la lista de palabras incluyendo, por ejemplo, la palabra «trabajo». Obviamente, a medida que añadimos palabras, las probabilidades de que aparezcan todas juntas se van reduciendo, pero es posible calcular la probabilidad de que aparezcan multiplicando las probabilidades de las palabras individuales. Supongamos que la palabra «trabajo» aparece en 5 de los mensajes que son spam y en 30 de los mensajes que no son spam de la muestra).
\(p (\text{trabajo | spam}) = 5/25 = 0,2\)
\(p (\text{trabajo | no spam}) = 30/75 = 0,4\)
Por tanto,
\(p (\text{spam │ barato y comprar y trabajo}) =\)
\(= \frac{0,8 × 0,6 × 0,2 × 0,25} {(0,8 × 0,6 × 0,2 × 0,25) + (0,066 × 0,133 × 0,4 × 0,75)} = 0,90\)
En este caso, según nuestro algoritmo, si un mensaje contiene simultáneamente las palabras «barato», «comprar» y «trabajo» hay un 90% de posibilidades de que sea spam —y un 10% de probabilidades de que no lo sea—. La razón de la reducción del porcentaje respecto al caso anterior es que «trabajo» es una palabra asociada a mensajes que no son de spam.
7.2 Clasificador bayesiano en R
Una escuela de aviación ha recopilado los datos del tiempo de 14 días y la indicación de si en cada uno de ellos se pudo volar o no. La siguiente tabla presenta los resultados. Con esta información quiere ser capaz de predecir, a partir de la previsión meteorológica, si será posible hacer prácticas de vuelo o no en un día determinado.
| Cielo | Temperatura | Humedad | Viento | Volar |
|---|---|---|---|---|
| sol | calor | alta | suave | no |
| sol | calor | alta | fuerte | no |
| nubes | calor | alta | suave | si |
| lluvia | templado | alta | suave | si |
| lluvia | frio | normal | suave | si |
| lluvia | frio | normal | fuerte | no |
| nubes | frio | normal | fuerte | si |
| sol | templado | alta | suave | no |
| sol | frio | normal | suave | si |
| lluvia | templado | normal | suave | si |
| sol | templado | normal | fuerte | si |
| nubes | templado | alta | fuerte | si |
| nubes | calor | normal | suave | si |
| lluvia | templado | alta | fuerte | no |
Comenzamos instalando y llamando al paquete e1071.
library(e1071)Importamos el fichero de datos con la información sobre el tiempo y las condiciones de vuelo. Cambiamos el formato de los datos de «character» a «factor».
vuelo <- read.table("vuelo.txt", header=T, sep = "\t")
library(dplyr)
vuelo <- vuelo%>%mutate_if(is.character, as.factor)El siguiente paso consiste en construir el modelo predictivo utilizando la función naiveBayes. Como argumentos, indicamos la variable que deseamos precedir y el origen de los datos.
VolarModel <- e1071::naiveBayes(Volar~., vuelo)Veamos las característica del modelo.
print(VolarModel)##
## Naive Bayes Classifier for Discrete Predictors
##
## Call:
## naiveBayes.default(x = X, y = Y, laplace = laplace)
##
## A-priori probabilities:
## Y
## no si
## 0.3571429 0.6428571
##
## Conditional probabilities:
## Cielo
## Y lluvia nubes sol
## no 0.4000000 0.0000000 0.6000000
## si 0.3333333 0.4444444 0.2222222
##
## Temperatura
## Y calor frio templado
## no 0.4000000 0.2000000 0.4000000
## si 0.2222222 0.3333333 0.4444444
##
## Humedad
## Y alta normal
## no 0.8000000 0.2000000
## si 0.3333333 0.6666667
##
## Viento
## Y fuerte suave
## no 0.6000000 0.4000000
## si 0.3333333 0.6666667Si nos fijamos, simplemente ha calculado la probabilidad global de volar (se pudo volar en 9 de los 14 días, es decir, 0,64) y de no volar (5 días de 14, 0,36), así como la probabilidad de cada uno de los valores de las variables en el modelo en función de que se haya podido volar o no (por ejemplo, la probabilidad de que llueva en un día de vuelo es de 0,33 porque de los 9 días que se ha podido volar, en 3 de ellos llovía).
\(p (\text{sí volar}) = 9/14 = 0,64\)
\(p (\text{no volar}) = 5/14 = 0,36\)
A partir de esta información, dado un nuevo caso, podemos predecir si en esas condiciones se podrá volar o no. Supongamos que la previsión para mañana es de tiempo soleado, temperatura fría, humedad alta y vientos fuertes. Por tanto:
\(p (sí | previsión)=\)
\(\frac{p (sol | sí) × p (frío | sí) × p (alta | sí) × p (fuerte | sí) × p (volar)} {[p (sol | sí) × p (frío | sí) × p (alta | sí) × p (fuerte | sí) × p (volar)] + [p (sol | no) × p (frío | no) × p (alta | no) × p (fuerte | no) × p (no.volar)]}=\)
\(\frac{0,22 × 0,33 × 0,33 × 0,33 × 0,64} {[0,22 × 0,33 × 0,33 × 0,33 × 0,64] + [0,60 × 0,20 × 0,80 × 0,60 × 0,36]}=\)
\(\frac{0,005} {0,005 + 0,020}=0,20\)
Concluimos que los más probable es que no se podrá volar.
Para ver el resultado en R, creamos un objeto con los datos de la previsión.
nuevo_caso <- data.frame(Cielo = "sol",
Temperatura = "frio",
Humedad = "alta",
Viento = "fuerte")Y pedimos al modelo que nos haga una predicción
predict(VolarModel, nuevo_caso)## [1] no
## Levels: no siSi queremos ver los porcentajes
predict(VolarModel, nuevo_caso, type = "raw")## no si
## [1,] 0.7954173 0.2045827Como vemos, la previsión es que no se podrá volar.
8 Clustering
En apartados anteriores hemos vimos cómo clasificar objetos en función de alguna característica de interés (por ejemplo, si se podrá volar o no en función de la previsión meteorológica). Ese tipo de tareas se denominan de segmentación «supervisada» porque, de entrada, ya tenemos una variable (en este caso la probabilidad de «volar») en función de la cual queremos clasificar los objetos.
En otras ocasiones quizá queramos encontrar grupos (por ejemplo, grupos de clientes), sin ninguna característica predefinida. Es decir, ¿pueden dividirse los clientes de manera natural en diversos grupos? Esta idea de agrupar los datos en categorías «no supervisadas» se denomina clustering. Existen dos métodos de clustering:
Jerárquicos:
Aglomerativos (bottom up): construyen los grupos por agregación a partir de los individuos considerados uno a uno.
Disociativos (top down): parten del conjunto total de individuos y los dividen en grupos más pequeños.
No jerárquicos: formación de nuevos clústers fusionando o dividiendo los existentes sin seguir un orden jerárquico. Por ejemplo, k-means.
8.1 Clustering jerárquico
La tabla muestra la frecuencia en la adquisición de cuatro productos por parte de cinco clientes. Utilizando esta información vamos a intentar crear segmentos de clientes que tengan un comportamiento similar.
| Producto_1 | Producto_2 | Producto_3 | Producto_4 | |
|---|---|---|---|---|
| Cliente_1 | 2 | 3 | 4 | 6 |
| Cliente_2 | 1 | 3 | 5 | 4 |
| Cliente_3 | 5 | 6 | 4 | 6 |
| Cliente_4 | 8 | 8 | 9 | 9 |
| Cliente_5 | 7 | 9 | 8 | 8 |
A partir de la tabla de datos, calculamos la distancia entre todos los pares de individuos utilizando la distancia euclidiana.
\[d(X, Y) = \sqrt{(X_{b} - X_{a})^2 + (Y_{b} - Y_{a})^2}\]
Por ejemplo, la distancia entre los clientes 1 y 2 será la siguiente:
\[d(C_{1}, C_{2}) = \sqrt{(2 - 1)^2 + (3 - 3)^2 + (4 - 5)^2 + (6 - 4)^2} = \sqrt{6} = 2,45\]
El resultado se muestra en una matriz simétrica.
data_dist = dist(data, method = "euclidean")
data_dist## Cliente_1 Cliente_2 Cliente_3 Cliente_4
## Cliente_2 2.449490
## Cliente_3 4.242641 5.477226
## Cliente_4 9.746794 10.723805 6.855655
## Cliente_5 9.000000 9.848858 5.744563 2.000000La distancia más reducida es la existente entre los clientes 4 y 5, que conforman el primer cluster con una distancia entre ellos de 2.
Primer cluster
Ahora hemos de recalcular la tabla de distancias. Las distancias entre los clientes 1, 2 y 3 continúan siendo las mismas, pero hemos de determinar la distancia entre estos clientes y el cluster {4,5}. Hay diversas maneras de hacerlo. Una forma es la denominada single link o nearest neigbor, según la cual elegimos la distancia al elemento más cercano dentro del cluster. Por ejemplo, la distancia entre el cliente 1 y el cluster {4,5} será de 9, la distancia que había entre los clientes 1 y 5, porque la distancia entre los clientes 1 y 4 era mayor, de 9,7.
Matrices de distancias
La gráfica muestra los sucesivos recalculos de las distancias hasta agrupar a todos los clientes. El código para hacerlo con R y generar el dendograma es el siguiente.
cluster_clients = hclust(data_dist, method="single")
plot(cluster_clients, hang=-1)
Veamos un segundo ejemplo. En este caso utilizaremos un dataset que viene en R con datos de delicuencia en Estados Unidos.
data("USArrests")
summary(USArrests)## Murder Assault UrbanPop Rape
## Min. : 0.800 Min. : 45.0 Min. :32.00 Min. : 7.30
## 1st Qu.: 4.075 1st Qu.:109.0 1st Qu.:54.50 1st Qu.:15.07
## Median : 7.250 Median :159.0 Median :66.00 Median :20.10
## Mean : 7.788 Mean :170.8 Mean :65.54 Mean :21.23
## 3rd Qu.:11.250 3rd Qu.:249.0 3rd Qu.:77.75 3rd Qu.:26.18
## Max. :17.400 Max. :337.0 Max. :91.00 Max. :46.00Generamos un dendrograma jerárquico aglomerativo con método de agrupación average. En este caso la distancia de un individuo a un cluster no es la distancia al individuo más cercano, sino la distancia media a todos los miembros del cluster. A diferencia del ejemplo anterior, en este caso no hemos calculado las distancias por separado, sino que todo el proceso está agrupado en una única función.
crim <- hclust(dist(USArrests), "average")
plot(crim, hang=-1)
Veamos, por ejemplo, los datos de Iowa y New Hampshire que aparecen como dos estados muy similares en la parte derecha del dendrograma.
## Murder Assault UrbanPop Rape
## Iowa 2.2 56 57 11.3
## New Hampshire 2.1 57 56 9.5O los datos de Illinois y New York, que también aparecen como estados similares.
## Murder Assault UrbanPop Rape
## Illinois 10.4 249 83 24.0
## New York 11.1 254 86 26.1Dividimos el dendrograma en cuatro clusters y los marcamos en color rojo.
plot(crim)
rect.hclust(crim, k = 4, border = "red")
8.2 Clustering no jerárquico (k-means)
Veremos un ejemplo de creación de clusters no jerárquicos con el método k-means. Utilizaremos el paquete cluster sobre un dataset disponible en R.
library(cluster) # Llamamos al paquete
data(ruspini) # Descargamos los datos
plot(ruspini) # Representamos los datos en un diagrama de dispersión
set.seed(7) # Inicio aleatorio
k <- kmeans(ruspini, centers = 4) # Creamos 4 clusters (podemos optar por otra cantidad)
g1 <- plot(ruspini, col = k$cluster) # Representamos los clusters con diversos colores
points(k$centers, cex=2, col=11, pch=19) # Representamos los centroides
Utilizaremos de nuevo los datos de delicuencia de un caso anterior para ilustrar un segundo ejemplo de clustering no jerárquico.
k2 <- kmeans(scale(USArrests), 4)
# Medias de las variables en cada cluster
aggregate(USArrests, by=list(cluster=k2$cluster), mean)## cluster Murder Assault UrbanPop Rape
## 1 1 3.60000 78.53846 52.07692 12.17692
## 2 2 5.65625 138.87500 73.87500 18.78125
## 3 3 10.81538 257.38462 76.00000 33.19231
## 4 4 13.93750 243.62500 53.75000 21.41250# Representación gráfica
library(factoextra)
fviz_cluster(k2, data = USArrests,
palette = c("#00AFBB","#2E9FDF", "#E7B800", "#FC4E07", "#E495A5"),
main = "Agrupación de estados por criminalidad")