Contenido

Introducción

Este texto fue creado para complementar la serie de videos en el canal DaniMedi: https://youtube.com/playlist?list=PLiR4mMxzSHWhM8Ek5Dy9qvSFGbcmSCWgv

Es importante poder visualizar los datos epidemiológicos obtenidos para poder crear hipótesis o para realizar comparaciones de forma más sencilla. El objetivo de este texto y de la serie de videos es brindar herramientas básicas para generar gráficos epidemiológicos simples usando el programa estadístico R (con RStudio).

A la fecha de crear este contenido estamos pasando por la pandemia COVID-19 y los datos recogidos de esta pandemia suelen almacenarse en forma de bases de datos mantenidas por organizaciones, instituciones o países; en esta situación es de especial importancia poder saber cómo trabajar con bases de datos si se está haciendo algún estudio relacionado al tema. Se usará la información epidemiológica sobre COVID-19 que se encuentra en algunas bases de datos para explicar los temas a tratar, espero que la información sea de utilidad.


Dedicatoria:

Es fácil concentrarse en las gráficas y la exploración de datos y en un momento olvidar que cada punto en el gráfico representa una persona, ya sea una persona enferma o quizá una persona que ya no está con nosotros. Esta serie de videos y este corto texto están dedicados a todas aquellas personas con las cuales se formaron estos gráficos.


Creando un proyecto en RStudio

Lo primero es tener el programa funcionando en la computadora, hice un video anteriormente sobre cómo instalar el programa: Descargando R, RStudio y Rtools.

En esta ocasión es importante que creemos un proyecto que nos permitirá tener nuestra información más ordenada en la carpeta que elijamos, para esto usaremos la opción de RStudio de crear un nuevo proyecto.

Al crear un nuevo proyecto creamos un ambiente de trabajo en la carpeta que elejimos en el que podemos tener acceso a los archivos que estén en la carpeta, podemos tener un historial con los cambios que hagamos y otras opciones más.

Obtención de información (bases de datos)

Hay diversas bases de datos de diferentes temas, incluso existen revistas en las que se publican bases de datos, como Scientific Data. En nuestro caso usaremos información sobre COVID-19 como se mencionó, y para esto usaremos dos fuentes, para mostrar bases de datos con diferentes presentaciones.

Usaremos la base de datos usada por Our World in Data y usaremos las bases de datos de la Plataforma Nacional de Datos Abiertos del Ministerio de Salud de Perú.

Podemos encontrar estas bases de datos en internet y para poder descargarlas directamente con el programa necesitamos el link con el cuál se descarga o se accede a la base de datos, por lo cual tenemos que buscar dónde tienen disponibles estas bases de datos. Para ejemplificar buscaremos los links de la base datos de Our World in Data y de la base de datos de los casos confirmados de COVID-19 en los Datos Abiertos del MINSA.

Para descargar estos archivos usando el programa se usará la función download.file(), pero no es importante que entiendan su funcionamiento, pueden simplemente reemplazar los links y elegir el nombre del archivo que quieran (no se olviden de colocar la extensión en el nombre del archivo, generalmente es .csv). Además, es de utilidad que al menos guarden la fecha de descarga, para tener evidencia de cuándo es que descargaron los datos, ya que suelen actualizarse todos los días, para esto usaremos la función Sys.time()

# datos abiertos sobre los positivos (Perú)
url <- "https://cloud.minsa.gob.pe/s/Y8w3wHsEdYQSZRp/download"
filename <- "peru_positivos_covid.csv"
download.file(url, destfile = filename)
# datos de Our World in Data (mundo)
url <- "https://raw.githubusercontent.com/owid/covid-19-data/master/public/data/owid-covid-data.csv"
filename <- "owid-covid-data.csv"
download.file(url, destfile = filename)
Sys.time()
## [1] "2021-01-18 15:23:01 -05"

Es sencillo ingresar estas bases de datos descargadas al programa, si están en el formato CSV (que es el caso) usamos la función read.csv() y podemos guardar esta lectura en una variable que guardará esta información (esto de guardar en una variable se explicará en el siguiente tema).

peru_positivos <- read.csv("peru_positivos_covid.csv")
owid_data <- read.csv("owid-covid-data.csv")

Al descargarse deben aparecer estos archivos en la carpeta del proyecto y para revisar que estén descargados de forma correcta pueden usar la función de RStudio View().

View(peru_positivos)
View(owid_data)

Ahora examinemos las bases de datos y comparemos. Podemos notar que en la base de datos de Our World in Data la base de datos cuenta con mucha información epidemiológica e incluso demográfica ordenada por localización y por fecha. Por otro lado, la base de datos de Perú tiene información mucho más limitada y solo presenta información de los casos positivos (porque tiene una base de datos similar separada para los fallecidos).

Más adelante usaremos estas dos bases de datos para explicar cómo podemos analizar la información que se encuentra en ellas.

Preparación de los datos y funciones básicas necesarias

Para poder comenzar a trabajar con los datos necesitamos conocer aspectos básicos del programa. Para esto, es más fácil aprender practicando, así que es recomendable que primero se familiaricen haciendo cosas simples para después entender lo que se muestre después.

Asignación y operaciones matemáticas

Pueden realizarse operaciones matemáticas básicas con el programa, además, podemos dar valores a nuestras variables usando = o <-.

1 + 2 * 3 - 4
## [1] 3
x <- 2
y = 3
x + y - 1
## [1] 4

Vectores

Podemos crear listas de valores (no necesariamente numéricos) usando c() separando los elementos por comas. Podemos aplicar funciones a estas listas para obtener algunos valores como la suma (usando sum()) o la media aritmética (usando mean()).

x <- c(1, 13, 6, 45, 7)
x + 1
## [1]  2 14  7 46  8
sum(x)
## [1] 72
mean(x)
## [1] 14.4

Data sets

Podemos leer archivos en el formato CSV usando la función read.csv y tenemos que guardar esto en una variable.

# esto ya se vio en el tema anterior
owid_data <- read.csv("owid-covid-data.csv")

Para visualizar un data set podemos usar la función View().

View(owid_data)

Entender esto de los vectores es importante porque las columnas o las filas de un data set se comportan como vectores, por lo que pueden realizarse fácilmente las mismas operaciones.

Filtrando y seleccionando

Es importante poder filtrar y seleccionar la información que tengamos en nuestros data sets, se entenderá mejor con ejemplos.

Creamos este vector,

x <- c(10, 2, 22, 14, 5, 14)

podemos seleccionar el número que está en la posición 4,

x[4]
## [1] 14

podemos seleccionar los números menores de 20,

x[x < 20]
## [1] 10  2 14  5 14

podemos seleccionar los números iguales a 14, esto lo hacemos usando el símbolo == que representa igualdad (recordemos que el símbolo = se usa para asignar).

x[x == 14]
## [1] 14 14

Los vectores tienen una sola dimensión, mientras que los data sets tienen 2 dimensiones, entonces para los data sets tenemos que usar una coma ‘,’ en el interior de los corchetes [] que separe las filas de las columnas.

Creamos una pequeña matriz (solo es por el ejemplo, este código no tiene importancia),

x <- matrix(c(10,2,22,14,5,13),ncol=2)
x
##      [,1] [,2]
## [1,]   10   14
## [2,]    2    5
## [3,]   22   13

ahora podemos seleccionar un elemento por sus coordenadas,

x[3,2]
## [1] 13

podemos seleccionar filas enteras si dejamos la parte de la columna vacía (pero con la coma),

x[1, ]
## [1] 10 14

o podemos seleccionar filas enteras si dejamos la parte de la fila vacía (pero con la coma).

x[ ,2]
## [1] 14  5 13

En el caso de nuestro data set podemos usar otro “truco” que es el uso del símbolo $ con el que podemos seleccionar las columnas por su nombre:

owid_data$location
owid_data$new_deaths_per_million

Y podemos usar las funciones de vectores como ya se mencionó.

# na.rm = TRUE es para remover los valores faltantes
sum(owid_data$total_cases, na.rm = TRUE)
## [1] 2331052900
mean(owid_data$new_deaths, na.rm = TRUE)
## [1] 41.51676

Combinando todo lo que ahora sabemos podemos obtener la información que solo corresponde a Perú (filas en las que location contenga la palabra Peru). Notar que la palabra Peru está entre comillas, esto es importante porque las comillas le indican al programa que la palabra no se trata de una variable y tampoco olvidar la coma.

peru_data <- owid_data[owid_data$location == "Peru", ]

NOTA: la forma mostrada para filtrar la información es la más intuitiva, pero quizá una forma que es más preferible para muchos es usando la función subset() que hace precisamente lo mismo, seleccionar filas en base a una condición, funciona del siguiente modo:

peru_data <- subset(owid_data, location == "Peru")

En el resto del texto se muestra la forma de los corchetes y en los videos también se usa esta forma, pero si prefieren esta forma pueden usarla, funciona del mismo modo.

Comprobemos que obtuvimos lo que queríamos.

View(peru_data)

Es importante entender los conceptos hasta ahora que lo usaremos constantemente para obtener mucha de la información que queremos, por ejemplo, podemos obtener las muertes totales a la fecha en el Perú (de acuerdo a la base de datos):

# na.rm es para eliminar los valores en blanco
sum(peru_data$new_deaths, na.rm = TRUE)
## [1] 25856
max(peru_data$total_deaths, na.rm = TRUE)
## [1] 25856

Podemos crear pequeños programas para hacer esto con otros países (no olvidar las comillas en los nombres):

country <- "Bolivia"
country_data <- owid_data[owid_data$location == country, ]
max(country_data$total_deaths, na.rm = TRUE)
## [1] 3884

Dejo en su creatividad las posibilidades de lo mostrado hasta ahora, con lo que se aprenderá en los temas siguientes podrán tener más posibilidades de mejorar sus habilidades en la exploración de datos.

Es importante que se sientan en cierto punto familiarizados con esto, así que pueden jugar con estas cosas hasta agarrar el hilo de cómo filtrar la información, porque lo usaremos en todo momento y pueden decidir si usar la forma con los corchetes [] o si prefieren usar la función subset().

P.D. no se olviden de la coma ;)

Frecuencias

En el caso de la epidemiología en una epidemia o pandemia es muy importante obtener frecuencias. En la base de datos de Our World in Data ya nos dan la frecuencia, por ejemplo, en la forma de muertes totales, muertes nuevas, casos totales y casos nuevos separados por fecha para cada localización, es decir, ya contamos con esos datos, pero este no siempre es el caso. Esta vez trabajaremos con la base de datos de casos positivos en la Plataforma Nacional de Datos Abiertos del Ministerio de Salud de Perú (la base de datos para los fallecidos es similar y se puede trabajar de manera similar también).

Recordemos cómo presenta esta información.

peru_positivos <- read.csv("peru_positivos_covid.csv")
View(peru_positivos)

En este caso solo nos muestra una lista de los eventos y no nos dice directamente el número de muertes diarias o totales, entonces es algo que tenemos que obtenerlo nosotros.

Hay una función muy útil para obtener la frecuencia y se llama table().

x <- c("mujer", "hombre", "mujer", "hombre","hombre")
table(x)
## x
## hombre  mujer 
##      3      2

Si lo pensamos, podemos aplicar esto de diversas maneras en nuestro data set, además del sexo, podemos obtener el número de eventos (casos positivos) por fecha de la siguiente manera:

table(peru_positivos$FECHA_RESULTADO)
## 
## 20200306 20200307 20200308 20200309 20200310 20200311 20200312 20200313 
##        1        5        2        3        1        8        8       10 
## 20200314 20200315 20200316 20200317 20200318 20200319 20200320 20200321 
##       19       28       20       27       56       62       56       30 
## 20200322 20200323 20200324 20200325 20200326 20200327 20200328 20200329 
##       39       33       64      100       52       34      136      142 
## 20200330 20200331 20200401 20200402 20200403 20200404 20200405 20200406 
##      130      250      117      222      292      833      444      647 
## 20200407 20200408 20200409 20200410 20200411 20200412 20200413 20200414 
##     1042     1083      817     1038     1404      738     1284     1041 
## 20200415 20200416 20200417 20200418 20200419 20200420 20200421 20200422 
##     1383     1329     1109     1407     1076     2039     2171     2104 
## 20200423 20200424 20200425 20200426 20200427 20200428 20200429 20200430 
##     2056     2239     2397     1422     3399     3367     4238     3372 
## 20200501 20200502 20200503 20200504 20200505 20200506 20200507 20200508 
##     2625     3369     1922     3990     3969     3634     3612     3297 
## 20200509 20200510 20200511 20200512 20200513 20200514 20200515 20200516 
##     2469      936     3601     4348     4441     3739     4304     4125 
## 20200517 20200518 20200519 20200520 20200521 20200522 20200523 20200524 
##     1785     5130     5198     5290     5514     6100     5399     1792 
## 20200525 20200526 20200527 20200528 20200529 20200530 20200531 20200601 
##     6968     6919     7371     6425     5745     4613     2262     4242 
## 20200602 20200603 20200604 20200605 20200606 20200607 20200608 20200609 
##     3774     2614     4029     4944     4764     2637     4743     5310 
## 20200610 20200611 20200612 20200613 20200614 20200615 20200616 20200617 
##     5726     5069     4661     4500     2441     4363     3376     3915 
## 20200618 20200619 20200620 20200621 20200622 20200623 20200624 20200625 
##     3436     3447     3526     1446     4335     3768     4109     4154 
## 20200626 20200627 20200628 20200629 20200630 20200701 20200702 20200703 
##     4331     3526     1598     2729     3748     3648     3349     3862 
## 20200704 20200705 20200706 20200707 20200708 20200709 20200710 20200711 
##     3518     2299     3783     4035     2598     2495     4913     4246 
## 20200712 20200713 20200714 20200715 20200716 20200717 20200718 20200719 
##     2380     3114     4194     4432     4535     5141     5066     2228 
## 20200720 20200721 20200722 20200723 20200724 20200725 20200726 20200727 
##     3756     4815     5972     5474     5960     5626     2950     7071 
## 20200728 20200729 20200730 20200731 20200801 20200802 20200803 20200804 
##     3017     6721     7248     7601     6697     3194     7818     7754 
## 20200805 20200806 20200807 20200808 20200809 20200810 20200811 20200812 
##     7508     8442     7262     6535     3329     7172     7590     5069

Vemos que ya tenemos el número de casos positivos nuevos por día en el Perú.

Ahora, para obtener los casos positivos totales simplemente tenemos que obtener una frecuencia acumulada, para esto podemos usar la función cumsum().

freqs <- table(peru_positivos$FECHA_RESULTADO)
cumsum(freqs)
## 20200306 20200307 20200308 20200309 20200310 20200311 20200312 20200313 
##        1        6        8       11       12       20       28       38 
## 20200314 20200315 20200316 20200317 20200318 20200319 20200320 20200321 
##       57       85      105      132      188      250      306      336 
## 20200322 20200323 20200324 20200325 20200326 20200327 20200328 20200329 
##      375      408      472      572      624      658      794      936 
## 20200330 20200331 20200401 20200402 20200403 20200404 20200405 20200406 
##     1066     1316     1433     1655     1947     2780     3224     3871 
## 20200407 20200408 20200409 20200410 20200411 20200412 20200413 20200414 
##     4913     5996     6813     7851     9255     9993    11277    12318 
## 20200415 20200416 20200417 20200418 20200419 20200420 20200421 20200422 
##    13701    15030    16139    17546    18622    20661    22832    24936 
## 20200423 20200424 20200425 20200426 20200427 20200428 20200429 20200430 
##    26992    29231    31628    33050    36449    39816    44054    47426 
## 20200501 20200502 20200503 20200504 20200505 20200506 20200507 20200508 
##    50051    53420    55342    59332    63301    66935    70547    73844 
## 20200509 20200510 20200511 20200512 20200513 20200514 20200515 20200516 
##    76313    77249    80850    85198    89639    93378    97682   101807 
## 20200517 20200518 20200519 20200520 20200521 20200522 20200523 20200524 
##   103592   108722   113920   119210   124724   130824   136223   138015 
## 20200525 20200526 20200527 20200528 20200529 20200530 20200531 20200601 
##   144983   151902   159273   165698   171443   176056   178318   182560 
## 20200602 20200603 20200604 20200605 20200606 20200607 20200608 20200609 
##   186334   188948   192977   197921   202685   205322   210065   215375 
## 20200610 20200611 20200612 20200613 20200614 20200615 20200616 20200617 
##   221101   226170   230831   235331   237772   242135   245511   249426 
## 20200618 20200619 20200620 20200621 20200622 20200623 20200624 20200625 
##   252862   256309   259835   261281   265616   269384   273493   277647 
## 20200626 20200627 20200628 20200629 20200630 20200701 20200702 20200703 
##   281978   285504   287102   289831   293579   297227   300576   304438 
## 20200704 20200705 20200706 20200707 20200708 20200709 20200710 20200711 
##   307956   310255   314038   318073   320671   323166   328079   332325 
## 20200712 20200713 20200714 20200715 20200716 20200717 20200718 20200719 
##   334705   337819   342013   346445   350980   356121   361187   363415 
## 20200720 20200721 20200722 20200723 20200724 20200725 20200726 20200727 
##   367171   371986   377958   383432   389392   395018   397968   405039 
## 20200728 20200729 20200730 20200731 20200801 20200802 20200803 20200804 
##   408056   414777   422025   429626   436323   439517   447335   455089 
## 20200805 20200806 20200807 20200808 20200809 20200810 20200811 20200812 
##   462597   471039   478301   484836   488165   495337   502927   507996

Y ya tenemos los casos positivos totales por fecha en el Perú.

Del mismo modo que en el tema anterior, podemos usar nuestra creatividad para crear programas que realicen acciones específicas, por ejemplo, podemos obtener los casos positivos totales por fecha para un departamento en específico:

dep <- "AREQUIPA"
sel <- peru_positivos[peru_positivos$DEPARTAMENTO == dep, ]
freqs <- table(sel$FECHA_RESULTADO)
cumsum(freqs)
## 20200307 20200311 20200319 20200320 20200322 20200323 20200324 20200325 
##        1        2        3        4        7        9       10       15 
## 20200329 20200330 20200331 20200402 20200404 20200405 20200406 20200407 
##       21       26       30       37       51       60       67       71 
## 20200408 20200409 20200410 20200411 20200412 20200413 20200414 20200415 
##       76       82       84       96      111      123      136      163 
## 20200416 20200417 20200418 20200419 20200420 20200421 20200422 20200423 
##      204      216      249      252      280      308      351      407 
## 20200424 20200425 20200426 20200427 20200428 20200429 20200430 20200501 
##      473      497      501      577      653      764      833      871 
## 20200502 20200503 20200504 20200505 20200506 20200507 20200508 20200509 
##      902      909      969     1026     1082     1129     1168     1227 
## 20200510 20200511 20200512 20200513 20200514 20200515 20200516 20200517 
##     1252     1303     1407     1519     1640     1736     1821     1829 
## 20200518 20200519 20200520 20200521 20200522 20200523 20200524 20200525 
##     1915     2037     2193     2337     2534     2680     2715     2907 
## 20200526 20200527 20200528 20200529 20200530 20200531 20200601 20200602 
##     3035     3221     3407     3579     3675     3695     3756     3855 
## 20200603 20200604 20200605 20200606 20200607 20200608 20200609 20200610 
##     3910     3961     4079     4217     4251     4355     4481     4625 
## 20200611 20200612 20200613 20200614 20200615 20200616 20200617 20200618 
##     4726     4845     5019     5038     5131     5221     5389     5457 
## 20200619 20200620 20200621 20200622 20200623 20200624 20200625 20200626 
##     5571     5645     5763     5877     6065     6296     6448     6614 
## 20200627 20200628 20200629 20200630 20200701 20200702 20200703 20200704 
##     6886     6944     7002     7136     7322     7437     7556     7766 
## 20200705 20200706 20200707 20200708 20200709 20200710 20200711 20200712 
##     7834     7912     8087     8255     8387     8720     9040     9094 
## 20200713 20200714 20200715 20200716 20200717 20200718 20200719 20200720 
##     9206     9408     9661     9869    10067    10390    10460    10619 
## 20200721 20200722 20200723 20200724 20200725 20200726 20200727 20200728 
##    11018    11541    11883    12278    12728    12829    13323    13471 
## 20200729 20200730 20200731 20200801 20200802 20200803 20200804 20200805 
##    13829    14361    14997    15650    15741    16411    17006    17830 
## 20200806 20200807 20200808 20200809 20200810 20200811 20200812 
##    18664    19323    19924    20031    20676    21337    21692

Podemos hacer lo mismo para los casos positivos nuevos por día y trabajar de forma similar con la base de datos de fallecidos o incluso separando por sexo, provincia, etc.

Aquí dejo una ayuda mostrando cómo obtener casos positivos nuevos por día en hombres y en mujeres por separado.

hombres <- peru_positivos[peru_positivos$SEXO == "MASCULINO", ]
table(hombres$FECHA_RESULTADO)
## 
## 20200306 20200307 20200308 20200309 20200310 20200311 20200312 20200313 
##        1        4        1        2        1        2        2        6 
## 20200314 20200315 20200316 20200317 20200318 20200319 20200320 20200321 
##       12       15       14       16       28       34       24       11 
## 20200322 20200323 20200324 20200325 20200326 20200327 20200328 20200329 
##       17       16       34       65       29       23       90       86 
## 20200330 20200331 20200401 20200402 20200403 20200404 20200405 20200406 
##       72      151       79      137      180      528      270      388 
## 20200407 20200408 20200409 20200410 20200411 20200412 20200413 20200414 
##      627      705      477      631      899      426      858      630 
## 20200415 20200416 20200417 20200418 20200419 20200420 20200421 20200422 
##      857      745      607      862      660     1323     1502     1279 
## 20200423 20200424 20200425 20200426 20200427 20200428 20200429 20200430 
##     1240     1372     1415      857     2112     2026     2707     1980 
## 20200501 20200502 20200503 20200504 20200505 20200506 20200507 20200508 
##     1502     1993     1102     2300     2299     2142     2157     1927 
## 20200509 20200510 20200511 20200512 20200513 20200514 20200515 20200516 
##     1473      522     2301     2643     2722     2210     2689     2580 
## 20200517 20200518 20200519 20200520 20200521 20200522 20200523 20200524 
##     1020     3271     3125     3227     3282     3690     3251     1001 
## 20200525 20200526 20200527 20200528 20200529 20200530 20200531 20200601 
##     4059     4044     4235     3764     3379     2677     1246     2355 
## 20200602 20200603 20200604 20200605 20200606 20200607 20200608 20200609 
##     2040     1492     2167     2657     2669     1464     2608     2960 
## 20200610 20200611 20200612 20200613 20200614 20200615 20200616 20200617 
##     3063     2667     2467     2348     1372     2432     1832     2172 
## 20200618 20200619 20200620 20200621 20200622 20200623 20200624 20200625 
##     1811     1815     1825      816     2303     2017     2281     2226 
## 20200626 20200627 20200628 20200629 20200630 20200701 20200702 20200703 
##     2304     1835      859     1484     2009     2020     1713     2066 
## 20200704 20200705 20200706 20200707 20200708 20200709 20200710 20200711 
##     1846     1217     2038     2173     1387     1277     2480     2230 
## 20200712 20200713 20200714 20200715 20200716 20200717 20200718 20200719 
##     1280     1798     2211     2362     2374     2732     2613     1153 
## 20200720 20200721 20200722 20200723 20200724 20200725 20200726 20200727 
##     2048     2698     3266     2810     3071     2758     1526     3529 
## 20200728 20200729 20200730 20200731 20200801 20200802 20200803 20200804 
##     1592     3492     3671     3776     3298     1611     3966     3908 
## 20200805 20200806 20200807 20200808 20200809 20200810 20200811 20200812 
##     3821     4247     3651     3367     1781     3487     3884     2625
mujeres <- peru_positivos[peru_positivos$SEXO == "FEMENINO", ]
table(mujeres$FECHA_RESULTADO)
## 
## 20200307 20200308 20200309 20200311 20200312 20200313 20200314 20200315 
##        1        1        1        6        6        4        7       13 
## 20200316 20200317 20200318 20200319 20200320 20200321 20200322 20200323 
##        6       11       28       28       32       19       22       17 
## 20200324 20200325 20200326 20200327 20200328 20200329 20200330 20200331 
##       30       35       23       11       46       56       58       99 
## 20200401 20200402 20200403 20200404 20200405 20200406 20200407 20200408 
##       38       85      112      305      174      259      415      378 
## 20200409 20200410 20200411 20200412 20200413 20200414 20200415 20200416 
##      340      407      505      312      426      411      526      584 
## 20200417 20200418 20200419 20200420 20200421 20200422 20200423 20200424 
##      502      545      416      716      669      825      816      867 
## 20200425 20200426 20200427 20200428 20200429 20200430 20200501 20200502 
##      982      565     1287     1341     1531     1392     1123     1376 
## 20200503 20200504 20200505 20200506 20200507 20200508 20200509 20200510 
##      820     1690     1670     1492     1455     1370      996      414 
## 20200511 20200512 20200513 20200514 20200515 20200516 20200517 20200518 
##     1300     1705     1719     1529     1615     1545      765     1859 
## 20200519 20200520 20200521 20200522 20200523 20200524 20200525 20200526 
##     2073     2063     2232     2410     2148      791     2909     2875 
## 20200527 20200528 20200529 20200530 20200531 20200601 20200602 20200603 
##     3136     2661     2366     1936     1016     1887     1734     1122 
## 20200604 20200605 20200606 20200607 20200608 20200609 20200610 20200611 
##     1862     2287     2095     1173     2135     2350     2663     2402 
## 20200612 20200613 20200614 20200615 20200616 20200617 20200618 20200619 
##     2194     2152     1069     1931     1544     1743     1625     1632 
## 20200620 20200621 20200622 20200623 20200624 20200625 20200626 20200627 
##     1701      630     2032     1751     1828     1928     2027     1691 
## 20200628 20200629 20200630 20200701 20200702 20200703 20200704 20200705 
##      739     1245     1739     1628     1636     1796     1672     1082 
## 20200706 20200707 20200708 20200709 20200710 20200711 20200712 20200713 
##     1745     1862     1211     1218     2433     2016     1100     1316 
## 20200714 20200715 20200716 20200717 20200718 20200719 20200720 20200721 
##     1983     2070     2161     2409     2453     1075     1708     2117 
## 20200722 20200723 20200724 20200725 20200726 20200727 20200728 20200729 
##     2706     2664     2889     2868     1424     3542     1425     3229 
## 20200730 20200731 20200801 20200802 20200803 20200804 20200805 20200806 
##     3577     3825     3399     1583     3852     3846     3687     4195 
## 20200807 20200808 20200809 20200810 20200811 20200812 
##     3611     3168     1548     3685     3706     2444

P.D. no se olviden de la coma x2 ;)

Transformación de datos: ajustando a la población

Es de utilidad también poder transformar las cifras que obtengamos, por ejemplo, para comparar valores en distintas poblaciones se suele ajustar las cifras a la población total.

En la base de datos de Our World in Data ya contamos con una columna de valores ajustados a la población del país, así que no habría necesidad de calcular algunas cosas. En el caso de las bases de datos de Datos Abiertos del MINSA no se tiene esta información así que hay que calcularla y es lo que mostraré a continuación.

Para ajustar a la población lo primero es definir la población y obtener la cifra. Usaré a Perú en este ejemplo.

# aprox
peru_pop <- 31.99 # millones

Ahora, simplemente es cuestión de dividir los casos positivos nuevos por día entre esta población (recordemos cómo lo hicimos en el tema anterior).

peru_positivos <- read.csv("peru_positivos_covid.csv")
freqs <- table(peru_positivos$FECHA_RESULTADO)
adj_freqs <- freqs / peru_pop
adj_freqs
## 
##     20200306     20200307     20200308     20200309     20200310     20200311 
##   0.03125977   0.15629884   0.06251954   0.09377931   0.03125977   0.25007815 
##     20200312     20200313     20200314     20200315     20200316     20200317 
##   0.25007815   0.31259769   0.59393560   0.87527352   0.62519537   0.84401375 
##     20200318     20200319     20200320     20200321     20200322     20200323 
##   1.75054705   1.93810566   1.75054705   0.93779306   1.21913098   1.03157237 
##     20200324     20200325     20200326     20200327     20200328     20200329 
##   2.00062520   3.12597687   1.62550797   1.06283214   4.25132854   4.43888715 
##     20200330     20200331     20200401     20200402     20200403     20200404 
##   4.06376993   7.81494217   3.65739294   6.93966865   9.12785245  26.03938731 
##     20200405     20200406     20200407     20200408     20200409     20200410 
##  13.87933729  20.22507033  32.57267896  33.85432948  25.53923101  32.44763989 
##     20200411     20200412     20200413     20200414     20200415     20200416 
##  43.88871522  23.06970928  40.13754298  32.54141919  43.23226008  41.54423257 
##     20200417     20200418     20200419     20200420     20200421     20200422 
##  34.66708346  43.98249453  33.63551110  63.73866833  67.86495780  65.77055330 
##     20200423     20200424     20200425     20200426     20200427     20200428 
##  64.27008440  69.99062207  74.92966552  44.45139106 106.25195374 105.25164114 
##     20200429     20200430     20200501     20200502     20200503     20200504 
## 132.47889966 105.40793998  82.05689278 105.31416068  60.08127540 124.72647702 
##     20200505     20200506     20200507     20200508     20200509     20200510 
## 124.07002188 113.59799937 112.91028446 103.06345733  77.18036887  29.25914348 
##     20200511     20200512     20200513     20200514     20200515     20200516 
## 112.56642701 135.91747421 138.82463270 116.88027509 134.54204439 128.94654580 
##     20200517     20200518     20200519     20200520     20200521     20200522 
##  55.79868709 160.36261332 162.48827759 165.36417631 172.36636449 190.68458893 
##     20200523     20200524     20200525     20200526     20200527     20200528 
## 168.77149109  56.01750547 217.81806815 216.28633948 230.41575492 200.84401375 
##     20200529     20200530     20200531     20200601     20200602     20200603 
## 179.58737105 144.20131291  70.70959675 132.60393873 117.97436699  81.71303532 
##     20200604     20200605     20200606     20200607     20200608     20200609 
## 125.94560800 154.54829634 148.92153798  82.43201000 148.26508284 165.98937168 
##     20200610     20200611     20200612     20200613     20200614     20200615 
## 178.99343545 158.45576743 145.70178181 140.66895905  76.30509534 136.38637074 
##     20200616     20200617     20200618     20200619     20200620     20200621 
## 105.53297906 122.38199437 107.40856518 107.75242263 110.22194436  45.20162551 
##     20200622     20200623     20200624     20200625     20200626     20200627 
## 135.51109722 117.78680838 128.44638950 129.85307909 135.38605814 110.22194436 
##     20200628     20200629     20200630     20200701     20200702     20200703 
##  49.95311035  85.30790872 117.16161300 114.03563614 104.68896530 120.72522663 
##     20200704     20200705     20200706     20200707     20200708     20200709 
## 109.97186621  71.86620819 118.25570491 126.13316661  81.21287902  77.99312285 
##     20200710     20200711     20200712     20200713     20200714     20200715 
## 153.57924351 132.72897781  74.39824945  97.34291966 131.10346983 138.54329478 
##     20200716     20200717     20200718     20200719     20200720     20200721 
## 141.76305095 160.70647077 158.36198812  69.64676461 117.41169115 150.51578618 
##     20200722     20200723     20200724     20200725     20200726     20200727 
## 186.68333854 171.11597374 186.30822132 175.86745858  92.21631760 221.03782432 
##     20200728     20200729     20200730     20200731     20200801     20200802 
##  94.31072210 210.09690528 226.57080338 237.60550172 209.34667083  99.84370116 
##     20200803     20200804     20200805     20200806     20200807     20200808 
## 244.38887152 242.38824633 234.69834323 263.89496718 227.00844014 204.28258831 
##     20200809     20200810     20200811     20200812 
## 104.06376993 224.19506096 237.26164426 158.45576743

Podemos comparar estos valores con la columna new_cases_per_million de la base de datos de Our World in Data y notar que difieren un poco, pero son muy similares.

owid_data <- read.csv("owid-covid-data.csv")
peru_positivos_owid <- owid_data[owid_data$location=="Peru", ]
peru_positivos_owid$new_cases_per_million
##   [1]      NA      NA      NA   0.030      NA   0.182   0.061   0.061   0.182
##  [10]   0.152   0.485   0.152   0.849   0.455   0.940   0.849   2.699   0.880
##  [19]   1.668   1.365   0.971   0.637   4.307   0.667   1.668   1.092   5.490
##  [28]   2.972   3.488   7.825   2.760   5.490   4.580  16.226   8.492  11.919
##  [37]  42.097  27.721  19.441  28.843  20.351   0.000  84.436  35.545  30.814
##  [46]  30.268  28.236  36.637  21.139  45.857  42.855  50.467  22.261 111.701
##  [55]  66.299  35.849  75.549  83.132  92.352 105.636  62.932 102.936  43.795
##  [64] 115.765 110.033 112.490 100.722  96.082  69.514  45.948  98.175 128.807
##  [73] 130.354 118.010 122.711 113.187  80.675 137.997 137.602 144.032  88.833
##  [82] 123.014 127.533 121.922 175.058 186.644 178.152 197.320 224.009 267.046
##  [91] 168.720   0.000 269.169 129.929 127.442 132.173 144.275  96.476 122.529
## [100] 154.283 180.912 180.791 132.932 139.634  98.751 126.290 113.794 105.545
## [109] 107.273 103.513 109.123   0.000 178.152 117.646 118.677 114.097 109.942
## [118] 104.028  89.349  86.377  98.994 106.970 109.032 105.575 110.337  90.532
## [127] 108.426 110.185 107.273  96.992  92.928 109.669 115.159 113.551 116.979
## [136] 117.130 119.830 120.193 124.045 124.076 133.629 135.358 137.875 147.550
## [145] 118.980 149.006 149.218 160.379 172.208 206.510 219.672 225.890 202.203
## [154] 128.898 205.933 234.564 236.111 256.765 216.457 212.666 154.950 198.563
## [163] 269.169 286.335 251.730

Entre los motivos de esta diferencia están que difieren en fechas, en el número usado como población total del Perú (en nuestro cálculo solo se usó un aproximado, en realidad hay más de 32 millones de habitantes) y diferencias en el número de casos reportados.

Solo se explicó la transformación de datos más usada en este contexto que es ajustar a la población total, pero hay muchas otras formas de transformar nuestros valores, por ejemplo, podemos obtener el número de casos ajustado al número total de camas UCI o algún otro valor que queramos considerar.

Fechas

Nos acercamos al final y solo falta cómo realizar gráficos básicos, pero antes de pasar a eso es importante explicar cómo trabajar con fechas.

Preparamos los datos como lo hemos estado haciendo.

peru_positivos <- read.csv("peru_positivos_covid.csv")
owid_data <- read.csv("owid-covid-data.csv")

Para este ejemplo usaremos los casos nuevos por día de Perú con la función table() de paso que recordamos.

peru_positivos_freqs <- table(peru_positivos$FECHA_RESULTADO)

El problema es que el formato en el que están las fechas no es el correcto y el programa no las reconoce como fechas. Esto es importante para realizar gráficos en los que usemos a la fecha en el eje ‘x’. En realidad hay diferentes formas de solucionar este problema, mostraré la que considero más sencilla.

En primer lugar convertiremos nuestras frecuencias (lo que hicimos con la función table()) en un data frame, con el que estamos ya familiarizados (es el formato en el que están los data sets). Para esto usamos la función as.data.frame().

peru_positivos_freqs <- as.data.frame(peru_positivos_freqs)

Podemos ver que ahora está en el formato que ya conocemos:

View(peru_positivos_freqs)

Notemos que tenemos dos columnas: una en la que tenemos las fechas y otra en la que tenemos las frecuencias de cada fecha (casos nuevos al día). Lo que debemos hacer es hacer que la columna de fechas sea reconocida por el programa como tales, para eso usaremos la función as.Date().

Para que funcione tenemos que precisar el formato en el que se encuentra la fecha, si observamos, nuestras fechas están como un número en el cual está el año, el mes y el día sin espacio entre ellos. Esto es un formato "%Y%m%d". Más información sobre esto se encuentra en el siguiente blog: Dates and Times in R Without Losing Your Sanity

Observemos cómo funciona.

as.Date(peru_positivos_freqs$Var1, format = "%Y%m%d")
##   [1] "2020-03-06" "2020-03-07" "2020-03-08" "2020-03-09" "2020-03-10"
##   [6] "2020-03-11" "2020-03-12" "2020-03-13" "2020-03-14" "2020-03-15"
##  [11] "2020-03-16" "2020-03-17" "2020-03-18" "2020-03-19" "2020-03-20"
##  [16] "2020-03-21" "2020-03-22" "2020-03-23" "2020-03-24" "2020-03-25"
##  [21] "2020-03-26" "2020-03-27" "2020-03-28" "2020-03-29" "2020-03-30"
##  [26] "2020-03-31" "2020-04-01" "2020-04-02" "2020-04-03" "2020-04-04"
##  [31] "2020-04-05" "2020-04-06" "2020-04-07" "2020-04-08" "2020-04-09"
##  [36] "2020-04-10" "2020-04-11" "2020-04-12" "2020-04-13" "2020-04-14"
##  [41] "2020-04-15" "2020-04-16" "2020-04-17" "2020-04-18" "2020-04-19"
##  [46] "2020-04-20" "2020-04-21" "2020-04-22" "2020-04-23" "2020-04-24"
##  [51] "2020-04-25" "2020-04-26" "2020-04-27" "2020-04-28" "2020-04-29"
##  [56] "2020-04-30" "2020-05-01" "2020-05-02" "2020-05-03" "2020-05-04"
##  [61] "2020-05-05" "2020-05-06" "2020-05-07" "2020-05-08" "2020-05-09"
##  [66] "2020-05-10" "2020-05-11" "2020-05-12" "2020-05-13" "2020-05-14"
##  [71] "2020-05-15" "2020-05-16" "2020-05-17" "2020-05-18" "2020-05-19"
##  [76] "2020-05-20" "2020-05-21" "2020-05-22" "2020-05-23" "2020-05-24"
##  [81] "2020-05-25" "2020-05-26" "2020-05-27" "2020-05-28" "2020-05-29"
##  [86] "2020-05-30" "2020-05-31" "2020-06-01" "2020-06-02" "2020-06-03"
##  [91] "2020-06-04" "2020-06-05" "2020-06-06" "2020-06-07" "2020-06-08"
##  [96] "2020-06-09" "2020-06-10" "2020-06-11" "2020-06-12" "2020-06-13"
## [101] "2020-06-14" "2020-06-15" "2020-06-16" "2020-06-17" "2020-06-18"
## [106] "2020-06-19" "2020-06-20" "2020-06-21" "2020-06-22" "2020-06-23"
## [111] "2020-06-24" "2020-06-25" "2020-06-26" "2020-06-27" "2020-06-28"
## [116] "2020-06-29" "2020-06-30" "2020-07-01" "2020-07-02" "2020-07-03"
## [121] "2020-07-04" "2020-07-05" "2020-07-06" "2020-07-07" "2020-07-08"
## [126] "2020-07-09" "2020-07-10" "2020-07-11" "2020-07-12" "2020-07-13"
## [131] "2020-07-14" "2020-07-15" "2020-07-16" "2020-07-17" "2020-07-18"
## [136] "2020-07-19" "2020-07-20" "2020-07-21" "2020-07-22" "2020-07-23"
## [141] "2020-07-24" "2020-07-25" "2020-07-26" "2020-07-27" "2020-07-28"
## [146] "2020-07-29" "2020-07-30" "2020-07-31" "2020-08-01" "2020-08-02"
## [151] "2020-08-03" "2020-08-04" "2020-08-05" "2020-08-06" "2020-08-07"
## [156] "2020-08-08" "2020-08-09" "2020-08-10" "2020-08-11" "2020-08-12"

Ahora, para incluir estas fechas en nuestro data set podemos reemplazar la columna Var1 o podemos crear una nueva. En este caso se creará una columna nueva con el nombre fechas de la siguiente manera (notar cómo creamos una columna):

peru_positivos_freqs$fechas <- as.Date(peru_positivos_freqs$Var1, format = "%Y%m%d")

Ahora observemos nuestro data set.

View(peru_positivos_freqs)

En el caso de nuestro data set de Our World in Data reemplazaremos la columna date por sus valores convertidos a fechas.

owid_data$date <- as.Date(owid_data$date, format = "%Y-%m-%d")

Es así como obtenemos alguna columna con fechas reconocidas por el programa, esto será de gran utilidad al realizar gráficos, como ya lo veremos.

Gráficos de barras

Los gráficos de barras se generan con la función barplot() y están estrechamente relacionados con la función table(), de hecho, podemos expresar gráficamente los resultados de la función table() usando barplot().

Por ejemplo, para obtener un gráfico de barras comparando los casos positivos de COVID-19 en el Perú entre hombres y mujeres es tan sencillo como generar la tabla y colocarla dentro de la función barplot().

peru_positivos <- read.csv("peru_positivos_covid.csv")
peru_sexos <- table(peru_positivos$SEXO)
peru_sexos
## 
##  FEMENINO MASCULINO 
##    227825    280171
barplot(peru_sexos)

En realidad se puede hacer esto con la frecuencia de lo que nos imaginemos, por ejemplo, edades.

peru_edades <- table(peru_positivos$EDAD)
barplot(peru_edades)

O departamentos:

peru_departamentos <- table(peru_positivos$DEPARTAMENTO)
# los nombres se sobreponen y no aparecen todos
barplot(peru_departamentos)

# arreglando los nombres de los departamentos
barplot(peru_departamentos, cex.names=0.7, las=2)

También podemos obtener gráficos de barras sin usar table(), para esto tenemos que definir las alturas de las barras.

x <- c(2, 8, 5, 4)
barplot(x)

Esta vez usaremos datos de países en el data set de Our World in Data para comparar las muertes totales ajustadas a la población de los siguiente países (hay muchos países): Perú, Estados Unidos, Brasil y China.

Primero filtramos los datos por países como ya sabemos (no se olviden de la comaaa).

owid_data <- read.csv("owid-covid-data.csv")
peru <- owid_data[owid_data$location == "Peru", ]
usa <- owid_data[owid_data$location == "United States", ]
brazil <- owid_data[owid_data$location == "Brazil", ]
china <- owid_data[owid_data$location == "China", ]

Luego podemos extraer el máximo valor que corresponde a muertes totales ajustadas a la poblacion (ya que es el actual) usando la función max(). Estas serán nuestras alturas. Usamos na.rm = TRUE para remover los valores faltantes.

peru_adj_deaths <- max(peru$total_deaths_per_million, na.rm = TRUE)
usa_adj_deaths <- max(usa$total_deaths_per_million, na.rm = TRUE)
brazil_adj_deaths <- max(brazil$total_deaths_per_million, na.rm = TRUE)
china_adj_deaths <- max(china$total_deaths_per_million, na.rm = TRUE)

Podemos guardar estos valores como un vector alturas y luego crear el gráfico de barras.

alturas <- c(peru_adj_deaths, usa_adj_deaths, brazil_adj_deaths, china_adj_deaths)
barplot(alturas)

# agregando nombres y ajustando eje y
barplot(alturas, names.arg=c("Perú","Estados Unidos","Brasil","China"), ylim=c(0,900))

Líneas de tiempo

Las líneas de tiempo que se suelen observar, por ejemplo, en el contexto del COVID-19, suelen representar la frecuencia de algún evento (casos confirmados totales/nuevos, muertes confirmadas totales/nuevas) a lo largo del tiempo. Es decir tenemos al tiempo en el eje ‘x’ y al número de eventos en el eje ‘y’.

Hay diferentes maneras de hacer esto en R, se explicará cómo hacerlo usando la función plot().

x <- c(1,2,3,4,5)
y <- c(2,2,4,1,5)
plot(x,y)

Como vimos, solo es cuestión de definir los ejes y usar plot(), pero para hacerlo correctamente debemos evocar muchas de las cosas aprendidas. Lo haremos con ejemplos.

Se presentarán los siguientes gráficos: casos positivos nuevos por día en Arequipa (con media móvil) y casos positivos totales en India (con función logarítmica). Lógicamente se puede realizar lo mismo para otros países y para otros valores, eso ya queda para que usen su creatividad.

Casos positivos nuevos por día en Arequipa con media móvil

Primero leemos los datos,

peru_positivos <- read.csv("peru_positivos_covid.csv")

filtramos los datos que corresponden a Arequipa (no olvidar la coma),

aqp_positivos <- peru_positivos[peru_positivos$DEPARTAMENTO == "AREQUIPA", ]

luego obtenemos las frecuencias usando table,

freqs <- table(aqp_positivos$FECHA_RESULTADO)

y hacemos que el programa reconozca las fechas.

freqs <- as.data.frame(freqs)
freqs$fechas <- as.Date(freqs$Var1, format="%Y%m%d")

Podemos visualizar que todo este correcto.

View(freqs)

Ahora definimos los ejes y graficamos.

x <- freqs$fechas
y <- freqs$Freq
plot(x,y)

Vemos que obtenemos puntos, para obtener una línea en lugar de puntos solo debemos agregar el argumento type = "l"

plot(x,y, type = "l")

En teoría ya tenemos nuestra gráfica lista, pero tiene un aspecto extraño, hay días en los que se reportan muchos casos y otros días que no. Una solución a esto es la de usar una media móvil, esto quiere decir obtener una lista con los valores promediando los valores cercanos. Puede leerse esto en este artículo de Wikipedia: Media móvil.

Podemos crear la función mav() que nos permite hacer esto. Notemos que \(n\) está predefinido con el valor de 5, lo que nos dice que se creará un promedio con 5 valores cercanos, podemos ajustar esto a nuestro gusto.

mav <- function(x,n=5){as.vector(stats::filter(x,rep(1/n,n),sides=2))}

Veamos cómo funciona:

z <- c(1,20,2,30,23,4,15,24,40,2,3,20)
mav(z)
##  [1]   NA   NA 15.2 15.8 14.8 19.2 21.2 17.0 16.8 17.8   NA   NA

Podemos usar esto en nuestros casos positivos nuevos para hacer que los cambios no sean tan notorios de un día a otro. Notemos cómo mejora nuestro gráfico.

Media móvil con 5 días:

y_mav5 <- mav(y)
plot(x,y_mav5, type="l")

Media móvil con 7 días:

y_mav7 <- mav(y, n=7)
plot(x,y_mav7, type="l")

Con \(n\) decimos los días alrededor de los cuales obtenemos un promedio. Our World in Data usa una media móvil con 7 días, queda en nosotros decidir cuántos días usamos.

Casos positivos totales en India con función logarítmica

Leemos,

owid_data <- read.csv("owid-covid-data.csv")

filtramos (no olviden la coma),

india <- owid_data[owid_data$location == "India", ]

transformamos las fechas,

india$date <- as.Date(india$date, format = "%Y-%m-%d")

visualizamos.

View(india)

Ahora, definimos los ejes y graficamos.

x <- india$date
y <- india$total_cases
plot(x,y, type="l")

Nuestra gráfica ya está lista, pero hay un detalle que suele considerarse al momento de hacer este tipo de gráficos, resulta que al inicio de una epidemia (o pandemia) el crecimiento del número de casos o el número de muertes es exponencial, como podemos observar en la gráfica que obtuvimos. Comparemos esto con la gráfica de Italia por ejemplo.

italia <- owid_data[owid_data$location=="Italy",]
italia$date <- as.Date(italia$date, format="%Y-%m-%d")
plot(italia$date,italia$total_cases, type="l")

Resulta que en estas situaciones suele preferirse usar el número de eventos en función logarítmica, para hacer esto tan solo debemos usar el argumento log y definir entre comillas el eje al cual aplicar la función logarítmica.

plot(x,y, type="l", log = "y")
## Warning in xy.coords(x, y, xlabel, ylabel, log): 30 y values <= 0 omitted from
## logarithmic plot

El programa nos da una advertencia de que hay valores que se han omitido. Esto no es mucho problema, pero si se desea se puede omitir los primeros días, en los cuales la función logarítmica no expresa los valores de forma adecuada.

Por ejemplo, al graficar solo a partir de que hubieron más de 10 casos podemos observar que la parte inicial del gráfico mejora.

y1 <- y[india$total_cases > 10]
x1 <- x[india$total_cases > 10]
plot(x1,y1, type="l", log="y")

Se puede obtener estos gráficos para diferentes países y de diferentes variables, todo depende del objetivo de su exploración y del objetivo del estudio que estén realizando.

Combinando gráficos

Este es el último tema con manejo de datos, así que es el último esfuerzo.

Puede ser de utilidad el comparar distintas poblaciones o distintos valores en una misma población, para esto puede ser conveniente visualizarlo en un gráfico.

Puede ser de utilidad comparar entre diferentes poblaciones, una forma muy sencilla de hacer esto es ajustar los valores a la población total. Ya vimos anteriormente cómo podemos hacer esto, es momento de aplicarlo para obtener gráficos.

Primero trabajaremos con la base de datos de positivos de la plataforma de Datos Abiertos del MINSA y luego desarrollaremos otro ejemplo para la base de datos de Our World in Data.

Casos totales: Arequipa vs. Perú

Esta vez graficaremos los casos totales en el tiempo entre Arequipa y Perú usando la base de datos de positivos de la plataforma de Datos Abiertos del MINSA. En las comparaciones entre diferentes poblaciones es de gran importancia poder ajustar los valores al tamaño de la población y esto es algo que ya sabemos hacer, así que ahora lo aplicaremos. Además, sabemos que a la fecha, la curva de los casos en Arequipa y el Perú aumentan exponencialmente, por lo que el gráfico quedaría mejor si representamos los casos en función logarítmica.

Esta combinación es la más difícil que se me ocurrió y es al propósito, para que, en teoría, si somos capaces de realizar este gráfico, somos capaces de hacer muchos otros gráficos más simples, como gráficos comparativos de casos/muertes nuevos por día o gráficos usando la base de datos de Our World in Data que es más fácil de manejar.

No perdamos de vista que para realizar nuestro gráfico simplemente tenemos que obtener nuestros ejes: el eje ‘x’ para el tiempo y el eje ‘y’ para los casos positivos ajustados a la población.

Haremos esto aplicando muchas de las cosas que hemos aprendido:

Preparamos los datos,

peru_positivos <- read.csv("peru_positivos_covid.csv")
aqp_positivos <- peru_positivos[peru_positivos$DEPARTAMENTO == "AREQUIPA", ]

obtenemos los casos por día con la función table(),

peru_pos_freq <- table(peru_positivos$FECHA_RESULTADO)
aqp_pos_freq <- table(aqp_positivos$FECHA_RESULTADO)

convertimos estas tablas de frecuencia en el formato que ya conocemos,

peru_pos_freq <- as.data.frame(peru_pos_freq)
aqp_pos_freq <- as.data.frame(aqp_pos_freq)

manejamos las fechas,

peru_pos_freq$fechas <- as.Date(peru_pos_freq$Var1, format = "%Y%m%d")
aqp_pos_freq$fechas <- as.Date(aqp_pos_freq$Var1, format = "%Y%m%d")

visualizamos.

View(peru_pos_freq)
View(aqp_pos_freq)

Ahora, recordemos a la función cumsum() que vimos en el tema de frecuencias. La usamos para obtener la frecuencia absoluta, si deseamos podemos guardarla en una nueva columna o en alguna otra variable. Esta vez se creará una columna freq_acum.

peru_pos_freq$freq_acum <- cumsum(peru_pos_freq$Freq)
aqp_pos_freq$freq_acum <- cumsum(aqp_pos_freq$Freq)

Visualizemos nuevamente para no perdernos.

View(peru_pos_freq)
View(aqp_pos_freq)

Podemos crear otra columna con estas frecuencias ajustadas a la población. Esta vez la nombraremos freq_acum_adj.

# Perú
peru_pop <- 31.99 # millones aprox.
peru_pos_freq$freq_acum_adj <- peru_pos_freq$freq_acum / peru_pop
# Arequipa
aqp_pop <- 1.008 # millones aprox.
aqp_pos_freq$freq_acum_adj <- aqp_pos_freq$freq_acum / aqp_pop

Visualizamos.

View(peru_pos_freq)
View(aqp_pos_freq)

Lo logramos, ya tenemos nuestros ejes con los cuales graficaremos. Nuestro eje ‘x’ del tiempo está en la columna fechas y nuestro eje ‘y’ está en la columna freq_acum_adj. Ahora grafiquemos usando una función logarítmica para el eje ‘y’.

# Perú
x_peru <- peru_pos_freq$fechas
y_peru <- peru_pos_freq$freq_acum_adj
plot(x_peru, y_peru, type="l", log="y")

# Arequipa
x_aqp <- aqp_pos_freq$fechas
y_aqp <- aqp_pos_freq$freq_acum_adj
plot(x_aqp, y_aqp, type="l", log="y")

Para combinar tenemos que crear primero un gráfico y luego agregar el otro encima con la función lines(). Notemos que ya no tenemos que precisar type ni log en la función lines().

plot(x_peru, y_peru, type="l", log="y")
lines(x_aqp, y_aqp)

Para agregar color solo hacemos lo siguiente (esto lo veremos en el siguiente tema).

plot(x_peru, y_peru, type="l", log="y", col = "red")
lines(x_aqp, y_aqp, col = "blue")

Comparando muertes nuevas al día entre Perú, Brasil y Estados Unidos

Si se pudo hacer el gráfico anterior este será mucho más sencillo.

Leemos,

owid_data <- read.csv("owid-covid-data.csv")

manejamos fechas (para no tener que hacerlo después para cada país),

owid_data$date <- as.Date(owid_data$date, format = "%Y-%m-%d")

filtramos (no olviden la coma),

peru <- owid_data[owid_data$location == "Peru", ]
brasil <- owid_data[owid_data$location == "Brazil", ]
usa <- owid_data[owid_data$location == "United States", ]

obtenemos nuestros ejes ‘x’,

x_peru <- peru$date
x_brasil <- brasil$date
x_usa <- usa$date

obtenemos nuestros ejes ‘y’,

y_peru <- peru$new_deaths
y_brasil <- brasil$new_deaths
y_usa <- usa$new_deaths

y graficamos.

plot(x_peru, y_peru, type="l", col="indianred")
lines(x_brasil, y_brasil, col="gray")
lines(x_usa, y_usa, col="lightseagreen")

Podemos usar la media móvil para mejorar el gráfico (\(n=7\)) y listo.

mav <- function(x,n=7){as.vector(stats::filter(x,rep(1/n,n),sides=2))}
y_peru_m <- mav(y_peru)
y_brasil_m <- mav(y_brasil)
y_usa_m <- mav(y_usa)
plot(x_usa, y_usa_m, type="l", col="lightseagreen")
lines(x_brasil, y_brasil_m, col="gray")
lines(x_peru, y_peru_m, col="indianred")

Comparemos este gráfico con el mostrado en la página web de Our World in Data a la fecha (15 de agosto de 2020, día de Arequipa u.u ):


Se parecen?

Link: https://ourworldindata.org/coronavirus-data-explorer?zoomToSelection=true&country=PER~USA~BRA&deathsMetric=true&interval=smoothed&smoothing=7&pickerMetric=location&pickerSort=asc

En realidad con esto sería todo por este texto y por la serie de videos, lo aprendido puede ser utilizado de muy diversas maneras, como es de imaginar. En el siguiente video simplemente se explicará algunos detalles sobre cómo poner nuestros gráficos bonitos, si queremos incluirlos en un reporte en un estudio.

Personalización de gráficos

En temas pasados vimos cómo usar los gráficos de barras y las líneas de tiempo, así que ahora se mostrará cómo personalizarlos. En realidad la mejor forma de aprender es buscando en internet y viendo las funciones de ayuda de las funciones gráficas. Las funciones de ayuda pueden obtenerse con la función help() o buscándolo en la ventana ‘Help’ en RStudio.

help(barplot)
help(plot)

Otros parámetros gráficos pueden encontrarse en help(par).

help(par)

Podemos observar que son muchos, así que solo explicaré los más básicos y que probablemente lleguen a usar.

Gráficos de barras

En el video de gráfico de barras vimos cómo obtener lo siguiente:

peru_positivos <- read.csv("peru_positivos_covid.csv")
peru_sexos <- table(peru_positivos$SEXO)
barplot(peru_sexos)

Podemos personalizar esto de la siguiente forma (no olviden las comas):

barplot(
  peru_sexos,
  main = "Casos positivos de COVID-19 en el Perú según sexo",
  col = c("palevioletred", "royalblue"),
  ylim = c(0,300000),
  names.arg = c("Mujeres", "Hombres"),
  cex.names = 1.2
)

También vimos este gráfico:

peru_departamentos <- table(peru_positivos$DEPARTAMENTO)
barplot(peru_departamentos)

Los nombres de las barras se sobreponen y tapan a otros, podemos solucionar esto seleccionando menos departamentos o ajustando estos nombres (no olviden las comas):

barplot(
  peru_departamentos,
  cex.names=0.7,
  las=2,
  col = c("lightsalmon","lightblue","palegreen","khaki"),
  main = "Número de casos positivos de COVID-19 por departamento",
  ylim = c(0,250000)
)

Vimos un gráfico con las edades:

peru_edades <- table(peru_positivos$EDAD)
barplot(peru_edades)

En las edades suele crearse categorías, para crear estas categorías podemos usar la función cut() con la que definimos los límites de las categorías con el argumento breaks y las etiquetas de las categorías con el argumento labels (por supuesto se pueden elegir las categorías que uno prefiera).

peru_edades_cat <- cut(
  peru_positivos$EDAD, breaks=c(0,20,40,60,80,200),
  labels=c("< 20","20-40","40-60","60-80","> 80")
)

Podemos visualizar cómo es que nos da el resultado:

peru_edades_cat

Ahora podemos utilizar nuestra conocida función table() para obtener la frecuencia de las categorías y luego usar barplot() como lo hemos estado haciendo (recordemos que table() y barplot() funcionan bien juntos).

barplot(table(peru_edades_cat))

Y podemos personalizarlo como ya se vio en los anteriores gráficos de barras, abajo hay un ejemplo con varios cambios realizados al gráfico anterior. En realidad no es necesario realizar todo esto, es solo demostrativo.

barplot(
  table(peru_edades_cat),
  horiz = TRUE,
  las = 1,
  col = "lightblue",
  main = "Casos positivos de COVID - 19 por grupo etario en el Perú",
  xlab = "Número de casos positivos",
  ylab = "Edades (años)",
  cex.names = 0.8,
  cex.axis = 0.8,
  cex.lab = 1.1,
  cex.main = 1.3,
  mgp = c(3,0.5,0)
)

Líneas de tiempo

Las líneas de tiempo se personalizan de forma similar y ya personalizamos un poco anteriormente algunos gráficos, así que solo se presentará un ejemplo.

Ya vimos cómo obtener este gráfico:

owid_data <- read.csv("owid-covid-data.csv")
owid_data$date <- as.Date(owid_data$date, format = "%Y-%m-%d")
peru <- owid_data[owid_data$location == "Peru", ]
brasil <- owid_data[owid_data$location == "Brazil", ]
usa <- owid_data[owid_data$location == "United States", ]
mav <- function(x,n=7){as.vector(stats::filter(x,rep(1/n,n),sides=2))}
x_peru <- peru$date
x_brasil <- brasil$date
x_usa <- usa$date
y_peru <- mav(peru$new_deaths)
y_brasil <- mav(brasil$new_deaths)
y_usa <- mav(usa$new_deaths)
plot(x_usa, y_usa, type="l", col="lightseagreen")
lines(x_brasil, y_brasil, col="gray")
lines(x_peru, y_peru, col="indianred")

Ahora podemos acomodarlo un poco más, e incluso podemos delimitar las fechas que queremos graficar.

plot(
  x_usa, y_usa,
  type = "l",
  col = "lightseagreen",
  lwd = 2,
  xlab = "",
  ylab = "",
  las = 1,
  xlim = c(as.Date("2020-03-01"), as.Date("2020-08-15"))
)
lines(
  x_brasil, y_brasil,
  col="gray",
  lwd = 2
)
lines(
  x_peru, y_peru,
  col="indianred",
  lwd = 2
)
grid()
title(
  main = "Muertes confirmadas nuevas diarias por COVID-19 \n con una media móvil de 7 días",
  line = 0.8
)
legend(
  x = "topright",
  legend = c("Estados Unidos","Brasil","Perú"),
  col = c("lightseagreen","gray","indianred"),
  lty = 1,
  bty = "n",
  lwd = 2
)

Hay muchas más cosas que se pueden hacer, pero por ahora es suficiente. Si tienen alguna duda pueden usar la función de ayuda help() como ya se mencionó o buscar su problema en internet, la comunidad de R es grande y hay grandes probabilidades que el problema que tengan ya haya sido contestado en algún lugar.

Guardando gráficos

Para guardar los gráficos hay diversas formas, pero la más sencilla es la usando la opción ‘Export’ en la ventana ‘Plots’ en RStudio y de esta manera guardar el gráfico que tengamos en la pantalla.

Cómo continuar aprendiendo

Los gráficos son una herramienta muy poderosa para explorar información de diferente naturaleza, en los gráficos mostrados se usó los gráficos usando las funciones ya incluidas en R, sin embargo, actualmente se suele usar diferentes paquetes para realizar gráficos, entre ellos el más popular es ggplot2 que permite generar gráficos de muy buena calidad, sin embargo es necesario aprender a usar este paquete, ya que sus funciones usan una estructura diferente basada en el concepto de “grammar of graphics” (más información en su página web).

Ya sea que se aprenda a usar paquetes gráficos o se aprenda a usar las funciones gráficas ya existentes en R, saber más de cómo obtener los gráficos que uno quiere usando R es de gran utilidad para la exploración de datos. Aquí hay un ejemplo de cómo se pueden replicar prácticamente cualquier gráfico usando R: You can replicate almost any plot with R.

Referencias