Todas las entradas de: rvaquerizo

Etiquetas en scatter plot. Muertes covid por millón de habitantes vs gasto en salud

He intentado permanecer ajeno a los datos del coronoavirus pero es imposible, sin embargo, me gustaría aprovechar para mostrar algunos trucos con R y Python. Esta vez en una sola entrada vamos a tratar las siguientes situaciones:

  • Importar la tabla de worldometer sobre datos de países.
  • Problemas con la librería OECD.
  • Importar Excel con R.
  • Not in con R.
  • Gráficos de ranking ordenados con ggplot.
  • Etiquetas en los scatter plot.

Importar la tabla de worldometer

Esta parte es probable que ya la haya puesto en otra ocasión, scraping sobre la web que tenemos todos en favoritos https://www.worldometers.info/coronavirus/

url ='https://www.worldometers.info/coronavirus/'

worldometers <- url %>%
  html()%>%
  html_nodes(xpath='//*[@id="main_table_countries_today"]') %>%
  html_table()
worldometers <- worldometers[[1]]

worldometers <- worldometers %>% rename(pais = `Country,Other` , muertes_habitante = `Deaths/1M pop`) %>%
  mutate(muertes = as.numeric(removePunctuation(TotalDeaths))) %>%
  select(pais,muertes_habitante, muertes)

Creamos un data frame con las variables necesarias para realizar un ranking de porcentaje de muertes sobre millón de habitantes.

Problemas con la librería OECD

En este punto más que mostrar como leer datos de la OECD casi pido ayuda para entender la API y la librería. Si ejecutamos:

library(OECD)
salud <- search_dataset('health', data = get_datasets(), ignore.case = TRUE)

No encuentro la tabla SH.XPD.CHEX.GD.ZS que contiene Current health expenditure (% of GDP), el gasto en sanidad sobre % de PIB; tampoco he llegado a comprender bien la API para bajarme el XML o el Json, a los 45 minutos de lucha me rendí. Así que opte por descargar el Excel, tratarlo manualmente y leerlo.

Importar Excel con R

En este caso hay muchas entradas y muchas posibilidades pero me gustaría expresar una opinión, el mejor paquete para leer Excel es readxl, nada de Java por medio y la sintaxis no puede ser más sencilla:

ub = 'C:\\Users\\ API_SH.XPD.CHEX.GD.ZS_DS2_en_excel_v3.xls'
pib_salud <- read_excel(ub, sheet = 'Data')
names(pib_salud) = c('pais', 'pib_salud')
pib_salud <- pib_salud %>% mutate(pais = case_when(
  pais == 'Iran, Islamic Rep.' ~ 'Iran',
  pais == 'United Kingdom' ~ 'UK',
  pais == 'United States' ~'USA' ,
  TRUE ~ pais))

Descargué los datos en formato Excel de la tabla con ID SH.XPD.CHEX.GD.ZS En el propio Excel quité lo que sobraba y así pude crear un data frame y adaptar los nombres con case_when que pueden darme problemas a la hora de cruzar con worldometer.

Not in con R

Un código al que creo que ya he hecho mención en otras ocasiones pero que periódicamente aparece en el blog.

`%notin%` = Negate(`%in%`)
pinta <- worldometers %>% filter(pais %notin% c('Total:','World')) %>% 
  arrange(desc(muertes)) %>% filter(row_number()<=20)

pinta <- left_join(pinta, pib_salud)

La tabla de worldometer tiene un total y un mundo que hay que eliminar, de ahí el uso del not in, por otro lado no vamos a pintar los casi 200 países que tenemos nos vamos a quedar con los 20 con mayor número de muertes y por ello ordenamos de forma descendente con arrange y en filter ponemos row_number <= 20. También unimos los datos de muertes por millón de habitantes con los datos de PIB.

Gráficos de ranking ordenados con ggplot

Esta cuestión es recurrente cuando se trabaja con ggplot, es necesario ordenar los datos a representar, por defecto ggplot ordena por la variable x y en este caso sería una ordenación lexicográfica y necesitamos una ordenación numérica:

pinta %>%
  mutate(País = fct_reorder(pais, muertes_habitante)) %>%
  ggplot( aes(x=País, y=muertes_habitante)) +
  geom_bar(stat="identity", fill="#908785", alpha=.5, width=.4) +
  coord_flip() +
  ylab("Muertes COVID por 1m. habitantes") +
  theme_classic()

Si seguís el blog siempre suelo crear el objeto con los datos a pintar, no suelo emplear ggplot para hacer ni agregación ni filtrado, esto siempre vendrá en un objeto previo que en este caso se llama pinta y sobre él creamos la variable País (sin problemas con las tildes) y este País con fct_reorder le damos el orden que requerimos, en este caso ordenamos por muertes por millón de habitantes y ya realizamos el gráfico que deseamos:

De igual manera podemos realizar el gráfico de % de gasto en salud sobre PIB:

 pinta %>%
  mutate(País = fct_reorder(pais, pib_salud)) %>%
  ggplot( aes(x=País, y=pib_salud)) +
  geom_bar(stat="identity", fill="#908785", alpha=.5, width=.4) +
  coord_flip() +
  ylab("% gasto en salud sobre PIB") +
  theme_classic()

Se puede opinar que el % de gasto en salud sobre el PIB no es un dato apropiado porque los grupos de poder en USA distorsionan este dato (y arruinan familias), pero no quiero dejar en mal lugar a mi país y hay algún dato que deja en muy mal lugar al actual Gobierno. Ya hablaré.

Etiquetas en los scatter plot

Que grandes gráficos de bivariables se hacen con ggplot. Además la sintaxis no puede ser más sencilla:

 base <- ggplot(pinta, aes(x = pib_salud, y = muertes_habitante))
base + geom_point()

Sobre una base vamos añadiendo de lo más sencillo a lo más complicado. Añadimos las etiquetas de la forma más sencilla creando un vector que las contenga y empleando geom_text

 paises <- pinta$pais
base + geom_point() +
  geom_text(aes(label = paises), size = 3, hjust = 0, nudge_x = 0.2)

Sin embargo, en ocasiones tenemos que complicar mucho la sintaxis para que esas etiquetas queden “elegantes” pero la librería ggrepel viene en nuestra ayuda:

#install.packages("ggrepel")
library(ggrepel)

base + geom_point() +
  geom_text_repel(aes(label = paises), size = 3)

base + geom_point() +
  geom_label_repel(aes(label = paises), size = 3) +
  ylab("Muertes por millón de habitantes") +
  xlab("gasto en salud sobre % de PIB") +
  theme_classic()

Que gran colocación de las etiquetas en el gráfico bivariable, podríamos añadirle más países y ponerle las cosas más difíciles.

Creo que es una serie de trucos muy útiles.

Leer archivos Excel con Python

Entrada sobre la importación de Excel con Python, un aporte que sirve para mi documentación y que es posible que sea de ayuda para muchos que se estén iniciando en el uso de Python y Pandas, aunque en este caso para la lectura del Excel usaremos tanto Pandas como la librería xlrd.

Lectura de Excel con Pandas

Lo más sencillo para importarnos en Python un Excel y crearnos un data frame de Pandas es:

import pandas as pd
archivo = 'C:/Users/Documents/ejemplo.xlsx'

df = pd.read_excel(archivo, sheetname='Hoja1')

df.describe()

La función read_excel será suficiente en el 80% de las ocasiones que realicemos esta tarea. Como es habitual en la ayuda tenéis perfectamente descritas sus posibilidades.

Lectura de Excel con xlrd

Es posible que necesitemos realizar tareas más complejas a la hora de leer archivos Excel y podemos usar xlrd. Vemos algunas de las posibilidades:

import xlrd 
 
archivo = 'C:/Users/rvaquerizo/Documents/ejemplo.xlsx'
  
wb = xlrd.open_workbook(archivo) 

hoja = wb.sheet_by_index(0) 
print(hoja.nrows) 
print(hoja.ncols) 
print(hoja.cell_value(0, 0))

open_workbook nos abre el Excel para trabajar con él. Seleccionamos hojas por índice (empezando por el 0) y con la hoja seleccionada podemos ver el número de filas (nrows) o columnas (ncols). Seleccionar una celda lo hacemos con cell_value mediante índices (empezando por el 0). Otras posibilidades:

archivo = 'C:/Users/rvaquerizo/Documents/ejemplo.xlsx'
  
wb = xlrd.open_workbook(archivo) 

hoja = wb.sheet_by_name('Hoja1') 
for i in range(0,hoja.nrows):
    print(hoja.cell_value(i,1))

Si por ejemplo deseamos saber las cabeceras, los nombres de las columnas:

archivo = 'C:/Users/rvaquerizo/Documents/ejemplo.xlsx'
  
wb = xlrd.open_workbook(archivo) 

hoja = wb.sheet_by_index(0) 
nombres = hoja.row(0)  
print(nombres)

Y mediante xlrd podemos crear data frames de pandas con lo que es posible realizar lecturas de rangos:

archivo = 'C:/Users/rvaquerizo/Documents/ejemplo.xlsx'
  
wb = xlrd.open_workbook(archivo) 

hoja = wb.sheet_by_index(0) 

# Creamos listas
filas = []
for fila in range(1,hoja.nrows):
    columnas = []
    for columna in range(0,2):
        columnas.append(hoja.cell_value(fila,columna))
    filas.append(columnas)

import pandas as pd
df = pd.DataFrame(filas)
df.head()

Hay alguna librería que lo hace de forma más elegante pero la importación de rangos de Excel con xlrd, una vez te familiarizas, me parece bastante sencilla. Espero que sea de utilidad

Datos agrupados en R con dplyr

Entrada rápida para ilustrar como crear un campo autonumérico por un factor, es una duda que me plantean, tienen datos de clientes y fechas y necesitan crear un autonumérico en R que les diga el número de orden de los eventos de una fecha. Algo parecido a lo que hacemos con el retain de R. Vamos a ilustrar la tarea con un ejemplo:

clientes = 100
id_cliente = rpois(clientes,10)
fecha = rpois(clientes, today()-rpois(clientes,5) )

eventos <- cbind.data.frame(id_cliente,fecha)

eventos$fecha <- as.Date(eventos$fecha, origin="1970-01-01")
eventos <- eventos %>% arrange(id_cliente,fecha) 

100 clientes que aparecen una o n veces con fechas asociadas, el primer paso que sugiero hacer es eliminar duplicados con dplyr:

eventos <- eventos %>% group_by(id_cliente, fecha) %>% 
  filter(row_number()==1) %>% as_tibble()

Agrupamos por cliente y fecha y nos quedamos con el primer registro, en este caso da igual quedarse con el primero que con el último. Ahora que no tenemos duplicados la agrupación ya no es por cliente y por fecha, como vamos a crear un valor agrupado por cliente haremos el group_by solo por cliente:

eventos <- eventos %>% group_by(id_cliente) %>% 
  mutate(evento=row_number()) %>% as_tibble()

Cada evento irá numerado de 1 a n gracias a row_number(), el mismo se reinicia a 1 cada vez que cambia el valor del group_by.

Los pilares de mi simulación de la extensión del COVID19

No debería publicar esta simulación de la extensión del CODVID10 o coronavirus porque puede disparar alarmas, provocar insultos, levantar ampollas,… el caso es que yo llevo 7 días de aislamiento más que el resto de España porque sólo había que ver los datos de Italia para saber lo que iba a pasar y no avisé a nadie para no disparar alarmas, provocar insultos, levantar ampollas… Y AL FINAL YO TENÍA RAZÓN. Así que os voy a exponer el motivo por el cual estoy muy asustado, bueno, hoy quiero mostraros el inicio de una simulación mala y sin fundamento que estoy realizando sobre la extensión en España del COVID19. Para hacerla vamos a emplear la siguiente información:

Y allá voy a comentaros que estoy montando. Se trata de poner a los 47 millones de españoles en una tabla, situarlos en unas coordenadas y dadas 5 personas iniciales ver como se propaga el virus municipio a municipio y, en 98 días, determinar cuantas personas pueden estar contagiadas, cuantas enfermas, cuantas sanas o cuantas desgraciadamente muertas. Esto no es que tenga lagunas, es que está inventado, pero no os creáis que las cifras oficiales son más fiables. Evidentemente, lo voy a hacer con R y dplyr. No lo subo a git porque el equipo que uso tiene un usuario de github que no es el adecuado, pero ya sabéis que el código está a vuestra disposición.

Creación de la tabla de personas edad

library(tidyverse)
library(pxR)
library(sqldf)

#detach("package:dplyr", unload=TRUE)

censo = "C:\\temp\\personales\\covid19\\0001.px"
datos <-  read.px(censo)
datos <- datos$DATA[[1]]
names(datos) = c("rango_edad", 'seccion', 'sexo', 'habitantes' )
datos <- data.frame(lapply(datos, as.character), stringsAsFactors=FALSE)

muestra <- datos %>% mutate(habitantes=round(as.numeric(habitantes)/10,0)) %>% 
  filter(seccion != "TOTAL" & sexo == "Ambos Sexos" & rango_edad != "Total") %>% select(-sexo) %>% 
  mutate(rango_edad = case_when( 
    rango_edad %in% c('0-4', '5-9', '10-14', '15-19', '20-24') ~ '<25',
    rango_edad %in% c('80-84', '85-89', '90-94', '95-99', '100 y más') ~ '80 >',
    TRUE ~rango_edad  )) 

muestra <- muestra %>% group_by(seccion,rango_edad) %>% summarise(habitantes=sum(habitantes))


espania <- muestra %>% group_by(seccion,rango_edad) %>% expand(count = seq(1:habitantes)) %>% as_tibble()

Nota: si no funciona la creación de la muestra hacéis detach de dplyr

Leemos los datos del censo que nos hemos descargado del INE, es un fichero px pero con el paquete pxR podemos manejarlo. Los datos que tenemos están a nivel de sección censal, rangos de edad, sexo y disponemos del número de habitantes. Con esta tabla de frecuencias generamos con expand una tabla que repite un registro tantas veces como digamos en una variable, es decir, repetirá la edad, el sexo, la sección censal tantas veces como habitantes tenga. Manejo una muestra del 10% porque el tema tiene un tiempo importante de procesamiento. Con esto también hago un llamamiento por si Amazon, Microsoft o Google pueden poner buenas máquinas en manos de los gestores de información (mal llamados científicos de datos ahora) de forma altruista. En fin, tenemos a todos los españoles, ahora vamos a ubicarlos con la cartografía por sección censal del INE.

Creación de la tabla de centroides municipal

library(maptools)
library(sf)
ub_shp = "C:\\temp\\mapas\\Seccion_censal\\SECC_CPV_E_20111101_01_R_INE.shp"
seccion_censal = readShapeSpatial(ub_shp)

mapa <- map_data(seccion_censal)

centroides <- mapa %>% group_by(OBJECTID = as.numeric(region)) %>% 
  summarise(centro_long=mean(long), centro_lat=mean(lat))

ggplot(data = centroides, aes(x = centro_long, y = centro_lat, group = 1)) +
  geom_polygon() 

secciones <- seccion_censal@data %>% mutate(seccion=as.character(CUSEC), municipio=as.character(CUMUN)) %>%
  select(OBJECTID,seccion,municipio)

municipios <- left_join(secciones,centroides) %>% group_by(municipio) %>% 
  summarise(centro_long=mean(long), centro_lat=mean(lat)) %>%
  select(municipio, centro_long, centro_lat)  

#Matriz de distancias
distancias <- sqldf(" select a.municipio, b.municipio as municipio2, 
                    a.centro_long, a.centro_lat, b.centro_long as centro_long2, b.centro_lat as centro_lat2
                    from municipios a , municipios b where a.municipio != b.municipio")

distancias <- distancias %>% mutate(distancia=sqrt((centro_long - centro_long2)**2 + (centro_lat-centro_lat2)**2))

Os habéis desgarcado el shapefile con las secciones censales de España y con ella calculamos el centroide de cada municipio, también he calculado una matriz de distancias porque, como veréis más adelante, la distancia de desplazamiento puede ser interesante para determinar como se mueve y como se expande el virus. En este punto está mi otra de mis reclamaciones, las compañías de telefonía podían ofrecer datos de movilidad para ayudarnos y controlar el movimiento de personas.

En fin, si cruzáis ambas tablas empieza la simulación (de mierda):

#Proceso
indices <- sample( 1:nrow( espania ), nrow(espania)/2)
espania2 <- espania[indices,]
espania2 <- espania2 %>% left_join(secciones) %>% 
  mutate(id_persona=row_number(),
         dia=0,contagiado=0, evolucion_dias=0)

sanos = espania2
contagiados = espania2[0,]
enfermos = espania2[0,]
curados = espania2[0,]
muertos = espania2[0,]

#Dia 1
#5 contagiados
dia <- sample_n(filter(espania2,seccion %in% c('2807908161','0810205003')) , 5)
contagiados <- inner_join(dia, select(sanos,id_persona)) %>% 
  mutate(contagiado=1)
lista_contagiados = unique(contagiados$id_persona)
sanos <- sanos %>% filter(id_persona %notin% lista_contagiados)

max_distancia =max(distancias$distancia,na.rm = T)

Tenemos una tabla con la población española por edad y ubicación, son 5 personas al azar de Igualada y Madrid las que empiezan todo… Veré si me atrevo a seguir contando porque lo que sigue me lo he inventado completamente.

Transformar todos los factores a carácter de mi data frame de R

En muchas ocasiones no quiero factores en mi dataframe cuando trabajo con R. Y estoy en mi derecho de poner una entrada sobre una de las tareas que más realizo y que siempre se me olvida el como la realizo, tardo menos en buscarlo en www.analisisydecision.es que entre mis programas:

df<- data.frame(lapply(df, as.character), stringsAsFactors=FALSE)

Todos los elementos factor ahora son character.

 

Mapa del COVID-19 por Comunidades Autónomas con R (más #rstats)

Estoy muy activo en twitter con el #covid-19 estos días y eso está dando lugar a algunas entradas en el blog. Sin embargo, he parado esa actividad porque el número de casos no me parece el indicador adecuado para medir la verdadera incidencia de la pandemia. Empiezo a tener posibles casos entre personas conocidas y no se realiza ningún test, permanecen en casa y son casos no informados. Sin embargo, quería que esta entrada sirviera de homenaje a la gente de Datadista que está recogiendo datos y realizan un seguimiento del número de camas ocupadas, uno de los mejores indicadores. Además sigo mi labor formativa con Rstats, hoy toca:

  • Mapa rápido y guarro de España con GADM
  • Homogeneización de textos con dplyr y tm
  • Complicar el web scraping con rvest

Esta entrada surge aquí:

Datadista pone a nuestra disposición datos actualizados por Comunidad Autónoma y con ellos podemos construir los mapas.

Mapa por Comunidad Autónoma con datos de Datadista

#Situación por Comunidad Autónoma
library(gganimate)
library(maptools)
library(raster)
library(maps)
library(tidyverse)

datadista = "https://raw.githubusercontent.com/datadista/datasets/master/COVID%2019/ccaa_covid19_casos.csv"

tabla_ccaa <- read.csv2(datadista, sep=',',encoding = 'UTF-8', check.names=FALSE)

Espania <- getData('GADM', country='Spain', level=1)
Espania$name = Espania$NAME_1
ccaa <- map_data(Espania)

pinta <- tabla_ccaa[,c(2,length(tabla_ccaa))]
names(pinta)=c("region","casos")

unique(ccaa$region)
unique(pinta$region)

ccaa <- ccaa %>% mutate(region=case_when(
  region == "Región de Murcia" ~ "Murcia",
  region == "Principado de Asturias" ~ "Asturias",
  region == "Comunidad de Madrid" ~ "Madrid",
  region == "Comunidad Foral de Navarra" ~ "Navarra",
  region == "Comunidad Valenciana" ~ "C. Valenciana",
  region == "Islas Canarias" ~ "Canarias",
  region == "Islas Baleares" ~ "Baleares",
  TRUE ~ region))

ccaa <- left_join(ccaa,pinta)

ggplot(data = ccaa, aes(x = long, y = lat, group = group)) +
  geom_polygon(aes(fill = casos)) +
  scale_fill_continuous(low="white",high="red") +
  labs(title = "Mapa del COVID-19 por Comunidad Autónoma") + 
  theme(panel.background =   
          element_rect(fill='#838596',colour='#838596'),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) + 
  theme(axis.line=element_blank(),axis.text.x=element_blank(),
        axis.text.y=element_blank(),axis.ticks=element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank())

Este código da lugar al mapa con el que se incia esta entrada. Como aspectos interesantes tiene descargar directamente el mapa con R de gadm o la lectura de cabeceras con formato fecha, algo que no conocía, nunca había usado check.names=FALSE. Por lo demás no es un código especialmente complicado. Pero me gustaría escribir sobre la relativización de los datos, no podemos decir que Madrid tiene 5 veces más casos que otra provincia si Madrid tiene 5 veces más habitantes que otra provincia, es necesario relativizar el número de casos y en este caso vamos a emplear el número de habitantes y además nos va a servir para hacer web scraping sobre una tabla de una página web.

Scraping sobre datosmacro. Mapa de casos por número de habitantes

El código empieza del siguiente modo:

library(rvest)
library(xml2)
library(tm)
numerea <- function(x) {as.numeric(sub(",",".",x)) }

url = 'https://datosmacro.expansion.com/demografia/poblacion/espana-comunidades-autonomas'

Si vais a la url indicada tenemos que extraer la tabla específica con el número de habitantes y para eso necesitamos saber en que lugar del código HTML se encuentra. En mi caso empleo Google Chrome, imagino que será análogo con otros navegadores. Hacemos lo siguiente:

Nos ubicamos sobre la tabla que deseamos scrapear (verbo regular de la primera conjugación) damos a inspeccionar y nos aparece la codificación, dentro de la codificación si pulsamos se marcará la tabla y Coppy + Copy XPath y con ello ya podemos crear un data frame con la tabla HTML:

poblacion <- url %>%
  html() %>%
  html_nodes(xpath='//*[@id="tb1"]') %>%
  html_table()
poblacion <- poblacion[[1]]

poblacion <- poblacion [,-4] %>% mutate(CCAA = removePunctuation(CCAA),
                                        CCAA = substr(CCAA,1,nchar(CCAA)-1),
                                        habitantes=numerea(removePunctuation(Población))) %>%
  rename(region=CCAA) %>%
  select(region, habitantes) %>% mutate(region=case_when(
    region == "Comunidad Valenciana" ~ "C. Valenciana",
    region == "Castilla La Mancha" ~ "Castilla-La Mancha",
    region == "Islas Baleares" ~ "Baleares",
    TRUE ~ region
  ))

en html_nodes hemos puesto el XPath y ya sabe que parte tiene que leer, como se genera una lista nos quedamos con el primer elemento de la lista y posteriormente se realiza la homogeneización de los nombres de las comunidades, eliminación de signos de puntuación con removePunctuation (que ha cambiado mi vida porque odio regex). Esta tabla puede ser cruzada con los datos de Datadista y crear un número de casos entre habitantes x 1000:

unique(poblacion$region)
unique(ccaa$region)

ccaa <- left_join(ccaa,poblacion)
ccaa$tasa_COVID <- (ccaa$casos/ccaa$habitantes)*1000


ggplot(data = ccaa, aes(x = long, y = lat, group = group)) +
  geom_polygon(aes(fill = tasa_COVID)) +
  scale_fill_continuous(low="white",high="red") +
  labs(title = "Mapa del COVID-19 por Comunidad Autónoma") + 
  theme(panel.background =   
          element_rect(fill='#838596',colour='#838596'),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) + 
  theme(axis.line=element_blank(),axis.text.x=element_blank(),
        axis.text.y=element_blank(),axis.ticks=element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank())

Y el resultado sigue siendo alarmante en Madrid pero la tonalidad del rojo cambia mucho en otras zonas de España, la importancia relativizar un dato.

Evolución del número de casos de coronavirus

Seguimos a vueltas con la (ya) pandemia y R y hoy quería traeros unos buenos ejemplos de uso de la librería dplyr para preparar datos. Se trata de ver una evolución del número de casos diarios para saber en qué punto tanto España como Italia pueden frenar el crecimiento de los casos de coronavirus, se trata de crear este gráfico:

Se observa como países como China o Korea vivieron un fuerte crecimiento que ahora se ha transformado en una caída del número de casos de coronavirus, pero parece que Irán ha estabilizado en 1000 casos diarios pero Italia y España siguen en fase de crecimiento por lo que no se espera que el comportamiento sea similar a China o Korea y es probable que el número de casos siga aumentando.

Para crear este gráfico estoy mejorando códigos que ya he venido utilizando y creo que son un buen ejemplo de uso de dplyr. El primer paso es crear el conjunto de datos inicial, código ya conocido:

library(lubridate)
library(ggplot2)
library(dplyr)
library(gridExtra)
datos <- read.csv2("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Confirmed.csv",
                   sep=',')


fechas <- seq(as.Date("2020/01/22"), as.Date(today()-1), "days")
fechas <- substr(as.character.Date(fechas),6,10)
names(datos) <- c("Provincia", "Pais","Latitud", "Longitud", fechas)

Para ver los datos de China es necesaria una agregación previa ya que China tiene los datos a nivel de región y no de país:

China <- filter(datos,Pais=="China")
China <- China %>% select(fechas) 
China <- data.frame(casos=as.numeric(mapply(sum,China)))
China <- China %>% mutate(fecha=row.names(China), dias=row_number())
China2 <- China %>% mutate(dias=dias+1) %>% rename(casos_menos1=casos) %>% select(-fecha)

China <- left_join(China,China2) %>% 
  mutate(incremento = casos-casos_menos1,
         incremento = case_when(is.na(incremento) ~0,
                                TRUE ~ incremento)) %>% 
  filter(dias<=50)

Aquí es interesante el uso de mapply para sumar todas las columnas del conjunto de datos. Por otro lado se trata de crear una variable incremento en función de los días que llevamos recogiendo datos, para ello lo que se hace es cruzar con los mismos datos pero sumamos un día al número de días de pandemia, de este modo, tenemos el dato del día y el dato del día anterior por lo que podemos crear una variable incremento. Para el resto de países, como la información no está a nivel de región, hacemos una función.

select_pais <- function(pais, numdias) {
  P <- filter(datos,Pais==pais) %>% select(fechas)
  P <- P %>% mutate(fecha=row.names(P), dias=row_number())
  P <- data.frame(casos=as.numeric(t(P)))
  P <- P %>% mutate(fecha=row.names(P), dias=row_number())
  P2 <- P %>% mutate(dias=dias+1) %>% rename(casos_menos1=casos) %>% select(-fecha)
  P <- left_join(P,P2) %>% 
    mutate(incremento = casos-casos_menos1,
           incremento = case_when(is.na(incremento) ~0,
                                  TRUE ~ incremento)) %>% 
    filter(dias<=numdias)
  return(P)}

Korea <- select_pais('Korea, South', 50) 
Iran <- select_pais('Iran', 50) 
Italia <- select_pais('Italy', 50) 
Espania <- select_pais('Spain', 50) 

Ya tenemos los datos con la forma deseada. Ahora nos toca realizar un gráfico para cada país:

evolucion <- function(pais) {
  df <- pais
  df$incremento <- ifelse(df$casos<=4,0,df$incremento)
  ggplot(df, aes(x=dias, y=incremento)) +
  geom_bar(stat="identity") + 
  scale_y_continuous(limits = c(0, 3000)) +
  ggtitle(paste0('Incremento de casos en ',deparse(substitute(pais)) )) +
  theme(plot.title = element_text(hjust = 0.5)) +
  xlab("Días desde 22/01/2020")}

p1 = evolucion(China)
p2 = evolucion(Korea)
p3 = evolucion(Espania)
p4 = evolucion(Italia)
p5 = evolucion(Iran)

Una función que realiza un gráfico de barras con ggplo2 para cada país y tiene un uso interesante de la función deparse que junto con sustitute nos permite poner en la función el nombre del data frame y no los datos que contiene. Es una forma sencilla de obtener el nombre de un data frame en una función. Si leéis que me reitero en algunas frases no os asustéis, sirve para facilitar búsquedas.

Ahora empleamos una librería que cambió mi vida gridExtra y podemos realizar el gráfico que abre esta entrada:

grid.arrange(p1, p2,p3,p4,p5,ncol=1)

En este punto sólo queda ejecutar este código todos los días y esperar a que Italia y España lleguen a ese máximo. Saludos.

Seguimiento del coronavirus en España por Comunidad Autónoma. Extraer información de un PDF con R

Una entrada anterior del blog ha dado lugar a una conversación interesante en twitter:

Es necesario obtener los datos del Ministerio y estos datos se hayan en un pdf (https://www.mscbs.gob.es/profesionales/saludPublica/ccayes/alertasActual/nCov-China/documentos/Actualizacion_41_COVID-19.pdf) Bien, tendremos que leer el pdf y crear un data frame para poder trabajar con estos datos. Para leer el pdf vamos a emplear el paquete de R tabulizer y la función extract_table pero necesitamos “algo de talento”.

La función extract_tables nos permite extraer información de un archivo pdf con R pero es necesario especificar la página que deseamos leer (fácil) y el área que deseamos leer. ¿El área de la hoja del pdf que deseamos leer? ¿Eso qué es? Pues ese es el “talento” que necesitamos, con la función locale_areas vamos a encontrar ese área. Veamos el código necesario:

library(tidyverse)
library(tabulizer)

ministerio = "https://www.mscbs.gob.es/profesionales/saludPublica/ccayes/alertasActual/nCov-China/documentos/Actualizacion_41_COVID-19.pdf"

area <- locate_areas(ministerio, pages = 2)

Al ejecutar ese código podemos seleccionar a mano alzada el área que deseamos seleccionar de la página específica del pdf:

Ya estamos en disposición de ver el área a leer:

area[[1]]

Un poco complicado, pero una vez sabemos el área crear un objeto con R que contenga la información actualizada por Comunidad Autónoma de los datos del coronavirus en España con R es así de sencillo:

pdf_lista <- extract_tables(
  ministerio,
  output = "data.frame",
  pages = c(2),
  area = list(
    c(337.89431,  90.69972, 684.23597, 510.25186)
  ),
  guess = FALSE,
  encoding = "UTF-8"
)


datos <- data.frame(pdf_lista[1])

Ahora ya tenéis los datos por Comunidad Autónoma actualizados, solo queda que alguien los tabule por día y haga representaciones gráficas. Ahora a esperar que el Ministerio no cambie el pdf. Mañana haremos algo con ellos.

Estimación de la evolución de casos del coronavirus en España

Ayer escribrí sobre la obtención de los datos del coronavirus con R y después me disponía ha escribir sobre modelos de regresión no lineal, hacer una estimación del coronavirus en España,… Pero estuve hablando con una amiga residente en Italia y allí el número de casos está dos semanas por delante de España, bueno, dos semanas exactamente no, 10 días:

library(lubridate)
library(ggplot2)
library(dplyr)
library(reshape2)
datos <- read.csv2("https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Confirmed.csv",
                   sep=',')

fechas <- seq(as.Date("2020/01/22"), as.Date(today()-1), "days")
fechas <- substr(as.character.Date(fechas),6,10)
names(datos) <- c("Provincia", "Pais","Latitud", "Longitud", fechas)

esp_ita <- data.frame(fecha=fechas)
esp_ita <- cbind.data.frame(esp_ita, Espania = t(datos %>% filter(Pais=="Spain") %>% select(fechas)))
esp_ita <- cbind.data.frame(esp_ita, Italia = t(datos %>% filter(Pais=="Italy") %>% select(fechas)))

p <- ggplot(esp_ita, aes(x=fecha)) + 
  geom_line(aes(y=Espania, group = 1, color="España")) + 
  geom_line(aes(y=Italia, group = 1, color="Italia")) + 
  scale_color_manual(values = c("España" = "red", "Italia" = "blue")) +
  xlab("") + ylab("")
p

A la vista de las dos evoluciones hace pensar que el número de casos en España siga la misma evolución que sigue en Italia, por eso en este caso en vez de manejar fechas vamos a manejar días de evolución y al dato de España vamos a quitarle 10 días:

esp_ita <- data.frame(fecha=fechas) 
esp_ita <- esp_ita %>% mutate(dia=row_number())
Italia <- data.frame(Italia=t(datos %>% filter(Pais=="Italy") %>% select(fechas)))
Italia <- Italia %>% mutate(dia = row_number())

Espania <- data.frame(Espania=t(datos %>% filter(Pais=="Spain") %>% select(fechas)))
Espania <- Espania %>% filter(row_number() >= 10) %>% mutate(dia=row_number())

esp_ita <- esp_ita %>% left_join(Espania) %>% left_join(Italia)

p <- ggplot(esp_ita, aes(x=fecha)) + 
  geom_line(aes(y=Espania, group = 1, color="España")) + 
  geom_line(aes(y=Italia, group = 1, color="Italia")) + 
  scale_color_manual(values = c("España" = "red", "Italia" = "blue")) +
  xlab("") + ylab("")
p

Mi estimación para los próximos 10 días sobre la evolución del coronavirus en España está clara, ahora solo queda que no salga el ejército a la calle como en Milán.