Introducción a la asignatura y a R base

Cuadernos prácticos de Series Temporales del Grado en Estadística Aplicada (curso 2024-2025)

Author

Javier Álvarez Liébana

1 Introducción a la asignatura

2 ¿Qué es una serie temporal?

Durante la carrera es probable que hayas tratado con multitud de datos pero hay uno muy especial que trataremos en esta asignatura de manera diferente: las series temporales.

 

Vamos a cargar el fichero retiro_temp.csv donde tenemos los datos de temperaturas diarios (AEMET) desde 1980 hasta 2024 de la estación instalada en El Retiro (Madrid).

Code
library(readr) # de tidyverse
# en tidyverse, read_ en lugar de read.
# tendremos datos en formato tibble en lugar de data.frame
retiro <- read_csv(file = "./datos/retiro_temp.csv")
Rows: 16314 Columns: 8
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (2): nombre, provincia
dbl  (5): id_station, altitud, tmed, tmin, tmax
date (1): fecha

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Code
retiro
# A tibble: 16,314 × 8
   fecha      id_station nombre         provincia altitud  tmed  tmin  tmax
   <date>          <dbl> <chr>          <chr>       <dbl> <dbl> <dbl> <dbl>
 1 2000-01-01       3195 MADRID, RETIRO MADRID        667   5.4   0.3  10.4
 2 2000-01-02       3195 MADRID, RETIRO MADRID        667   5     0.3   9.6
 3 2000-01-03       3195 MADRID, RETIRO MADRID        667   3.5   0.1   6.9
 4 2000-01-04       3195 MADRID, RETIRO MADRID        667   4.3   1.4   7.2
 5 2000-01-05       3195 MADRID, RETIRO MADRID        667   0.6  -0.4   1.6
 6 2000-01-06       3195 MADRID, RETIRO MADRID        667   3.8  -1.1   8.8
 7 2000-01-07       3195 MADRID, RETIRO MADRID        667   6.2   0.6  11.7
 8 2000-01-08       3195 MADRID, RETIRO MADRID        667   5.4  -0.1  11  
 9 2000-01-09       3195 MADRID, RETIRO MADRID        667   5.5   3     8  
10 2000-01-10       3195 MADRID, RETIRO MADRID        667   4.8   1.8   7.8
# ℹ 16,304 more rows

¿Qué analizar de estos datos?

 

Podemos por ejemplo visualizar un boxplot de las temperaturas medias de cada día durante estos últimos 44 años…

Code
library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ purrr     1.0.2
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Code
ggplot(retiro) +
  geom_boxplot(aes(y = tmed)) +
  scale_y_continuous(labels =
                       scales::label_number(suffix = "ºC")) +
  theme_minimal() +
  labs(title = "Temperatura desde 1980 hasta 2024",
       x = "Cuatrimestre", y = "Temperatura media diaria")
Warning: Removed 19 rows containing non-finite outside the scale range
(`stat_boxplot()`).

… la densidad de la temperatura durante todo ese tiempo…

Code
ggplot(retiro) +
  geom_density(aes(x = tmed)) +
  scale_x_continuous(labels =
                       scales::label_number(suffix = "ºC")) +
  theme_minimal() +
  labs(title = "Temperatura desde 1980 hasta 2024",
       x = "Temperatura media diaria")
Warning: Removed 19 rows containing non-finite outside the scale range
(`stat_density()`).

… pero también podríamos querer relacionar la temperatura media con el mes (por ejemplo con una regresión)…

Code
ggplot(retiro |> 
         mutate(mes = as_factor(lubridate::month(fecha))) |> 
         summarise(mean_temp = mean(tmed, na.rm = TRUE),
                   .by = "mes")) +
  geom_col(aes(x = mes, y = mean_temp)) +
  theme_minimal() +
  labs(title = "Temperatura media por mes",
       x = "Mes", y = "ºC (media)")

… o analizar cómo la temperatura media va incrementándose en cada década

Code
ggplot(retiro |> 
         mutate(periodo =
                  if_else(fecha < as_date("1990-01-01"),
                          "1980-1990",
                          if_else(fecha < as_date("2000-01-01"),
                                  "1990-2000",
                                  if_else(fecha < as_date("2010-01-01"),
                                          "2000-2010",
                                          if_else(fecha < as_date("2020-01-01"),
                                          "2010-2020", "después de 2020")))))) +
  geom_boxplot(aes(x = periodo, y = tmed)) +
  theme_minimal() +
  labs(title = "Temperatura media según periodo",
       x = "periodo", y = "ºC (media)")
Warning: Removed 19 rows containing non-finite outside the scale range
(`stat_boxplot()`).

En todos ejemplos anteriores hemos analizado una variable continua (temperatura) en función de una variable discreta o de grupo (periodo, década, etc).

 

¿Pero y si queremos relacionarla con una variable temporal “continua” como es la propia fecha?

Code
ggplot(retiro) +
  geom_line(aes(x = fecha, y = tmed), linewidth = 0.3, alpha = 0.7) +
  theme_minimal() +
  labs(title = "Temperatura media como SERIE TEMPORAL",
       x = "t (fecha)", y = "ºC (media)")

Fíjate bien…¿qué elementos detectas?

Code
ggplot(retiro) +
  geom_line(aes(x = fecha, y = tmed), linewidth = 0.3, alpha = 0.7) +
  theme_minimal() +
  labs(title = "Temperatura media como SERIE TEMPORAL",
       x = "t (fecha)", y = "ºC (media)")

  • Tendencia: lo que ajustarías con un modelo clásico (por ejemplo, una regresión lineal) y representa el comportamiento global de la serie, algo así como un nivel base respecto al que la serie oscila.

(en nuestro caso: la temperatura global aumenta con el paso de los años)

Code
ggplot(retiro, aes(x = fecha, y = tmed)) +
  geom_line(linewidth = 0.3, alpha = 0.7) +
  geom_smooth(method = "lm", se = FALSE) +
  theme_minimal() +
  labs(title = "Temperatura media como SERIE TEMPORAL",
       x = "t (fecha)", y = "ºC (media)")
`geom_smooth()` using formula = 'y ~ x'
Warning: Removed 19 rows containing non-finite outside the scale range
(`stat_smooth()`).

  • Estacionalidad: al margen de esa tendencia general, si hacemos zoom, en muchas series podemos observar un patrón que se repite cada x unidades temporales. En el caso de la temperatura, hay un patrón anual: diciembre hace más frío que en agosto.
Code
ggplot(retiro |> 
         filter(between(fecha, as_date("2020-01-01"), as_date("2023-12-31"))),
                aes(x = fecha, y = tmed)) +
  geom_line(linewidth = 0.3, alpha = 0.7) +
  geom_smooth(method = "loess") +
  theme_minimal() +
  labs(title = "Temperatura media diaria de 2020 a 2023",
       x = "t (fecha)", y = "ºC (media)")
`geom_smooth()` using formula = 'y ~ x'
Warning: Removed 11 rows containing non-finite outside the scale range
(`stat_smooth()`).

  • Atípicos: como sucede siempre en estadística será importantísimo analizar y tratar los datos atípicos muy alejados de lo esperado. Por ejemplo, en nuestro caso, Filomena.
Code
ggplot(retiro |> 
         filter(between(fecha, as_date("2020-01-01"), as_date("2023-12-31"))) |>
         mutate(filomena = between(fecha, as_date("2020-12-25"), as_date("2021-01-22"))),
                aes(x = fecha, y = tmed)) +
  geom_line(linewidth = 0.3, alpha = 0.7) +
  geom_point(aes(alpha = filomena), color = "#991545") +
  scale_alpha_manual(values = c(0, 1)) +
  guides(alpha = "none") +
  theme_minimal() +
  labs(title = "Temperatura media diaria de 2020 a 2023",
       x = "t (fecha)", y = "ºC (media)")
Warning: Removed 11 rows containing missing values or values outside the scale range
(`geom_point()`).

  • Intervenciones: incluso podría suceder que la serie tuviese un corte o salto en su comportamiento. Por ejemplo, imagina que de repente el aparato de medición empieza a medir +25 grados de la temperatura real.
Code
ggplot(retiro |> 
         filter(between(fecha, as_date("2020-01-01"), as_date("2023-12-31"))) |>
         mutate(tmed = if_else(fecha <= "2021-12-31", tmed, tmed + 25))) +
  geom_line(aes(x = fecha, y = tmed), linewidth = 0.3, alpha = 0.7) +
  guides(alpha = "none") +
  theme_minimal() +
  labs(title = "Temperatura media diaria de 2020 a 2023",
       subtitle = "Error de +25ºC a partir de 2022",
       x = "t (fecha)", y = "ºC (media)")

2.1 Estacionariedad

En esta asignatura será fundamental un concepto: estacionariedad. Diremos que una serie es estacionaria si oscila de manera estable con una media y varianza constante.

3 Introducción a R (base)

3.1 Conceptos básicos

3.1.1 Instalación de R

El lenguaje R será nuestra gramática y ortografía (nuestras reglas de juego)

  • Paso 1: entra en https://cran.r-project.org/ y selecciona tu sistema operativo.

  • Paso 2: para Mac basta con que hacer click en el archivo .pkg, y abrirlo una vez descargado. Para sistemas Windows, debemos clickar en install R for the first time y después en Download R for Windows. Una vez descargado, abrirlo como cualquier archivo de instalación.

  • Paso 3: abrir el ejecutable de instalación.

Warning

Siempre que tengas que descargar algo de CRAN (ya sea el propio R o un paquete), asegúrate de tener conexión a internet.

3.1.2 Instalación de R Studio

RStudio será el Word que usaremos para escribir (lo que se conoce como un IDE: entorno integrado de desarrollo).

  • Paso 1: entra la web oficial de RStudio (ahora llamado Posit) y selecciona la descarga gratuita.

  • Paso 2: selecciona el ejecutable que te aparezca acorde a tu sistema operativo.

  • Paso 3: tras descargar el ejecutable, hay que abrirlo como otro cualquier otro y dejar que termine la instalación.

3.1.3 Scripts (documentos .R)

Un script será el documento en el que programamos, nuestro archivo .doc (aquí con extensión .R) donde escribiremos las órdenes. Para abrir nuestro primero script, haz click en el menú en File < New File < R Script.

Cuidado

Es importante no abusar de la consola: todo lo que no escribas en un script, cuando cierres, lo habrás perdido.

Cuidado

R es case-sensitive: es sensible a mayúsculas y minúsculas por lo que x y X representa variables distintas.

3.1.4 Ejecutando el primer script

Ahora tenemos una cuarta ventana: la ventana donde escribiremos nuestros códigos. ¿Cómo ejecutarlo?

  1. Escribimos el código a ejecutar.

  2. Guardamos el archivo .R haciendo click en Save current document.

  3. El código no se ejecuta salvo que se lo indiquemos. Tenemos tres opciones de ejecutar un script:

  • Copiar y pegar en consola.
  • Seleccionar líneas y Ctrl+Enter
  • Activar Source on save a la derecha de guardar: no solo guarda sino que ejecuta el código completo.

3.1.5 Sé organizado: proyectos

De la misma manera que en el ordenador solemos trabajar de manera ordenada por carpetas, en RStudio podemos hacer lo mismo para trabajar de manera eficaz creando proyectos.

Un proyecto será una «carpeta» dentro de RStudio, de manera que nuestro directorio raíz automáticamente será la propia carpeta de proyecto (pudiendo pasar de un proyecto a otro con el menu superior derecho).

Podemos crear uno en una carpeta nueva o en una carpeta ya existente.

3.1.6 Buenas prácticas

  • Tip 1: asignar, evaluar y comparar no es lo mismo. Si te has fijado en R estamos usando <- para asignar valores a variables. Usaremos = para evaluar argumentos en funciones y == para saber si dos elementos son iguales.
x <- 1 # asignar
x = 1 # evaluar
x == 1 # comparar
  • Tip 2: programa como escribes. Al igual que cuando redactas en castellano, acostúmbrate a incorporar espacios y saltos de línea paranoquedarteciego (es una buena práctica y no un requisito porque R no procesa los espacios)
x <- 1 # óptimo
x<-1 # regu
x<- 1 # peor (decídete)
  • Tip 3: no seas caótico, estandariza nombres, acostúmbrate siempre a hacerlo igual. El único requisito es que debe empezar siempre por una letra (y sin tildes). La forma más recomendable es la conocida como snake_case
variable_en_modo_snake_case
otraFormaMasDificilDeLeer
hay.gente.que.usa.esto
Incluso_Haygente.Caotica_que.NoMereceNuestraATENCION
  • Tip 4: facilita la lectura y escritura, pon márgenes. En Tools < Global Options puedes personalizar algunas opciones de RStudio. En Code < Display podemos indicarle en Show margin (no interacciona con el código).

  • Tip 5: el tabulador es tu mejor amigo. En RStudio tenemos una herramienta maravillosa: si escribes parte del nombre de una variable o función y tabulas, RStudio te autocompleta

  • Tip 6: ni un paréntesis soltero. Siempre que abras un paréntesis deberás cerrarlo. Para facilitar esta tarea entra en Tools < Global Options < Code < Display y activa la opción Rainbow parentheses

  • Tip 7: fíjate en el lateral izquierdo. No solo podrás ver la línea de código por la que vas sino que, en caso de estar cometiendo un error de sintaxis, el propio RStudio te avisará.

  • Tip 8: intenta trabajar siempre por proyectos (para esta clase, crea un script clase2.R en el proyecto que creamos en la anterior clase)

 

Ver más tips en https://r4ds.had.co.nz/workflow-basics.html#whats-in-a-name

3.2 Tipos de datos

¿Existen variables más allá de los números en la ciencia de datos? Piensa por ejemplo en los datos que podrías guardar de una persona:

  • La edad o el peso será un número.
edad <- 33
  • Su nombre será una cadena de texto (conocida como string o char).
nombre <- "javi"
  • A la pregunta «¿estás matriculado en la Facultad?» la respuesta será lo que llamamos una variable lógica (TRUE si está matriculado o FALSE en otro caso).
matriculado <- TRUE
  • Su fecha de nacimiento será precisamente eso, una fecha, un tipo de variable crucial en esta asignatura

Puedes ver más detalles básicos de variables en R en https://javieralvarezliebana.es/docencia/R-datascience/diapos/#/title-slide

3.2.1 Variables de fecha

Un tipo de datos muy especial: los datos de tipo fecha.

fecha_char <- "2021-04-21"

Parece una simple cadena de texto pero debería representar un instante en el tiempo. ¿Qué debería suceder si sumamos un 1 a una fecha?

fecha_char + 1
Error in fecha_char + 1: non-numeric argument to binary operator

Las fechas NO pueden ser texto: debemos convertir la cadena de texto a fecha.

 

Para trabajar con fechas usaremos el paquete {lubridate}, que deberemos instalar antes de poder usarlo.

install.packages("lubridate")

Una vez instalado, de todos los paquetes (libros) que tenemos, le indicaremos que nos cargue ese concretamente.

library(lubridate) # instala si no lo has hecho

Para convertir a tipo fecha usaremos la función as_date() del paquete {lubridate} (por defecto en formato yyyy-mm-dd)

 

# ¡no es una fecha, es un texto!
fecha_char + 1
Error in fecha_char + 1: non-numeric argument to binary operator
class(fecha_char)
[1] "character"
fecha <- as_date("2023-03-28")
fecha + 1
[1] "2023-03-29"
class(fecha)
[1] "Date"

En as_date() el formato de fecha por defecto es yyyy-mm-dd así si la cadena de texto no se introduce de manera adecuada…

as_date("28-03-2023")
Warning: All formats failed to parse. No formats found.
[1] NA

Para cualquier otro formato debemos especificarlo en el argumento opcional format = ... tal que %d representa días, %m meses, %Y en formato de 4 años y %y en formato de 2 años.

as_date("28-03-2023", format = "%d-%m-%Y")
[1] "2023-03-28"
as_date("28-03-23", format = "%d-%m-%y")
[1] "2023-03-28"
as_date("03-28-2023", format = "%m-%d-%Y")
[1] "2023-03-28"
as_date("28/03/2023", format = "%d/%m/%Y")
[1] "2023-03-28"

En dicho paquete tenemos funciones muy útiles para manejar fechas:

  • Con today() podemos obtener directamente la fecha actual.
today()
[1] "2024-09-05"
  • Con now() podemos obtener la fecha y hora actual
now()
[1] "2024-09-05 16:58:18 CEST"
  • Con year(), month() o day() podemos extraer el año, mes y día
fecha <- today()
year(fecha)
[1] 2024
month(fecha)
[1] 9

Amplia contenido

Tienes un resumen en pdf de los paquetes más importantes en la carpeta correspondiente en el campus

3.3 Vectores: concatenar

Cuando trabajamos con datos normalmente tendremos columnas que representan variables: llamaremos vectores a una concatenación de celdas (valores) del mismo tipo (lo que sería una columna de una tabla).

La forma más sencilla es con el comando c() (c de concatenar), y basta con introducir sus elementos entre paréntesis y separados por comas

edades <- c(32, 27, 60, 61)
edades
[1] 32 27 60 61
Tip

Un número individual x <- 1 (o bien x <- c(1)) es en realidad un vector de longitud uno –> todo lo que sepamos hacer con un número podemos hacerlo con un vector de ellos.

3.4 💻 Tu turno

Intenta realizar los siguientes ejercicios sin mirar las soluciones

📝 Define el vector x como la concatenación de los 5 primeros números impares. Calcula la longitud del vector

Code
# Dos formas
x <- c(1, 3, 5, 7, 9)
x <- seq(1, 9, by = 2)

length(x)

📝 Accede al tercer elemento de x. Accede al último elemento (sin importar la longitud, un código que pueda ejecutarse siempre). Elimina el primer elemento.

Code
x[3]
x[length(x)]
x[-1]

📝 Obtén los elementos de x mayores que 4. Calcula el vector 1/x y guárdalo en una variable.

Code
x[x > 4]
z <- 1/x
z

📝 Crea un vector que represente los nombres de 5 personas, de los cuales uno es desconocido.

Code
nombres <- c("Javi", "Sandra", NA, "Laura", "Carlos")
nombres

📝 Encuentra del vector x de ejercicios anteriores los elementos mayores (estrictos) que 1 Y ADEMÁS menores (estrictos) que 7. Encuentra una forma de averiguar si todos los elementos son o no positivos.

Code
x[x > 1 & x < 7]
all(x > 0)

📝 Dado el vector x <- c(1, -5, 8, NA, 10, -3, 9), ¿por qué su media no devuelve un número sino lo que se muestra en el código inferior?

x <- c(1, -5, 8, NA, 10, -3, 9)
mean(x)
[1] NA

📝 Dado el vector x <- c(1, -5, 8, NA, 10, -3, 9), extrae los elementos que ocupan los lugares 1, 2, 5, 6.

Code
x <- c(1, -5, 8, NA, 10, -3, 9)
x[c(1, 2, 5, 6)]
x[-2]

📝 Dado el vector x del ejercicio anterior, ¿cuales tienen un dato ausente? Pista: las funciones is.algo() comprueban si el elemento es tipo algo (tabula)

Code
is.na(x)

📝 Define el vector x como la concatenación de los 4 primeros números pares. Calcula el número de elementos de x menores estrictamente que 5.

Code
x[x < 5] 
sum(x < 5)

📝 Calcula el vector 1/x y obtén la versión ordenada (de menor a mayor) de las dos formas posibles

Code
z <- 1/x
sort(z)
z[order(z)]

📝 Encuentra del vector x los elementos mayores (estrictos) que 1 y menores (estrictos) que 6. Encuentra una forma de averiguar si todos los elementos son o no negativos.

Code
x[x > 1 & x < 7]
all(x > 0)

3.5 Primera base de datos

Cuando analizamos datos solemos tener varias variables de cada individuo: necesitamos una «tabla» que las recopile. La opción más inmediata son las matrices: concatenación de variables del mismo tipo e igual longitud.

Imagina que tenemos estaturas y pesos de 4 personas. ¿Cómo crear un dataset con las dos variables? La opción más habitual es usando cbind(): concatenamos (bind) vectores en forma de columnas (c)

estaturas <- c(150, 160, 170, 180)
pesos <- c(63, 70, 85, 95)
datos_matriz <- cbind(estaturas, pesos)
datos_matriz
     estaturas pesos
[1,]       150    63
[2,]       160    70
[3,]       170    85
[4,]       180    95

3.5.1 Primer intento: matrices

También podemos construir la matriz por filas con la función rbind() (concatenar - bind - por filas - rows), aunque lo recomendable es tener cada variable en columna e individuo en fila como luego veremos.

rbind(estaturas, pesos) # Construimos la matriz por filas
          [,1] [,2] [,3] [,4]
estaturas  150  160  170  180
pesos       63   70   85   95
  • Podemos comprobar las dimensiones con dim(), nrow() y ncol(): las matrices son un tipo de datos tabulados (organizados en filas y columnas)
dim(datos_matriz)
[1] 4 2
nrow(datos_matriz)
[1] 4
ncol(datos_matriz)
[1] 2

3.5.2 Segundo intento: data.frame

Las matrices tienen el mismo problema que los vectores: si juntamos datos de distinto tipo, se perturba la integridad del dato ya que los convierte (fíjate en el código inferior: las edades y los TRUE/FALSE los ha convertido a texto)

edades <- c(14, 24, NA)
soltero <- c(TRUE, NA, FALSE)
nombres <- c("javi", "laura", "lucía")
matriz <- cbind(edades, soltero, nombres)
matriz
     edades soltero nombres
[1,] "14"   "TRUE"  "javi" 
[2,] "24"   NA      "laura"
[3,] NA     "FALSE" "lucía"

De hecho al no ser números ya no podemos realizar operaciones aritméticas

matriz + 1
Error in matriz + 1: non-numeric argument to binary operator

Para poder trabajar con variables de distinto tipo tenemos en R lo que se conoce como data.frame: concatenación de variables de igual longitud pero que pueden ser de tipo distinto.

tabla <- data.frame(edades, soltero, nombres)
class(tabla)
[1] "data.frame"
tabla
  edades soltero nombres
1     14    TRUE    javi
2     24      NA   laura
3     NA   FALSE   lucía

Dado que un data.frame es ya un intento de «base de datos» las variables no son meros vectores matemáticos: tienen un significado y podemos (debemos) ponerles nombres que describan su significado

library(lubridate)
tabla <-
  data.frame("edad" = edades, "estado" = soltero, "nombre" = nombres,
             "f_nacimiento" = as_date(c("1989-09-10", "1992-04-01", "1980-11-27")))
tabla
  edad estado nombre f_nacimiento
1   14   TRUE   javi   1989-09-10
2   24     NA  laura   1992-04-01
3   NA  FALSE  lucía   1980-11-27

3.5.3 Intento final: tibble

Las tablas en formato data.frame tienen algunas limitaciones. La principal es que no permite la recursividad: imagina que definimos una base de datos con estaturas y pesos, y queremos una tercera variable con el IMC

data.frame("estatura" = c(1.7, 1.8, 1.6), "peso" = c(80, 75, 70),
           "IMC" = peso / (estatura^2))
Error in data.frame(estatura = c(1.7, 1.8, 1.6), peso = c(80, 75, 70), : object 'peso' not found

En adelante usaremos el formato tibble (data.frame mejorado) del paquete {tibble}

library(tibble)
datos_tb <- 
  tibble("estatura" = c(1.7, 1.8, 1.6), "peso" = c(80, 75, 70), "IMC" = peso / (estatura^2))
class(datos_tb)
[1] "tbl_df"     "tbl"        "data.frame"
datos_tb
# A tibble: 3 × 3
  estatura  peso   IMC
     <dbl> <dbl> <dbl>
1      1.7    80  27.7
2      1.8    75  23.1
3      1.6    70  27.3
datos_tb <-
  tibble("estatura" = c(1.7, 1.8, 1.6), "peso" = c(80, 75, 70), "IMC" = peso / (estatura^2))
datos_tb
# A tibble: 3 × 3
  estatura  peso   IMC
     <dbl> <dbl> <dbl>
1      1.7    80  27.7
2      1.8    75  23.1
3      1.6    70  27.3

Las tablas en formato tibble nos permitirá una gestión más ágil, eficiente y coherente de los datos, con 4 ventajas principales:

  • Metainformación: si te fijas en la cabecera, nos dice ya automáticamente el número de filas y columnas, y el tipo de cada variable

  • Recursividad: permite definir las variables secuencialmente (como hemos visto)

  • Consistencia: si accedes a una columna que no existe avisa con un warning

datos_tb$invent
Warning: Unknown or uninitialised column: `invent`.
NULL
  • Por filas: crear por filas (copiar y pegar de una tabla) con tribble()
tribble(~colA, ~colB,
        "a",   1,
        "b",   2)
# A tibble: 2 × 2
  colA   colB
  <chr> <dbl>
1 a         1
2 b         2
Tip

El paquete {datapasta} nos permite copiar y pegar tablas de páginas web y documentos sencillos

3.6 💻 Tu turno (tb/df)

Intenta realizar los siguientes ejercicios sin mirar las soluciones

📝 Carga del paquete {datasets} el conjunto de datos airquality (variables de la calidad del aire de Nueva York desde mayo hasta septiembre de 1973). ¿Es el conjunto de datos airquality de tipo tibble? En caso negativo, conviértelo a tibble (busca en la documentación del paquete en https://tibble.tidyverse.org/index.html).

Code
library(tibble)
class(datasets::airquality)
airquality_tb <- as_tibble(datasets::airquality)

📝 Una vez convertido a tibble obtén el nombre de las variables y las dimensiones del conjunto de datos. ¿Cuántas variables hay? ¿Cuántos días se han medido?

Code
names(airquality_tb)
ncol(airquality_tb)
nrow(airquality_tb)

📝 Filtra solo los datos de la quinta observación

Code
airquality_tb[Month == 8, ]

📝 Filtra solo los datos del mes de agosto. ¿Cómo indicarle que queremos solo las filas que cumplan una condición concreta? (pista: en realidad todo son vectores “formateados”)

Code
airquality_tb[Month == 8, ]

📝 Selecciona aquellos datos que no sean ni de julio ni de agosto.

Code
airquality_tb[Month != 7 & Month != 8, ]
airquality_tb[!(Month %in% c(7, 8)), ]

📝 Modifica el siguiente código para quedarte solo con las variable de ozono y temperatura (sin importar qué posición ocupen)

airquality_tb[, 3]

📝 Selecciona los datos de temperatura y viento de agosto.

Code
airquality_tb[Month == 8, c("Temp", "Wind")]

📝 Traduce a castellano el nombre de las variables.

Code
names(airquality_tb) <- c("ozono", "rad_solar", "viento", "temp", "mes", "dia") 

3.7 Comunicar: rmd y Quarto

Una de las principales fortalezas de R es la facilidad para generar informes, libros, webs, apuntes y hasta diapositivas (este mismo material por ejemplo). Para ello instalaremos antes

  • el paquete {rmarkdown} (para generar archivos .rmd)
install.packages("rmarkdown")
  • instalar Quarto (si ya conocías R, el «nuevo» .rmd ahora como .qmd)

Hasta ahora solo hemos programado en scripts (archivos .R) dentro de proyectos, pero en muchas ocasiones no trabajaremos solos y necesitaremos comunicar los resultados en diferentes formatos:

  • apuntes (para nosotros mismos)
  • diapositivas
  • web
  • informes

Para todo ello usaremos Quarto (ver más en https://ivelasq.quarto.pub/intro-to-quarto/)

Los archivos de extensión .qmd (o .rmd antes) nos permitirán fácilmente combinar:

  • Markdown: lenguaje tipado que nos permite crear contenido simple (tipo wordpress, con texto, negritas, cursivas, etc) con un diseño legible.

  • Matemáticas (latex): lenguaje para escribir notación matemática como \(x^2\) o \(\sqrt{y}\) o \(\int_{a}^{b} f(x) dx\)

  • Código y salidas: podremos no solo mostrar el paso final sino el código que has ido realizando (en R, Python, C++, Julia, …), con cajitas de código llamadas CHUNKS.

  • Imágenes, gráficas, tablas, estilos (css, js), etc.

La principal ventaja de realizar este tipo de material en Quarto/Rmarkdown es que, al hacerlo desde RStudio, puedes generar un informe o una presentación sin salirte del entorno de programación en el que estás trabajando

De esta forma podrás analizar los datos, resumirlos y a la vez comunicarlos con la misma herramienta. Recientemente el equipo de RStudio desarrolló Quarto, una versión mejorada de Rmarkdown (archivos .qmd), con un formato un poco más estético y simple. Tienes toda la documentación y ejemplos en https://quarto.org/

3.7.1 Uso de Quarto

Imágenes obtenidas de https://ivelasq.quarto.pub/intro-to-quarto/#/working-with-the-rstudio-visual-editor

3.7.2 Nuestro primer informe

Vamos a crear el primer fichero rmarkdown con Quarto con extensión .qmd. Para ello solo necesitaremos hacer click en

File << New File << Quarto Document

Tras hacerlo nos aparecerán varias opciones de formatos de salida:

  • archivo .pdf
  • archivo .html (recomendable): documento dinámico, permite la interacción con el usuario, como una «página web».
  • archivo .doc (nada recomendable)

De momento dejaremos marcado el formato HTML que viene por defecto, y escribiremos el título de nuestro documento. Tras ello tendremos nuestro archivo .qmd (ya no es un script .R como los que hemos abierto hasta ahora).

Deberías tener algo similar a la captura de la imagen con dos modos de edición: Source (con código, la opción recomendada hasta que lo domines) y Visual (más parecido a un blog)

Para ejecutar TODO el documento debes clickar Render on Save y darle a guardar.

Deberías haber obtenido una salida en html similar a esta (y se te ha generado en tu ordenador un archivo html)

Como se indicaba, tienes dos formas de trabajar: con código puro y algo parecido a un Notion (blog)

Imagen obtenida de https://ivelasq.quarto.pub/intro-to-quarto/#/working-with-the-rstudio-visual-editor

Un fichero .qmd se divide básicamente en tres partes:

  • Cabecera: la parte que tienes al inicio entre ---.

  • Texto: que podremos formatear y mejorar con negritas (escrito como negritas, con doble astérisco al inicio y final), cursivas (cursivas, con barra baja al inicio y final) o destacar nombres de funciones o variables de R. Puedes añadir ecuaciones como \(x^2\) (he escrito $x^2$, entre dólares).

  • Código R

La cabecera están en formato YAML y contiene los metadatos del documento

  • title y subtitle: el título/subtítulo del documento
  • author: autor del mismo
  • format: formato de salida (podremos personalizar)
    • theme: si tienes algún archivo de estilos
    • toc: si quieres índice o no
    • toc-location: posición del índice
    • toc-title: título del índice
  • editor: si estás en modo visual o source.
---
title: "prueba"
author: "javier álvarez liébana"
format:
  html:
    style: style.css
    toc: true
    toc-location: left
    toc-title: Índice
editor: visual
---

Respecto a la escritura solo hay una cosa importante: salvo que indiquemos lo contrario, TODO lo que vamos a escribir es texto (normal). No código R.

Vamos a empezar escribiendo una sección al inicio (# Intro y detrás por ej. la frase

Este material ha sido diseñado por el profesor Javier Álvarez Liébana, docente en la Universidad Complutense de Madrid

Además al Running Code le añadiremos una almohadilla #: las almohadillas FUERA DE CHUNKS nos servirán para crear epígrafes (secciones) en el documento

Para que el índice capture dichas secciones modificaremos la cabecera del archivo como se observa en la imagen (puedes cambiar la localización del índice y el título si quieres para probar).

Vamos a personalizar un poco el texto haciendo lo siguiente:

  • Vamos a añadir negrita al nombre (poniendo ** al inicio y al final).

  • Vamos añadir cursiva a la palabra material (poniendo _ al inicio y al final).

  • Vamos añadir un enlace https://www.ucm.es, asociándolo al nombre de la Universidad. Para ello el título lo ponemos entre corchetes y justo detrás el enlace entre paréntesis [«Universidad Complutense de Madrid»](https://www.ucm.es)

Para añadir código R debemos crear nuestras cajas de código llamadas chunks: altos en el camino en nuestro texto markdown donde podremos incluir código de casi cualquier lenguaje (y sus salidas).

 

Para incluir uno deberá de ir encabezado de la siguiente forma tienes un atajo Command + Option + I (Mac) o Ctrl + Shift + I (Windows)

Dentro de dicha cajita (que tiene ahora otro color en el documento) escribiremos código R como lo veníamos haciendo hasta ahora en los scripts.

Vamos por ejemplo a definir dos variables y su suma de la siguiente manera, escribiendo dicho código en nuestro .qmd (dentro de ese chunk)

# Código R
x <- 1
y <- 2
x + y
[1] 3

Los chunks pueden tener un nombre o etiqueta, de forma que podamos referenciarlos de nuevo para no repetir código.

En cada chunk aparecen dos botones:

  • botón de play: activa la ejecución y salida de ese chunk particular (lo puedes visualizar dentro de tu propio RStudio)

  • botón de rebobinar: activa la ejecución y salida de todos los chunk hasta ese (sin llegar a él)

 

Además podemos incluir código R dentro de la línea de texto (en lugar de mostrar el texto x ejecuta el código R mostrando la variable).

Los chunks podemos personalizarlos con opciones al inicio del chunk precedido de #|:

  • #| echo: false: ejecuta código y se muestra resultado pero no visualiza código en la salida.

  • #| include: false: ejecuta código pero no muestra resultado y no visualiza código en la salida.

  • #| eval: false: no ejecuta código, no muestra resultado pero sí visualiza código en la salida.

  • #| message: false: ejecuta código pero no muestra mensajes de salida.

  • #| warning: false: ejecuta código pero no muestra mensajes de warning.

  • #| error: true: ejecuta código y permite que haya errores mostrando el mensaje de error en la salida.

Estas opciones podemos aplicarlas chunk a chunk o fijar los parámetros de forma global con knitr::opts_chunk$set() al inicio del documento (dentro de un chunk).

Si queremos que aplique la opción a todos los chunks por defecto debemos incluirlo al final de la cabecera, como opciones de ejecución

---
title: "¡Hola!"
format: html
editor: visual
execute:
  echo: false
---

Además de texto y código podemos introducir lo siguiente:

  • Ecuaciones: puedes añadir además ecuaciones como \(x^2\) (he escrito $x^2$, la ecuación entre dólares).

  • Listas: puedes itemizar elementos poniendo *

* Paso 1: ...

* Paso 2: ...

  • Cross-references: puedes etiquetar partes del documento (la etiqueta se construye con {#nombre-seccion}) y llamarlas luego con [Sección](@nombre-seccion)

Por último, también podemos añadir pies de gráficas o imágenes añadiendo #| fig-cap: "..."

Fíjate que el caption está en el margen (por ejemplo). Puedes cambiarlo introduciendo ajustes en la cabecera (todo lo relativo a figuras empieza por fig-, y puedes ver las opciones tabulando). Tienes más información en https://quarto.org/

Por último puedes añadir un tema personalizado incluyendo un archivo de estilos (archivo en formato .scss o .css). Te he dejado uno en https://github.com/dadosdelaplace/docencia-R-master-bio-2324/tree/main/material.

Importante

El archivo de estilos debe estar en la misma carpeta que el archivo .qmd

También puedes hacerlo de manera sencilla añadiendo a los textos un poco de HTML. Por ejemplo, para personalizar el color de un texto va entre corchetes y justo tras el texto, entre llaves, las opciones de estilo

Esta palabra es [roja]{style="color:red;"} ...
... y esta [verde y en negrita]{style="color:green; font-weight: bold;"}

Esta palabra es roja

… y esta verde y en negrita

 

También puedes usar los bloques de llamada que por defecto son note, tip, warning, caution e important (aunque los puedes crear y personalizar). Para ello basta con usar :::{.callout-tipo} y el tipo que quieras

:::{.callout-tip}

Note that there are five types of callouts, including: 
`note`, `tip`, `warning`, `caution`, and `important`.

:::
Tip

Recuerda que los 5 tipos son note, tip, warning, caution e important.

Caution

Úsalos con cabeza, a veces mucho recursos estético puede marear.

Además {reticulate} nos permite crear chunks de python dentro de un Quarto en R (ver https://quarto.org/docs/computations/python.html para crear jupyter notebooks directamente desde Quarto)

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

install_python("3.9.12") # Instalar python en PC sino lo tienes

# Instalar paquetes de Python
reticulate::py_install("numpy")
reticulate::py_install("matplotlib")
import numpy as np
import matplotlib.pyplot as plt
r = np.arange(0, 2, 0.05)
theta = 2 * np.pi * r
fig, ax = plt.subplots(
  subplot_kw = {'projection': 'polar'} 
)
ax.plot(theta, r)
plt.show()

3.7.3 Ejemplo de entrega

Vamos a realizar un pequeño simulacro antes de la entrega usando el dataset starwars del paquete {dplyr}

library(dplyr)
starwars
# A tibble: 87 × 14
   name     height  mass hair_color skin_color eye_color birth_year sex   gender
   <chr>     <int> <dbl> <chr>      <chr>      <chr>          <dbl> <chr> <chr> 
 1 Luke Sk…    172    77 blond      fair       blue            19   male  mascu…
 2 C-3PO       167    75 <NA>       gold       yellow         112   none  mascu…
 3 R2-D2        96    32 <NA>       white, bl… red             33   none  mascu…
 4 Darth V…    202   136 none       white      yellow          41.9 male  mascu…
 5 Leia Or…    150    49 brown      light      brown           19   fema… femin…
 6 Owen La…    178   120 brown, gr… light      blue            52   male  mascu…
 7 Beru Wh…    165    75 brown      light      blue            47   fema… femin…
 8 R5-D4        97    32 <NA>       white, red red             NA   none  mascu…
 9 Biggs D…    183    84 black      light      brown           24   male  mascu…
10 Obi-Wan…    182    77 auburn, w… fair       blue-gray       57   male  mascu…
# ℹ 77 more rows
# ℹ 5 more variables: homeworld <chr>, species <chr>, films <list>,
#   vehicles <list>, starships <list>

En él tenemos diferentes variables de los personajes de Star Wars, con características de su pelo, piel, altura, nombre, etc.

Crea un documento .qmd con nombre, título, formato e índice. Cada ejercicio posterior será una subsección del documento. Ejecuta los chunks que consideres y comenta las salidas para responder a cada pregunta

Ejercicio 1. ¿Cuántos personajes hay guardados en la base de datos? ¿Cuántas características se han medido de cada uno?

Ejercicio 2. Extrae en dos variables distintas nombres y edades las variables correspondientes de la tabla. ¿De qué tipo es la variable nombre? ¿Y la variable birth_year?

Ejercicio 3. Obtén el vector de nombres de los personajes ordenados de mayores a jóvenes.

Ejercicio 4. Busca ayuda de la función unique(). Úsala para saber que modalidades tiene la variable cualitativa correspondiente al color de ojos. ¿Cuántos distintos hay?

Ejercicio 5. ¿Existe ALGÚN valor ausente en la variable de color ojos?

Ejercicio 6. Calcula la media y desviación típica de las variables de estatura y peso (cuidado con los ausentes). Define un nuevo tibble con esas dos variables e incorpora una tercera variable que se llame “IMC” que calcule el índice de masa corporal. Incorpora con $ $ la fórmula usada para el IMC.

3.8 Estructuras de control

Una estructura de control se compone de una serie de comandos orientados a decidir el camino que tu código debe recorrer

  • Si se cumple la condición A, ¿qué sucede?

  • ¿Y si sucede B?

  • ¿Cómo puedo repetir una misma expresión (dependiendo de una variable)?

Si has programado antes, quizás te sea familiar las conocidas como estructuras condicionales tales como if (blabla) {...} else {...} o bucles for/while (a evitar siempre que podamos).

Una de las estructuras de control más famosas son las conocidas como estructuras condicionales if.

SI (IF) un conjunto de condiciones se cumple (TRUE), entonces ejecuta lo que haya dentro de las llaves

Por ejemplo, la estructura if (x == 1) { código A } lo que hará será ejecutar el código A entre llaves pero SOLO SI la condición entre paréntesis es cierta (solo si x es 1). En cualquier otro caso, no hará nada.

Por ejemplo, definamos un vector de edades de 8 personas

edad <- c(14, 17, 24, 56, 31, 20, 87, 73)
edad < 18
[1]  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE

Nuestra estructura condicional hará lo siguiente: si existe algún menor de edad, imprimirá por pantalla un mensaje.

if (any(edad < 18)) { 
  
  print("Existe alguna persona menor de edad")
  
}
[1] "Existe alguna persona menor de edad"
if (any(edad < 18)) { 
  
  print("Existe alguna persona menor de edad")
  
}

En caso de que las condiciones no sean ciertas dentro de if() (FALSE), no sucede nada

if (all(edad >= 18)) { 
  
  print("Todos son mayores de edad")
  
}

No obtenemos ningún mensaje porque la condición all(edad >= 18) no es TRUE, así que no ejecuta nada.

La estructura if (condicion) { código A } puede combinarse con un else { código B }: cuando la condición no está verificada, se ejecutará el código alternativo B dentro de else { }, permitiéndonos decidir que sucede cuando se cumple y cuando no.

Por ejemplo, if (x == 1) { código A } else { código B } ejecutará A si x es igual a 1 y B en cualquier otro caso.

if (all(edad >= 18)) { 
  
  print("Todos son mayores de edad")
  
} else {
  
  print("Existe alguna persona menor de edad")
}
[1] "Existe alguna persona menor de edad"

Esta estructura if - else puede ser anidada: imagina que queremos ejecutar un código si todos son menores; si no sucede, pero todos son mayores de 16, hacer otra cosa; en cualquier otra cosa, otra acción.

if (all(edad >= 18)) { 
  
  print("Todos son mayores de edad")
  
} else if (all(edad >= 16)) {
  
  print("Hay algún menor de edad pero todos con 16 años o más")
  
} else { print("Hay alguna persona con menos de 16 años") }
[1] "Hay alguna persona con menos de 16 años"
Truco

Puedes colapsar las estructuras haciendo click en la flecha a la izquierda que aparece en tu script.

Esta estructura condicional se puede vectorizar (en una sola línea) con if_else() (del paquete {dplyr}), cuyos argumentos son

  • la condición a evaluar
  • lo que sucede cuando se cumple y cuando no
  • un argumento opcional para cuando la condición a evaluar es NA

Vamos a etiquetar sin son mayores/menores y un “desconocido” cuando no conocemos

library(dplyr)
edad <- c(NA, edad)
if_else(edad >= 18, "mayor", "menor", missing = "desconocido")
[1] "desconocido" "menor"       "menor"       "mayor"       "mayor"      
[6] "mayor"       "mayor"       "mayor"       "mayor"      

En R base existe ifelse(): no deja especificar que hacer con los ausentes pero permite especificar distintos tipos de datos en TRUE y en FALSE.

3.9 Bucles

Aunque en la mayoría de ocasiones se pueden reemplazar por otras estructuras más eficientes y legibles, es importante conocer una de las expresiones de control más famosas: los bucles.

  • for { }: permite repetir el mismo código en un número prefijado y conocido de veces.

  • while { }: permite repetir el mismo código pero en un número indeterminado de veces (hasta que una condición deje de cumplirse).

3.9.1 Bucles for

Un bucle for es una estructura que permite repetir un conjunto de órdenes un número finito, prefijado y conocido de veces dado un conjunto de índices.

Vamos a definir un vector x <- c(0, -7, 1, 4) y otra variable vacía y. Tras ello definiremos un bucle for con for () { }: dentro de los paréntesis indicaremos un índice y unos valores a recorrer, dentro de las llaves el código a ejecutar en cada iteración (en este caso, rellenar y como x + 1)

x <- c(0, -7, 1, 4)
y <- c()

for (i in 1:4) {
  y[i] <- x[i] + 1
}

Fíjate que debido a que R funciona de manera vectorial por defecto, el bucle es lo mismo que hacer x + 1 directamente.

x <- c(0, -7, 1, 4)
y <- c()

for (i in 1:4) {
  y[i] <- x[i] + 1
}
y
[1]  1 -6  2  5
y2 <- x + 1
y2
[1]  1 -6  2  5

Otra opción habitual es indicar los índices de manera «automática»: desde el primero 1 hasta el último (que corresponde con la longitud de x length(x))

x <- c(0, -7, 1, 4)
y <- c()

for (i in 1:length(x)) {
  y[i] <- x[i] + 1
}
y
[1]  1 -6  2  5

Así la estructura general de un bucle for será siempre la siguiente

for (índice in conjunto) { 
  código (dependiente de i)
}

SIEMPRE sabemos cuántas iteraciones tenemos (tantas como elementos haya en el conjunto a indexar). Como ya hemos aprendido con el paquete{microbenchmark} podemos chequear como los bucles suelen ser muy ineficientes (de ahí que debamos evitarlos en la mayoría de ocasiones

library(microbenchmark)
x <- 1:1000
microbenchmark(y <- x^2, 
               for (i in 1:100) { y[i] <- x[i]^2 },
               times = 500)
Warning in microbenchmark(y <- x^2, for (i in 1:100) {: less accurate
nanosecond times to avoid potential integer overflows
Unit: nanoseconds
                                    expr    min       lq       mean   median
                                y <- x^2    984   1148.0   1412.286   1230.0
 for (i in 1:100) {     y[i] <- x[i]^2 } 805240 827646.5 896646.302 846465.5
       uq     max neval
   1353.0    6478   500
 878076.5 3847317   500

Podemos ver otro ejemplo de bucle combinando números y textos: definimos un vector de edades y de nombres, e imprimimos el nombre y edad i-ésima.

nombres <- c("Javi", "Sandra", "Carlos", "Marcos", "Marta")
edades <- c(33, 27, 18, 43, 29)
library(glue)
for (i in 1:5) { 
  
  print(glue("{nombres[i]} tiene {edades[i]} años")) 
  
}
Javi tiene 33 años
Sandra tiene 27 años
Carlos tiene 18 años
Marcos tiene 43 años
Marta tiene 29 años

Aunque normalmente se suelen indexar con vectors numéricos, los bucles pueden ser indexados sobre cualquier estructura vectorial, da igual de que tipo sea el conjunto

library(stringr)
week_days <- c("monday", "tuesday", "wednesday", "thursday",
               "friday", "saturday", "sunday")

for (days in week_days) {
  
  print(str_to_upper(days))
}
[1] "MONDAY"
[1] "TUESDAY"
[1] "WEDNESDAY"
[1] "THURSDAY"
[1] "FRIDAY"
[1] "SATURDAY"
[1] "SUNDAY"

Vamos a combinar las estructuras condicionales y los bucles: usando el conjunto swiss del paquete {datasets}, vamos a asignar NA si los valores de fertilidad son mayores de 80.

for (i in 1:nrow(swiss)) {
  
  if (swiss$Fertility[i] > 80) { 
    
    swiss$Fertility[i] <- NA
    
  }
}

Esto es exactamente igual a un if_else() vectorizado

data("swiss")
swiss$Fertility <- if_else(swiss$Fertility > 80, NA, swiss$Fertility)

3.9.2 Bucles while

Otra forma de crear un bucle es con la estructura while { }, que nos ejecutará un bucle un número desconocido de veces, hasta que una condición deje de cumplirse (de hecho puede que nunca termine). Por ejemplo, vamos a inializar una variable ciclos <- 1, que incrementaremos en cada paso, y no saldremos del bucle hasta que ciclos > 4.

ciclos <- 1
while(ciclos <= 4) {
  
  print(glue("No todavía, vamos por el ciclo {ciclos}")) 
  ciclos <- ciclos + 1
  
}
No todavía, vamos por el ciclo 1
No todavía, vamos por el ciclo 2
No todavía, vamos por el ciclo 3
No todavía, vamos por el ciclo 4

Un bucle while será siempre como sigue

while(condición) {
  
  código a hacer mientras la condición sea TRUE
  # normalmente aquí se actualiza alguna variable
  
}

¿Qué sucede cuando la condición nunca es FALSE? Pruébalo tu mismo

while (1 > 0) {
  
  print("Presiona ESC para salir del bucle")
  
}

 

Cuidado

Un bucle while { } puede ser bastante «peligroso» sino controlamos bien cómo pararlo.

Contamos con dos palabras reservadas para abortar un bucle o forzar su avance:

  • break: permite abortar un bucle incluso si no se ha llegado a su final
for(i in 1:10) {
  if (i == 3) {
    
    break # si i = 3, abortamos bucle
    
  }
  print(i)
}
[1] 1
[1] 2

Contamos con dos palabras reservadas para abortar un bucle o forzar su avance:

  • next: fuerza un bucle a avanzar a la siguiente iteración
for(i in 1:5) {
  if (i == 3) {
    
    next # si i = 3, la obvia y continua al siguiente
    
  }
  print(i)
}
[1] 1
[1] 2
[1] 4
[1] 5

3.10 💻 Tu turno

Intenta realizar los siguientes ejercicios sin mirar las soluciones

📝 ¿Cuál es la salida del siguiente código?

if_else(sqrt(9) < 2, sqrt(9), 0)
Code
La salida es 0 ya que sqrt(9) es igual 3, y dado que no es menor que 2, devuelve el segundo argumento que es 0

📝 ¿Cuál es la salida del siguiente código?

x <- c(1, NA, -1, 9)
if_else(sqrt(x) < 2, 0, 1)
Code
La salida es el vector c(0, NA, NA, 1) ya que sqrt(1) sí es menor que 2, sqrt(9) no lo es, y tanto en el caso de sqrt(NA) (raíz de ausente) como sqrt(-1) (devuelve NaN, not a number), su raíz cuadrada no puede verificarse si es menor que 2 o no, así que la salida es NA.

📝 Modifica el código inferior para que, cuando no se pueda verificar si la raíz cuadrada de un número es menor que 2, devuelva -1

x <- c(1, NA, -1, 9)
if_else(sqrt(x) < 2, 0, 1)
Code
x <- c(1, NA, -1, 9)
if_else(sqrt(x) < 2, 0, 1, missing = -1)

📝 ¿Cuál es son los valores de x e y del código inferior para z <- 1, z <- -1 y z <- -5?

z <- -1
if (z > 0) {
  
  x <- z^3
  y <- -sqrt(z)
  
} else if (abs(z) < 2) {
  
  x <- z^4
  y <- sqrt(-z)
  
} else {
  
  x <- z/2
  y <- abs(z)
  
}
Code
En primero caso x = 1 e y = -1. En el segundo caso x = 1 e y = 1. En el tercer caso -1 y 2

📝 ¿Qué sucederá si ejecutamos el código inferior?

z <- "a"
if (z > 0) {
  
  x <- z^3
  y <- -sqrt(z)
  
} else if (abs(z) < 2) {
  
  x <- z^4
  y <- sqrt(-z)
  
} else {
  
  x <- z/2
  y <- abs(z)
  
}
Code
# dará error ya que no es un argumento numérico
Error in z^3 : non-numeric argument to binary operator

📝 Del paquete {lubridate}, la función hour() nos devuelve la hora de una fecha dada, y la función now() nos devuelve fecha y hora del momento actual. Con ambas funciones haz que se imprima por pantalla (cat()) “buenas noches” solo a partir de las 21 horas.

Code
# Cargamos librería
library(lubridate)

# Fecha-hora actual
fecha_actual <- now()

# Estructura if
if (hour(fecha_actual) > 21) {
  
  cat("Buenas noches") # print/cat dos formas de imprimir por pantalla
}

📝 Modifica el código inferior para que se imprima un mensaje por pantalla si y solo si todos los datos de airquality son con mes distinto a enero

library(datasets)
months <- airquality$Month

if (months == 2) {
  print("No hay datos de enero")
}
Code
library(datasets)
months <- airquality$Month

if (all(months != 1)) {
  print("No hay datos de enero")
}

📝 Modifica el código inferior para guardar en una variable llamada temp_alta un TRUE si alguno de los registros tiene una temperatura superior a 90 grados Farenheit y FALSE en cualquier otro caso

temp <- airquality$Temp

if (temp == 100) {
  print("Algunos de los registros tienen temperaturas superiores a 90 grados Farenheit")
}
Code
# Option 1
temp <- airquality$Temp
temp_alta <- FALSE
if (any(temp > 90)) {
   temp_alta <- TRUE
}

# Option 2
temp_alta <- any(airquality$Temp > 90)

📝 Modifica el código inferior para diseñar un bucle for de 5 iteraciones que solo recorra los primeros 5 impares (y en cada paso del bucle los imprima)

for (i in 1:5) {
  
  print(i)
}
Code
for (i in c(1, 3, 5, 7, 9)) {
  
  print(i)
}

📝 Modifica el código inferior para diseñar un bucle while que empiece con un contador count <- 1 y pare cuando llegue a 6

count <- 1
while (count == 2) {
  
  print(count)
}
Code
count <- 1
while (count < 6) {
  
  print(count)
  count <- count + 1
  
}

3.11 Funciones

No solo podemos usar funciones predeterminadas que vienen ya cargadas en paquetes, además podemos crear nuestras propias funciones para automatizar tareas. ¿Cómo crear nuestra propia función? Veamos su esquema básico:

  • Nombre: por ejemplo name_fun (sin espacios ni caracteres extraños). Al nombre le asignamos la palabra reservada function().

  • Definir argumentos de entrada (dentro de function()).

  • Cuerpo de la función dentro de { }.

  • Finalizamos la función con los argumentos de salida con return().

name_fun <- function(arg1, arg2, ...) {
  
  código a ejecutar
  
  return(var_salida)
  
}
  • arg1, arg2, ...: serán los argumentos de entrada, los argumentos que toma la función para ejecutar el código que tiene dentro

  • código: líneas de código que queramos que ejecute la función.

  • return(var_salida): se introducirán los argumentos de salida.

name_fun <- function(arg1, arg2, ...) {
  
  # Código que queramos ejecutar
  código
  
  # Salida
  return(var_salida)
  
}
Importante

Todas las variables que definamos dentro de la función son variables LOCALES: solo existirán dentro de la función salvo que especifiquemos lo contrario.

Veamos un ejemplo muy simple de función para calcular el área de un rectángulo.

Dado que el área de un rectángulo se calcula como el producto de sus lados, necesitaremos precisamente eso, sus lados: esos serán los argumentos de entrada y el valor a devolver será justo su área (\(lado_1 * lado_2\)).

# Definición del nombre de función y argumentos de entrada
calcular_area <- function(lado_1, lado_2) {
  
  area <- lado_1 * lado_2
  return(area)
  
}

También podemos hacer una definición directa de las variables sin almacenar por el camino.

# Definición del nombre de función y argumentos de entrada
calcular_area <- function(lado_1, lado_2) {
  
  return(lado_1 * lado_2)
  
}

¿Cómo aplicar la función?

calcular_area(5, 3) # área de un rectángulo 5 x 3 
[1] 15
calcular_area(1, 5) # área de un rectángulo 1 x 5
[1] 5
Tip

Aunque no sea necesario, es recomendable hacer explícita la llamada de los argumentos, especificando en el código qué valor es para cada argumento para que no dependa de su orden, haciendo el código más legible

calcular_area(lado_1 = 5, lado_2 = 3) # área de un rectángulo 5 x 3 
[1] 15
calcular_area(lado_2 = 3, lado_1 = 5) # área de un rectángulo 5 x 3 
[1] 15

3.11.1 Argumentos por defecto

Imagina ahora que nos damos cuenta que el 90% de las veces usamos dicha función para calcular por defecto el área de un cuadrado (es decir, solo necesitamos un lado). Para ello, podemos definir argumentos por defecto en la función: tomarán dicho valor salvo que le asignemos otro.

¿Por qué no asignar lado_2 = lado_1 por defecto, para ahorrar líneas de código y tiempo?

calcular_area <- function(lado_1, lado_2 = lado_1) {
  
  # Cuerpo de la función
  area <- lado_1 * lado_2
  
  # Resultado que devolvemos
  return(area)
  
}
calcular_area <- function(lado_1, lado_2 = lado_1) {
  
  # Cuerpo de la función
  area <- lado_1 * lado_2
  
  # Resultado que devolvemos
  return(area)
  
}

Ahora por defecto el segundo lado será igual al primero (si se lo añadimos usará ambos).

calcular_area(lado_1 = 5) # cuadrado
[1] 25
calcular_area(lado_1 = 5, lado_2 = 7) # rectángulo
[1] 35

3.11.2 Salida múltiple

Compliquemos un poco la función y añadamos en la salida los valores de cada lado, etiquetados como lado_1 y lado_2, empaquetando la salida en una vector.

# Definición del nombre de función y argumentos de entrada
calcular_area <- function(lado_1, lado_2 = lado_1) {
  
  # Cuerpo de la función
  area <- lado_1 * lado_2
  
  # Resultado
  return(c("area" = area, "lado_1" = lado_1, "lado_2" = lado_2))
  
}

Podemos complicar un poco más la salida añadiendo una cuarta variable que nos diga, en función de los argumentos, si rectángulo o cuadrado, teniendo que añadir en la salida una variable que de tipo caracter (o lógica).

# Definición del nombre de función y argumentos de entrada
calcular_area <- function(lado_1, lado_2 = lado_1) {
  
  # Cuerpo de la función
  area <- lado_1 * lado_2
  
  # Resultado
  return(c("area" = area, "lado_1" = lado_1, "lado_2" = lado_2,
           "tipo" = if_else(lado_1 == lado_2, "cuadrado", "rectángulo")))
  
}
calcular_area(5, 3)
        area       lado_1       lado_2         tipo 
        "15"          "5"          "3" "rectángulo" 

Problema: al intentar juntar números y texto, lo convierte todo a números. Podríamos guardarlo todo en un tibble() como hemos aprendido o en un objeto conocido en R como listas

3.11.3 Orden de los argumentos

Antes nos daba igual el orden de los argumentos pero ahora el orden de los argumentos de entrada importa, ya que en la salida incluimos lado_1 y lado_2.

Recomendación

Como se comentaba, altamente recomendable hacer la llamada a la función indicando explícitamente los argumentos para mejorar legibilidad e interpretabilidad.

# Equivalente a calcular_area(5, 3)
calcular_area(lado_1 = 5, lado_2 = 3)
        area       lado_1       lado_2         tipo 
        "15"          "5"          "3" "rectángulo" 

3.11.4 Variables locales vs globales

Un aspecto importante sobre el que reflexionar con las funciones: ¿qué sucede si nombramos a una variable dentro de una función a la que se nos ha olvidado asignar un valor dentro de la misma?

Debemos ser cautos al usar funciones en R, ya que debido a la «regla lexicográfica», si una variable no se define dentro de la función, R buscará dicha variable en el entorno de variables.

x <- 1
funcion_ejemplo <- function() {
    
  print(x) # No devuelve nada, solo realiza la acción 
}
funcion_ejemplo()
[1] 1

Si una variable ya está definida fuera de la función (entorno global), y además es usada dentro de cambiando su valor, el valor solo cambia dentro pero no en el entorno global.

x <- 1
funcion_ejemplo <- function() {
    
  x <- 2
  print(x) # lo que vale dentro
}
# lo que vale dentro
funcion_ejemplo() #<<
[1] 2
# lo que vale fuera
print(x) #<<
[1] 1

Si queremos que además de cambiar localmente lo haga globalmente deberemos usar la doble asignación (<<-).

x <- 1
y <- 2
funcion_ejemplo <- function() {
  
  # no cambia globalmente, solo localmente
  x <- 3 
  # cambia globalmente
  y <<- 0 #<<
  
  print(x)
  print(y)
}

funcion_ejemplo() # lo que vale dentro
[1] 3
[1] 0
x # lo que vale fuera
[1] 1
y # lo que vale fuera
[1] 0

3.12 💻 Tu turno

Intenta realizar los siguientes ejercicios sin mirar las soluciones

📝 Modifica el código inferior para definir una función llamada funcion_suma, de forma que dados dos elementos, devuelve su suma.

nombre <- function(x, y) {
  suma <- # código a ejecutar
  return()
}
# Aplicamos la función
suma(3, 7)
Code
funcion_suma <- function(x, y) {
  suma <- x + y
  return(suma)
}
funcion_suma(3, 7)

📝 Modifica el código inferior para definir una función llamada funcion_producto, de forma que dados dos elementos, devuelve su producto, pero que por defecto calcule el cuadrado

nombre <- function(x, y) {
  producto <- # código de la multiplicación
  return()
}
producto(3)
producto(3, -7)
Code
funcion_producto <- function(x, y = x) {
  producto <- x * y
  return(producto)
}
funcion_producto(3)
funcion_producto(3, -7)

📝 Define una función llamada igualdad_nombres que, dados dos nombres, nos diga si son iguales o no. Hazlo considerando importantes las mayúsculas, y sin que importen las mayúsculas. Usa el paquete {stringr}.

Code
# Distinguiendo mayúsculas
igualdad_nombres <- function(persona_1, persona_2) {
  return(persona_1 == persona_2)
}
igualdad_nombres("Javi", "javi")
igualdad_nombres("Javi", "Lucía")

# Sin importar mayúsculas
igualdad_nombres <- function(persona_1, persona_2) {
  return(toupper(persona_1) == toupper(persona_2))
}
igualdad_nombres("Javi", "javi")
igualdad_nombres("Javi", "Lucía")

📝 Crea una función llamada calculo_IMC que, dados dos argumentos (peso y estatura en metros) y un nombre, devuelva una lista con el IMC (\(peso/(estatura_m^2)\)) y el nombre.

Code
calculo_IMC <- function(nombre, peso, estatura) {
  
  return(list("nombre" = nombre, "IMC" = peso/(estatura^2)))
}

📝 Repite el ejercicio anterior pero con otro argumento opcional que se llame unidades (por defecto, unidades = "metros"). Desarrolla la función de forma que haga lo correcto si unidades = "metros" y si unidades = "centímetros".

Code
calculo_IMC <- function(nombre, peso, estatura, unidades = "metros") {
  
  return(list("nombre" = nombre,
              "IMC" = peso/(if_else(unidades == "metros", estatura, estatura/100)^2)))
}

📝 Crea un tibble ficticio de 7 personas, con tres variables (inventa nombre, y simula peso, estatura en centímetros), y aplica la función definida de forma que obtengamos una cuarta columna con su IMC.

Code
datos <-
  tibble("nombres" = c("javi", "sandra", "laura",
                       "ana", "carlos", "leo", NA),
         "peso" = rnorm(n = 7, mean = 70, sd = 1),
         "estatura" = rnorm(n = 7, mean = 168, sd = 5))

datos |> 
  mutate(IMC = calculo_IMC(nombres, peso, estatura, unidades = "centímetros")$IMC)

📝 Crea una función llamada atajo que tenga dos argumentos numéricos x e y. Si ambos son iguales, debes devolver "iguales" y hacer que la función acaba automáticamente (piensa cuándo una función sale). OJO: x e y podrían ser vectores. Si son distintos (de igual de longitud) calcula la proporción de elementos diferentes. Si son distintos (por ser distinta longitud), devuelve los elementos que no sean comunes.

Code
atajo <- function(x, y) {
  
  if (all(x == y) & length(x) == length(y)) { return("iguales") }
  else {
   
    if (length(x) == length(y)) {
      
      n_diff <- sum(x != y) / length(x)
      return(n_diff)
      
    } else {
      
      diff_elem <- unique(c(setdiff(x, y), setdiff(y, x)))
      return(diff_elem)
    }
    
  }
}

3.13 R base vs Tidyverse

Hasta ahora todo lo que hemos repasado en R lo hemos realizado en el paradigma de programación conocido como R base. Y es que cuando R nació como lenguaje, muchos de los que programaban en él imitaron formas y metodologías heredadas de otros lenguajes, basado en el uso de

  • Bucles for

  • Bucles while

  • Estructuras if-else

Y aunque conocer dichas estructuras puede sernos en algunos casos interesantes, en la mayoría de ocasiones han quedado caducas y vamos a poder evitarlas (en especial los bucles) ya que R está especialmente diseñado para trabajar de manera funcional (en lugar de elemento a elemento).

En ese contexto de programación funcional, hace una década nacía {tidyverse}, un «universo» de paquetes para garantizar un flujo de trabajo eficiente, coherente y lexicográficamente sencillo de entender, basado en la idea de que nuestros datos están limpios y ordenados (tidy)

3.13.1 Filosofía base: tidy data

Tidy datasets are all alike, but every messy dataset is messy in its own way (Hadley Wickham, Chief Scientist en RStudio)

TIDYVERSE

El universo de paquetes {tidyverse} se basa en la idea introducida por Hadley Wickham (el Dios al que rezamos) de estandarizar el formato de los datos para

  • sistematizar la depuración
  • hacer más sencillo su manipulación.
  • código legible

Lo primero por tanto será entender qué son los conjuntos tidydata ya que todo {tidyverse} se basa en que los datos están estandarizados.

  1. Cada variable en una única columna

  2. Cada individuo en una fila diferente

  3. Cada celda con un único valor

  4. Cada dataset en un tibble

  5. Si queremos cruzar múltiples tablas debemos tener una columna común

3.13.2 Tubería (pipe)

En {tidyverse} será clave el operador pipe (tubería) definido como |> (ctrl+shift+M): será una tubería que recorre los datos y los transforma.

En R base, si queremos aplicar tres funciones first(), second() y third() en orden, sería

third(second(first(datos)))

En {tidyverse} podremos leer de izquierda a derecha y separar los datos de las acciones

datos |> first() |> second() |> third()
Apunte importante

Desde la versión 4.1.0 de R disponemos de |>, un pipe nativo disponible fuera de tidyverse, sustituyendo al antiguo pipe %>% que dependía del paquete {magrittr} (bastante problemático).

La principal ventaja es que el código sea muy legible (casi literal) pudiendo hacer grandes operaciones con los datos con apenas código.

datos |>
  limpio(...) |>
  filtro(...) |>
  selecciono(...) |>
  ordeno(...) |>
  modifico(...) |>
  renombro(...) |>
  agrupo(...) |>
  cuento(...) |>
  resumo(...) |>
  pinto(...)

3.13.3 Datos SUCIOS: messy data

¿Pero qué aspecto tienen los datos no tidy? Vamos a cargar la tabla table4a del paquete {tidyr} (ya lo tenemos cargado del entorno tidyverse).

library(tidyr)
table4a
# A tibble: 3 × 3
  country     `1999` `2000`
  <chr>        <dbl>  <dbl>
1 Afghanistan    745   2666
2 Brazil       37737  80488
3 China       212258 213766

¿Qué puede estar fallando?

table4a
# A tibble: 3 × 3
  country     `1999` `2000`
  <chr>        <dbl>  <dbl>
1 Afghanistan    745   2666
2 Brazil       37737  80488
3 China       212258 213766

❎ Cada fila representa dos observaciones (1999 y 2000) → las columnas 1999 y 2000 en realidad deberían ser en sí valores de una variable y no nombres de columnas.

Incluiremos una nueva columna que nos guarde el año y otra que guarde el valor de la variable de interés en cada uno de esos años. Y lo haremos con la función pivot_longer(): pivotaremos la tabla a formato long:

table4a |> 
  pivot_longer(cols = c("1999", "2000"), names_to = "year", values_to = "cases")
# A tibble: 6 × 3
  country     year   cases
  <chr>       <chr>  <dbl>
1 Afghanistan 1999     745
2 Afghanistan 2000    2666
3 Brazil      1999   37737
4 Brazil      2000   80488
5 China       1999  212258
6 China       2000  213766
table4a |> 
  pivot_longer(cols = c("1999", "2000"),
               names_to = "year",
               values_to = "cases")
# A tibble: 6 × 3
  country     year   cases
  <chr>       <chr>  <dbl>
1 Afghanistan 1999     745
2 Afghanistan 2000    2666
3 Brazil      1999   37737
4 Brazil      2000   80488
5 China       1999  212258
6 China       2000  213766

  • cols: nombre de las variables a pivotar
  • names_to: nombre de la nueva variable a la quemandamos la cabecera de la tabla (los nombres).
  • values_to: nombre de la nueva variable a la que vamos a mandar los datos.

Veamos otro ejemplo con la tabla table2

table2
# A tibble: 12 × 4
   country      year type            count
   <chr>       <dbl> <chr>           <dbl>
 1 Afghanistan  1999 cases             745
 2 Afghanistan  1999 population   19987071
 3 Afghanistan  2000 cases            2666
 4 Afghanistan  2000 population   20595360
 5 Brazil       1999 cases           37737
 6 Brazil       1999 population  172006362
 7 Brazil       2000 cases           80488
 8 Brazil       2000 population  174504898
 9 China        1999 cases          212258
10 China        1999 population 1272915272
11 China        2000 cases          213766
12 China        2000 population 1280428583

¿Qué puede estar fallando?

# A tibble: 12 × 4
   country      year type            count
   <chr>       <dbl> <chr>           <dbl>
 1 Afghanistan  1999 cases             745
 2 Afghanistan  1999 population   19987071
 3 Afghanistan  2000 cases            2666
 4 Afghanistan  2000 population   20595360
 5 Brazil       1999 cases           37737
 6 Brazil       1999 population  172006362
 7 Brazil       2000 cases           80488
 8 Brazil       2000 population  174504898
 9 China        1999 cases          212258
10 China        1999 population 1272915272
11 China        2000 cases          213766
12 China        2000 population 1280428583

❎ Cada observación está dividido en dos filas → los registros con el mismo año deberían ser el mismo

Lo que haremos será lo opuesto: con pivot_wider() ensancharemos la tabla

table2 |>  pivot_wider(names_from = type, values_from = count)
# A tibble: 6 × 4
  country      year  cases population
  <chr>       <dbl>  <dbl>      <dbl>
1 Afghanistan  1999    745   19987071
2 Afghanistan  2000   2666   20595360
3 Brazil       1999  37737  172006362
4 Brazil       2000  80488  174504898
5 China        1999 212258 1272915272
6 China        2000 213766 1280428583

Veamos otro ejemplo con la tabla table3

table3
# A tibble: 6 × 3
  country      year rate             
  <chr>       <dbl> <chr>            
1 Afghanistan  1999 745/19987071     
2 Afghanistan  2000 2666/20595360    
3 Brazil       1999 37737/172006362  
4 Brazil       2000 80488/174504898  
5 China        1999 212258/1272915272
6 China        2000 213766/1280428583

¿Qué puede estar fallando?

table3
# A tibble: 6 × 3
  country      year rate             
  <chr>       <dbl> <chr>            
1 Afghanistan  1999 745/19987071     
2 Afghanistan  2000 2666/20595360    
3 Brazil       1999 37737/172006362  
4 Brazil       2000 80488/174504898  
5 China        1999 212258/1272915272
6 China        2000 213766/1280428583

❎ Cada celda contiene varios valores

Lo que haremos será hacer uso de la función separate() para mandar separar cada valor a una columna diferente.

table3 |> separate(rate, into = c("cases", "pop"))
# A tibble: 6 × 4
  country      year cases  pop       
  <chr>       <dbl> <chr>  <chr>     
1 Afghanistan  1999 745    19987071  
2 Afghanistan  2000 2666   20595360  
3 Brazil       1999 37737  172006362 
4 Brazil       2000 80488  174504898 
5 China        1999 212258 1272915272
6 China        2000 213766 1280428583
table3 |> separate(rate, into = c("cases", "pop"))
# A tibble: 6 × 4
  country      year cases  pop       
  <chr>       <dbl> <chr>  <chr>     
1 Afghanistan  1999 745    19987071  
2 Afghanistan  2000 2666   20595360  
3 Brazil       1999 37737  172006362 
4 Brazil       2000 80488  174504898 
5 China        1999 212258 1272915272
6 China        2000 213766 1280428583

Fíjate que los datos, aunque los ha separado, los ha mantenido como texto cuando en realidad deberían ser variables numéricas. Para ello podemos añadir el argumento opcional convert = TRUE

table3 |> separate(rate, into = c("cases", "pop"), convert = TRUE)
# A tibble: 6 × 4
  country      year  cases        pop
  <chr>       <dbl>  <int>      <int>
1 Afghanistan  1999    745   19987071
2 Afghanistan  2000   2666   20595360
3 Brazil       1999  37737  172006362
4 Brazil       2000  80488  174504898
5 China        1999 212258 1272915272
6 China        2000 213766 1280428583

Veamos el último ejemplo con la tabla table5

table5
# A tibble: 6 × 4
  country     century year  rate             
  <chr>       <chr>   <chr> <chr>            
1 Afghanistan 19      99    745/19987071     
2 Afghanistan 20      00    2666/20595360    
3 Brazil      19      99    37737/172006362  
4 Brazil      20      00    80488/174504898  
5 China       19      99    212258/1272915272
6 China       20      00    213766/1280428583

¿Qué puede estar fallando?

table5
# A tibble: 6 × 4
  country     century year  rate             
  <chr>       <chr>   <chr> <chr>            
1 Afghanistan 19      99    745/19987071     
2 Afghanistan 20      00    2666/20595360    
3 Brazil      19      99    37737/172006362  
4 Brazil      20      00    80488/174504898  
5 China       19      99    212258/1272915272
6 China       20      00    213766/1280428583

❎ Tenemos mismos valores divididos en dos columnas

Usaremos unite() para unir los valores de siglo y año en una misma columna

table5 |> unite(col = year_completo, century, year, sep = "")
# A tibble: 6 × 3
  country     year_completo rate             
  <chr>       <chr>         <chr>            
1 Afghanistan 1999          745/19987071     
2 Afghanistan 2000          2666/20595360    
3 Brazil      1999          37737/172006362  
4 Brazil      2000          80488/174504898  
5 China       1999          212258/1272915272
6 China       2000          213766/1280428583

3.13.4 Ejemplo: relig_income

Vamos a realizar un ejemplo juntos con la tabla relig_income del paquete {tidyr}. Como se indica en la ayuda ? relig_income, la tabla representa la cantidad de personas que hay en cada tramo de ingresos anuales (20k = 20 000$) y en cada religión.

relig_income
# A tibble: 18 × 11
   religion `<$10k` `$10-20k` `$20-30k` `$30-40k` `$40-50k` `$50-75k` `$75-100k`
   <chr>      <dbl>     <dbl>     <dbl>     <dbl>     <dbl>     <dbl>      <dbl>
 1 Agnostic      27        34        60        81        76       137        122
 2 Atheist       12        27        37        52        35        70         73
 3 Buddhist      27        21        30        34        33        58         62
 4 Catholic     418       617       732       670       638      1116        949
 5 Don’t k…      15        14        15        11        10        35         21
 6 Evangel…     575       869      1064       982       881      1486        949
 7 Hindu          1         9         7         9        11        34         47
 8 Histori…     228       244       236       238       197       223        131
 9 Jehovah…      20        27        24        24        21        30         15
10 Jewish        19        19        25        25        30        95         69
11 Mainlin…     289       495       619       655       651      1107        939
12 Mormon        29        40        48        51        56       112         85
13 Muslim         6         7         9        10         9        23         16
14 Orthodox      13        17        23        32        32        47         38
15 Other C…       9         7        11        13        13        14         18
16 Other F…      20        33        40        46        49        63         46
17 Other W…       5         2         3         4         2         7          3
18 Unaffil…     217       299       374       365       341       528        407
# ℹ 3 more variables: `$100-150k` <dbl>, `>150k` <dbl>,
#   `Don't know/refused` <dbl>

¿Es tidydata?

No lo es ya que en realidad solo deberíamos tener una variable de ingresos y la tenemos dividida en 11: todas ellas es la misma variable solo que adopta un valor diferente. ¿Cómo convertirla a tidy data?

La idea es pivotar todas las columnas de ingresos para que acaben en una sola columna llamada income, y los valores (el número de personas) en otra llamada people (por ejemplo). La tabla la haremos más larga y menos ancha así que…

relig_tidy <-
  relig_income |>
  pivot_longer(cols = "<$10k":"Don't know/refused", names_to = "income",
               values_to = "people")
relig_tidy 
# A tibble: 180 × 3
   religion income             people
   <chr>    <chr>               <dbl>
 1 Agnostic <$10k                  27
 2 Agnostic $10-20k                34
 3 Agnostic $20-30k                60
 4 Agnostic $30-40k                81
 5 Agnostic $40-50k                76
 6 Agnostic $50-75k               137
 7 Agnostic $75-100k              122
 8 Agnostic $100-150k             109
 9 Agnostic >150k                  84
10 Agnostic Don't know/refused     96
# ℹ 170 more rows

Vamos a hilar más fino: ahora mismo en la variable income en realidad tenemos dos valores, el límite inferior y el superior de la renta. Vamos a separar dicha variable e ingresos en dos, llamadas income_inf y income_sup

relig_tidy 
# A tibble: 180 × 3
   religion income             people
   <chr>    <chr>               <dbl>
 1 Agnostic <$10k                  27
 2 Agnostic $10-20k                34
 3 Agnostic $20-30k                60
 4 Agnostic $30-40k                81
 5 Agnostic $40-50k                76
 6 Agnostic $50-75k               137
 7 Agnostic $75-100k              122
 8 Agnostic $100-150k             109
 9 Agnostic >150k                  84
10 Agnostic Don't know/refused     96
# ℹ 170 more rows

Vamos a hilar más fino: ahora mismo en la variable income en realidad tenemos dos valores, el límite inferior y el superior de la renta. Vamos a separar dicha variable e ingresos en dos, llamadas income_inf y income_sup

relig_tidy |>
  # Separamos por -
  separate(income, into = c("income_inf", "income_sup"), sep = "-")
Warning: Expected 2 pieces. Missing pieces filled with `NA` in 54 rows [1, 9, 10, 11,
19, 20, 21, 29, 30, 31, 39, 40, 41, 49, 50, 51, 59, 60, 61, 69, ...].
# A tibble: 180 × 4
   religion income_inf         income_sup people
   <chr>    <chr>              <chr>       <dbl>
 1 Agnostic <$10k              <NA>           27
 2 Agnostic $10                20k            34
 3 Agnostic $20                30k            60
 4 Agnostic $30                40k            81
 5 Agnostic $40                50k            76
 6 Agnostic $50                75k           137
 7 Agnostic $75                100k          122
 8 Agnostic $100               150k          109
 9 Agnostic >150k              <NA>           84
10 Agnostic Don't know/refused <NA>           96
# ℹ 170 more rows

¿Está ya bien? Fíjate bien…

relig_tidy |>
  # Separamos por -
  separate(income, into = c("income_inf", "income_sup"), sep = "-")
Warning: Expected 2 pieces. Missing pieces filled with `NA` in 54 rows [1, 9, 10, 11,
19, 20, 21, 29, 30, 31, 39, 40, 41, 49, 50, 51, 59, 60, 61, 69, ...].
# A tibble: 180 × 4
   religion income_inf         income_sup people
   <chr>    <chr>              <chr>       <dbl>
 1 Agnostic <$10k              <NA>           27
 2 Agnostic $10                20k            34
 3 Agnostic $20                30k            60
 4 Agnostic $30                40k            81
 5 Agnostic $40                50k            76
 6 Agnostic $50                75k           137
 7 Agnostic $75                100k          122
 8 Agnostic $100               150k          109
 9 Agnostic >150k              <NA>           84
10 Agnostic Don't know/refused <NA>           96
# ℹ 170 more rows

Si te fijas la primera columna el "$10k" debería ser una cota superior, no inferior. ¿Cómo indicarle que separe bien ese caso? Le indicaremos que separe si encuentra "-" o "<" (usamos | para separar ambas opciones)

relig_tidy |>
  # Separamos por -
  separate(income, into = c("income_inf", "income_sup"), sep = "-|<")
Warning: Expected 2 pieces. Missing pieces filled with `NA` in 36 rows [9, 10, 19, 20,
29, 30, 39, 40, 49, 50, 59, 60, 69, 70, 79, 80, 89, 90, 99, 100, ...].
# A tibble: 180 × 4
   religion income_inf           income_sup people
   <chr>    <chr>                <chr>       <dbl>
 1 Agnostic ""                   $10k           27
 2 Agnostic "$10"                20k            34
 3 Agnostic "$20"                30k            60
 4 Agnostic "$30"                40k            81
 5 Agnostic "$40"                50k            76
 6 Agnostic "$50"                75k           137
 7 Agnostic "$75"                100k          122
 8 Agnostic "$100"               150k          109
 9 Agnostic ">150k"              <NA>           84
10 Agnostic "Don't know/refused" <NA>           96
# ℹ 170 more rows
relig_tidy <-
  relig_tidy |>
  # Separamos por -
  separate(income, into = c("income_inf", "income_sup"), sep = "-|<")
Warning: Expected 2 pieces. Missing pieces filled with `NA` in 36 rows [9, 10, 19, 20,
29, 30, 39, 40, 49, 50, 59, 60, 69, 70, 79, 80, 89, 90, 99, 100, ...].
relig_tidy
# A tibble: 180 × 4
   religion income_inf           income_sup people
   <chr>    <chr>                <chr>       <dbl>
 1 Agnostic ""                   $10k           27
 2 Agnostic "$10"                20k            34
 3 Agnostic "$20"                30k            60
 4 Agnostic "$30"                40k            81
 5 Agnostic "$40"                50k            76
 6 Agnostic "$50"                75k           137
 7 Agnostic "$75"                100k          122
 8 Agnostic "$100"               150k          109
 9 Agnostic ">150k"              <NA>           84
10 Agnostic "Don't know/refused" <NA>           96
# ℹ 170 more rows

Piensa ahora como podemos convertir los límites de ingresos a numéricas (eliminando símbolos, letras, etc)

Para ello usaremos el paquete {stringr}, en concreto la función str_remove_all(), a la que le podemos pasar los caracteres que queremos eliminar (fíjate que $ al ser un caracter reservado en R hay que indicárselo con \\$)

relig_tidy$income_inf <-
  str_remove_all(relig_tidy$income_inf, "\\$|>|k")
relig_tidy$income_sup <-
  str_remove_all(relig_tidy$income_sup, "\\$|>|k")

relig_tidy
# A tibble: 180 × 4
   religion income_inf          income_sup people
   <chr>    <chr>               <chr>       <dbl>
 1 Agnostic ""                  10             27
 2 Agnostic "10"                20             34
 3 Agnostic "20"                30             60
 4 Agnostic "30"                40             81
 5 Agnostic "40"                50             76
 6 Agnostic "50"                75            137
 7 Agnostic "75"                100           122
 8 Agnostic "100"               150           109
 9 Agnostic "150"               <NA>           84
10 Agnostic "Don't now/refused" <NA>           96
# ℹ 170 more rows

Fíjate que tenemos "Don't now/refused". ¿Qué deberíamos tener?

Debería ser un dato ausente así que usaremos if_else(): si contiene dicha frase, NA, en caso contrario su valor (consejo: str_detect() para detectar patrones en textos, y evitar tener que escribir toda la palabra sin errores)

relig_tidy$income_inf <-
  if_else(str_detect(relig_tidy$income_inf, "refused"), NA, relig_tidy$income_inf)
relig_tidy$income_sup <-
  if_else(str_detect(relig_tidy$income_sup, "refused"), NA, relig_tidy$income_sup)
relig_tidy
# A tibble: 180 × 4
   religion income_inf income_sup people
   <chr>    <chr>      <chr>       <dbl>
 1 Agnostic ""         10             27
 2 Agnostic "10"       20             34
 3 Agnostic "20"       30             60
 4 Agnostic "30"       40             81
 5 Agnostic "40"       50             76
 6 Agnostic "50"       75            137
 7 Agnostic "75"       100           122
 8 Agnostic "100"      150           109
 9 Agnostic "150"      <NA>           84
10 Agnostic  <NA>      <NA>           96
# ℹ 170 more rows

En la primera línea, ese "" también debería ser `NA``

relig_tidy$income_inf <-
  if_else(relig_tidy$income_inf == "", NA, relig_tidy$income_inf)
relig_tidy$income_suop <-
  if_else(relig_tidy$income_sup == "", NA, relig_tidy$income_sup)

Además si te fijas los números son en realidad caracteres, así que vamos a convertirlos a números

relig_tidy$income_inf <- as.numeric(relig_tidy$income_inf)
relig_tidy$income_sup <- as.numeric(relig_tidy$income_sup)
relig_tidy
# A tibble: 180 × 5
   religion income_inf income_sup people income_suop
   <chr>         <dbl>      <dbl>  <dbl> <chr>      
 1 Agnostic         NA         10     27 10         
 2 Agnostic         10         20     34 20         
 3 Agnostic         20         30     60 30         
 4 Agnostic         30         40     81 40         
 5 Agnostic         40         50     76 50         
 6 Agnostic         50         75    137 75         
 7 Agnostic         75        100    122 100        
 8 Agnostic        100        150    109 150        
 9 Agnostic        150         NA     84 <NA>       
10 Agnostic         NA         NA     96 <NA>       
# ℹ 170 more rows

¿Se te ocurre alguna forma de «cuantificar numéricamente» los valores ausentes que tenemos en este caso? Si te fijas en realidad cuando hay ausente en el límite inferior en realidad podríamos poner un 0 (nadie puede ganar menos de eso) y cuando lo tenemos en el límite superior sería Inf

relig_tidy$income_inf <-
  if_else(is.na(relig_tidy$income_inf), 0, relig_tidy$income_inf)
relig_tidy$income_sup <-
  if_else(is.na(relig_tidy$income_sup), Inf, relig_tidy$income_sup)
relig_tidy
# A tibble: 180 × 5
   religion income_inf income_sup people income_suop
   <chr>         <dbl>      <dbl>  <dbl> <chr>      
 1 Agnostic          0         10     27 10         
 2 Agnostic         10         20     34 20         
 3 Agnostic         20         30     60 30         
 4 Agnostic         30         40     81 40         
 5 Agnostic         40         50     76 50         
 6 Agnostic         50         75    137 75         
 7 Agnostic         75        100    122 100        
 8 Agnostic        100        150    109 150        
 9 Agnostic        150        Inf     84 <NA>       
10 Agnostic          0        Inf     96 <NA>       
# ℹ 170 more rows

Aunque nos haya llevado un rato este es el código completo resumido

relig_tidy <-
  relig_income |>
  pivot_longer(cols = "<$10k":"Don't know/refused", names_to = "income",
               values_to = "people") |>
  separate(income, into = c("income_inf", "income_sup"), sep = "-|<")

relig_tidy$income_inf <- str_remove_all(relig_tidy$income_inf, "\\$|>|k")
relig_tidy$income_sup <- str_remove_all(relig_tidy$income_sup, "\\$|>|k")

relig_tidy$income_inf <-
  if_else(str_detect(relig_tidy$income_inf, "refused") |
            relig_tidy$income_inf == "", 0, as.numeric(relig_tidy$income_inf))
relig_tidy$income_sup <-
  if_else(str_detect(relig_tidy$income_sup, "refused") |
            relig_tidy$income_sup == "", Inf, as.numeric(relig_tidy$income_sup))

¿Por qué era importante tenerlo en tidydata? Lo veremos más adelante al visualizar los datos pero esto ya nos permite realizar filtros muy rápidos con muy poco código.

Por ejemplo: ¿cuántas personas agnósticas con ingresos superiores (o iguales) a 30 tenemos?

# una línea de código
sum(relig_tidy$people[relig_tidy$religion == "Agnostic" & relig_tidy$income_inf >= 30])
[1] 609

3.14 💻 Tu turno

Intenta realizar los siguientes ejercicios sin mirar las soluciones

📝 Usa el dataset original relig_income y trata de responder a la última pregunta: ¿cuántas personas agnósticas con ingresos superiores (o iguales) a 30 tenemos? Compara el código a realizar cuando tenemos tidydata a cuando no. ¿Cuál es más legible si no supieses R? ¿Cuál tiene mayor probabilidad de error?

Code
sum(relig_income[relig_income$religion == "Agnostic",
             c("$30-40k", "$40-50k", "$50-75k", "$75-100k", "$100-150k", ">150k")])

📝 Usando relig_tidy determina quién tiene más ingresos medios, ¿católicos (Catholic) o agnósticos (Agnostic)? Crea antes una variable avg_income (ingresos medios por intervalo): si hay 5 personas entre 20 y 30, y 3 personas entre 30 y 50, la media sería \((25*5 + 40*3)/8\) (si es Inf por arriba, NA)

Code
relig_tidy$avg_income <- 
  if_else(is.infinite(relig_tidy$income_sup), NA, (relig_tidy$income_sup + relig_tidy$income_inf)/2)

# Agnosticos vs catolicos
sum((relig_tidy$avg_income[relig_tidy$religion == "Agnostic"] * relig_tidy$people[relig_tidy$religion == "Agnostic"]), na.rm = TRUE) /
  sum(relig_tidy$people[relig_tidy$religion == "Agnostic"], na.rm = TRUE)

sum((relig_tidy$avg_income[relig_tidy$religion == "Catholic"] * relig_tidy$people[relig_tidy$religion == "Catholic"]), na.rm = TRUE) /
  sum(relig_tidy$people[relig_tidy$religion == "Catholic"], na.rm = TRUE)

📝 Si debemos elegir budismo (Buddhist) e hinduismo (Hindu), ¿cuál de las dos es la religión mayoritaria entre los que ganan más de 50 000$ anuales?

Code
greatest_income <-
  relig_tidy[relig_tidy$income_inf >= 50 & relig_tidy$religion %in% c("Buddhist", "Hindu"), ]

sum(greatest_income$people[greatest_income$religion == "Buddhist"], na.rm = TRUE)
sum(greatest_income$people[greatest_income$religion == "Hindu"], na.rm = TRUE)

📝 Echa un vistazo a la tabla table4b del paquete {tidyr}. ¿Es tidydata? En caso negativo, ¿qué falla? ¿Cómo convertirla a tidy data en caso de que no lo sea ya?

Code
table4b |>
  pivot_longer(cols = "1999":"2000", names_to = "year",
               values_to = "cases")

📝 Echa un vistazo a la tabla billboard del paquete {tidyr}. ¿Es tidydata? En caso negativo, ¿qué falla? ¿Cómo convertirla a tidy data en caso de que no lo sea ya?

Code
billboard |>
  pivot_longer(cols = "wk1":"wk76",
               names_to = "week",
               names_prefix = "wk",
               values_to = "position",
               values_drop_na = TRUE)

4 🐣 Caso práctico I: tibble

Del paquete {Biostatistics} usaremos el conunto de datos pinniped, que guarda los datos de peso de cuerpo y cerebro (desagregado por sexo y mono/poligamia) de 33 especies de mamíferos marinos.

Biostatistics::pinniped
                       Species Male_brain_g Female_brain_g Male_mass_Kg
1       Monachus schauinslandi        370.0             NA        173.0
2            Monachus monachus        480.0          480.0        260.0
3      Mirounga angustirostris        700.0          640.0       2275.0
4             Mirounga leonina       1431.3          898.8       3510.0
5       Leptonychotes weddelli        535.0          637.5        450.0
6            Ommatophoca rossi        425.0          530.0        153.8
7        Lobodon carcinophagus        578.2          538.8        220.5
8            Hydrurga leptonyx        765.0          660.0        324.0
9          Cystophora cristata        480.0          430.0        343.2
10         Erignathus barbatus           NA          460.0        312.5
11          Halichoerus grypus        342.5          272.5        233.0
12          Phoca groenlandica        297.5          252.5        145.0
13              Phoca fasciata        257.5          240.0         94.8
14                Phoca largha        257.5          250.0         97.0
15               Phoca caspica        165.0          160.0         70.5
16              Phoca sibirica        185.0          190.0         89.5
17               Phoca hispida        229.3          220.0         84.0
18              Phoca vitulina        362.3          265.0         97.1
19      Zalophus californianus        405.0          361.5        244.5
20          Eumetopias jubatus        747.5          575.0       1000.0
21              Otaria byronia        546.3          470.0        300.0
22            Neophoca cinerea        440.0          337.5        300.0
23          Phocarctos hookeri        417.5          370.0        364.0
24         Callorhinus ursinus        355.0          302.5        140.0
25     Arctocephalus townsendi           NA             NA        112.0
26     Arctocephalus philippii        415.0             NA        140.0
27 Arctocephalus galapagoensis        302.5          280.0         64.5
28     Arctocephalus australis        350.0          265.0         91.0
29      Arctocephalus forsteri        340.0          300.0        125.0
30       Arctocephalus gazella        360.0          320.0        155.0
31    Arctocephalus tropicalis        322.5          330.0        152.5
32      Arctocephalus pusillus        401.3          337.5        263.0
33           Odobenus rosmarus       1303.0         1340.5       1233.0
   Female_mass_Kg Mate_type
1           272.2      mono
2           275.0      mono
3           488.0      poly
4           565.7      poly
5           447.0      poly
6           164.0      mono
7           224.0      mono
8           367.0      mono
9           222.5      mono
10          326.0      mono
11          205.8      poly
12          139.0      mono
13           80.4      mono
14           71.3      mono
15           55.0      mono
16           85.0      mono
17           81.2      mono
18           85.2      mono
19           81.0      poly
20          287.5      poly
21          144.0      poly
22           78.6      poly
23          114.7      poly
24           33.3      poly
25           49.6      poly
26           48.1      poly
27           27.4      poly
28           48.5      poly
29           38.1      poly
30           45.0      poly
31           50.0      poly
32           64.1      poly
33          811.5      poly

Del paquete {Biostatistics} usaremos el conunto de datos pinniped, que guarda los datos de peso de cuerpo y cerebro (desagregado por sexo y mono/poligamia) de 33 especies de mamíferos marinos.

Biostatistics::pinniped
                       Species Male_brain_g Female_brain_g Male_mass_Kg
1       Monachus schauinslandi        370.0             NA        173.0
2            Monachus monachus        480.0          480.0        260.0
3      Mirounga angustirostris        700.0          640.0       2275.0
4             Mirounga leonina       1431.3          898.8       3510.0
5       Leptonychotes weddelli        535.0          637.5        450.0
6            Ommatophoca rossi        425.0          530.0        153.8
7        Lobodon carcinophagus        578.2          538.8        220.5
8            Hydrurga leptonyx        765.0          660.0        324.0
9          Cystophora cristata        480.0          430.0        343.2
10         Erignathus barbatus           NA          460.0        312.5
11          Halichoerus grypus        342.5          272.5        233.0
12          Phoca groenlandica        297.5          252.5        145.0
13              Phoca fasciata        257.5          240.0         94.8
14                Phoca largha        257.5          250.0         97.0
15               Phoca caspica        165.0          160.0         70.5
16              Phoca sibirica        185.0          190.0         89.5
17               Phoca hispida        229.3          220.0         84.0
18              Phoca vitulina        362.3          265.0         97.1
19      Zalophus californianus        405.0          361.5        244.5
20          Eumetopias jubatus        747.5          575.0       1000.0
21              Otaria byronia        546.3          470.0        300.0
22            Neophoca cinerea        440.0          337.5        300.0
23          Phocarctos hookeri        417.5          370.0        364.0
24         Callorhinus ursinus        355.0          302.5        140.0
25     Arctocephalus townsendi           NA             NA        112.0
26     Arctocephalus philippii        415.0             NA        140.0
27 Arctocephalus galapagoensis        302.5          280.0         64.5
28     Arctocephalus australis        350.0          265.0         91.0
29      Arctocephalus forsteri        340.0          300.0        125.0
30       Arctocephalus gazella        360.0          320.0        155.0
31    Arctocephalus tropicalis        322.5          330.0        152.5
32      Arctocephalus pusillus        401.3          337.5        263.0
33           Odobenus rosmarus       1303.0         1340.5       1233.0
   Female_mass_Kg Mate_type
1           272.2      mono
2           275.0      mono
3           488.0      poly
4           565.7      poly
5           447.0      poly
6           164.0      mono
7           224.0      mono
8           367.0      mono
9           222.5      mono
10          326.0      mono
11          205.8      poly
12          139.0      mono
13           80.4      mono
14           71.3      mono
15           55.0      mono
16           85.0      mono
17           81.2      mono
18           85.2      mono
19           81.0      poly
20          287.5      poly
21          144.0      poly
22           78.6      poly
23          114.7      poly
24           33.3      poly
25           49.6      poly
26           48.1      poly
27           27.4      poly
28           48.5      poly
29           38.1      poly
30           45.0      poly
31           50.0      poly
32           64.1      poly
33          811.5      poly

4.1 Pregunta 1

Comprueba si los datos están en formato tibble. En caso negativo conviértelo.

Code
# chequeamos si es tibble
library(tibble)
is_tibble(Biostatistics::pinniped)

# Convertimos a tibble
pinniped_tb <- as_tibble(Biostatistics::pinniped)

4.2 Pregunta 2

¿Cuántos registros hay? ¿Cuántas variables? ¿De qué tipo es cada una? ¿Cuáles son sus nombres?

Code
nrow(pinniped_tb)
ncol(pinniped_tb)
names(pinniped_tb)

4.3 Pregunta 3

Incorpora una variable nueva llamada phoca que sea de tipo lógico y que nos diga si una especie es de la categoría Phoca o no.

Code
pinniped_tb$phoca <- pinniped_tb$Species == "Phoca"

4.4 Pregunta 4

¿A qué sexo le pesa más el cerebro: a las hembras o a los machos?

Code
# ¿a quién le pesa más el cerebro?
mean(pinniped_tb$Male_brain_g, na.rm = TRUE) >
  mean(pinniped_tb$Female_brain_g, na.rm = TRUE)

4.5 Pregunta 5

¿A quienes les pesa más el cuerpo a los monógamos o a los polígamos? Recuerda que tienes los pesos divididos por sexos en variables distintas que tendrás que juntar de alguna forma

Code
# ¿a quién le pesa más el cerebro?
mean(c(pinniped_tb$Male_mass_Kg[pinniped_tb$Mate_type == "mono"],
       pinniped_tb$Female_mass_Kg[pinniped_tb$Mate_type == "mono"])) >
  mean(c(pinniped_tb$Male_mass_Kg[pinniped_tb$Mate_type == "poly"],
         pinniped_tb$Female_mass_Kg[pinniped_tb$Mate_type == "poly"]))

4.6 Pregunta 6

Incopora una nueva variable llamada dif_m_f que represente la diferencia entre el peso del cerebro entre machos y hembras (machos - hembras) para cada especie.

Code
pinniped_tb$dif_m_f <- pinniped_tb$Male_brain_g - pinniped_tb$Female_brain_g
pinniped_tb

5 🐣 Caso práctico II: bucles y condicionales

Para practicar estructuras de control vamos a realizar un ejercicio de simulación

5.1 Pregunta 1

Define una varible llamada importe que empiece en 100. Diseña un bucle de 20 iteraciones donde en cada iteración, importe se reduzca a la mitad de su valor. Piensa que tipo de estructura de bucle deberías usar. El valor final de importe deberia ser 0.000095367 (aprox)

Code
# Usamos un for ya que sabemos el número de iteraciones 
# de manera prefijada (y no depende de nada)

# definimos inicialmente importe en 100
importe <- 100 

# para el bucle usamos por ejemplo i como índice, que va de 1 a 20
for (i in 1:20) {
  
  # el código fíjate que es el mismo y no depende de i
  importe <- importe/2
}
importe

5.2 Pregunta 2

Diseña una estructura de bucle de manera que encuentres la iteración en la que importe es menor que 0.001 por primera vez. Una vez encontrado guárdalo en iter y para el bucle.

Code
# dos formas de hacerlo: for y while

# con for
importe <- 100 

# ya sabemos que en 20 es menor que 0.001 así que podemos poner
# dicha cantidad como tope sabiendo que no llegará
for (i in 1:20) {
  
  # si todavía no es menor, seguimos dividiendo
  if (importe >= 0.001) {
    
    importe <- importe/2
    
  } else {
    
    # si ya es menor, guardamos la iteración (piensa por qué i - 1)
    iter <- i - 1 
    
    # y paramos
    break
  }
  
}

# con while
importe <- 100 

iter <- 0 # debemos inicializar las iteraciones

# no sabemos cuantas iteraciones, solo que debe parar cuando
# importe esté por debajo de dicha cantidad
while (importe >= 0.001) {
  
  importe <- importe/2
    
  # estructura clásica de while: si se corre iteración
  # actualizamos un valor (en este caso que cuente una iteración)
  iter <- iter + 1
}

iter

5.3 Pregunta 3

En R tenemos la función %%: si ponemos a %% b nos devuelve el resto que daría la división \(a/b\). Por ejemplo, 4 %% 2 da 0 ya que 4 es un número par (es decir, su resto al dividir entre 2 es 0). Si ponemos 13 %% 5 nos devuelve 3, ya que el resto de dividir 13 entre 5 es 3.

# Resto al dividir entre 2
3 %% 2
[1] 1
4 %% 2
[1] 0
5 %% 2
[1] 1
6 %% 2
[1] 0
# Resto al dividir entre 3
9 %% 3
[1] 0
10 %% 3
[1] 1
11 %% 3
[1] 2
12 %% 3
[1] 0

Empezando en un importe inicial importe_inicial de 100 (euros), diseña un bucle que te sume 3€ más la iteración por la que estés si el importe actual es par y te reste 5€ menos la iteración por la que estés si es impar, SALVO que el importe ya esté igual o por debajo de 0 (en ese caso no debe sumar ni retar). Ejemplo: si importe tiene 50 euros y estás en la iteración 13, sumará 3 + 13 (66 en total); si importe tiene 51 euros y estás en la iteración 13, restará 5 + 13 (33 en total); si importe tiene -2 euros y estás en la iteración 13, sumará 3 + 13 (14 en total); si importe tiene -1 euros y estás en la iteración 13, no hará nada. Guarda los importes resultantes de cada iteración (máximo de 150 iteraciones). Empieza a partir de la iteración 2

Code
importe_inicial <- 100
importe <- c(importe_inicial, rep(NA, 149))
for (i in 2:150) {
  
  if (importe[i - 1] %% 2 == 0) {
    
    importe[i] <- importe[i - 1] + 3 + i
    
  } else if (importe[i - 1] > 0) {
    
    importe[i] <- importe[i - 1] - (5 + i)
    
  } else {
    
    importe[i] <- importe[i - 1]
    
  }
}

¿Qué ha pasado?

5.4 Pregunta 4

En R la función sample(x = ..., size = ...) va sernos muy útil: de una colección de elementos x, selecciona un número size al azar de ellos.

Por ejemplo, si queremos simular 3 veces el lanzamiento de un dado tenemos 6 elementos posibles (x = 1:6) y lo seleccionamos 3 veces (size = 3)

sample(x = 1:6, size = 3)
[1] 6 3 2

Al ser aleatorio, cada vez que lo ejecutas saldrá algo distinto

sample(x = 1:6, size = 3)
[1] 1 3 2

¿Y si queremos tirar 10 veces?

sample(x = 1:6, size = 10)
Error in sample.int(length(x), size, replace, prob): cannot take a sample larger than the population when 'replace = FALSE'

Al tener solo 6 elementos posibles y elegir 10, no puede, así que le tenemos que indicar que queremos un sample (muestreo) con reemplazamiento (como sucede en el dado, cada cara puede repetirse al volver a tirarlo)

sample(x = 1:6, size = 10, replace = TRUE)
 [1] 5 3 4 3 6 6 3 5 1 1

Con lo anterior, imagina que estás en un concurso de televisión donde te dan a elegir 3 puertas: en una hay un premio millonarios y en las otras 2 una galleta oreo. Diseña el estudio de simulación con bucles for para aproximar la probabilidad de que te toque el premio (obviamente tiene que darte aprox 0.3333333). Realiza el experimento para 10, 50 intentos, 100 intentos, 500 intentos, 1000 intentos, 10 000 intentos y 25 000 intentos (pista: necesitas un bucle dentro de otro). ¿Qué observas?

Code
# Definimos las posibilidades
puertas <- c(1, 2, 3)

# Definimos los intentos
intentos <- c(10, 50, 100, 500, 1000, 10000, 25000)

# Para cada escenario de intentos, definimos las veces que hemos ganado
# (al inicio empieza en 0 claro)
exitos <- rep(0, length(intentos))

# primer bucle: cantidad de intentos permitidos
for (i in 1:length(intentos)) {
  
  # segundo bucle: para cada intento, simulaciones una eleccion de 
  # puerta y un premio
  for (j in 1:intentos[i]) {
    
    # premio: de 3 puertas, solo está en una
    premio <- sample(x = puertas, size = 1)
    
    # puerta que seleccionas como concursante: de 3 puertas, te quedas con una
    eleccion <- sample(x = puertas, size = 1)
    
    # si la puerta seleccionada coincide con la que tiene el premio
    # sumas un éxito, sino te quedas como estás
    exitos[i] <- if_else(eleccion == premio, exitos[i] + 1, exitos[i])
    
  }
  # Tras jugar, lo dividimos entre el número de veces que has jugado
  # para tener una proporción
  exitos[i] <- exitos[i] / intentos[i]
}
exitos

¿Y si en cada ronda, te abriesen una de las puertas no premiadas que no has elegido, cambiarías de puerta o te mantendrías? Simula ambos casos y descubre cuál es la estrategia correcta (este problema se conoce como problema de Monty Hall y aparece incluso en películas como 21 Black Jack)

Code
puertas <- c(1, 2, 3)
intentos <- c(10, 50, 100, 500, 1000, 10000, 25000)
exitos_mantengo <- exitos_cambio <- rep(0, length(intentos))

for (i in 1:length(intentos)) {
  for (j in 1:intentos[i]) {
    
    # puerta que seleccionas como concursante: de 3 puertas, te quedas con una
    eleccion_inicial <- sample(x = puertas, size = 1)
    
    # premio: de 3 puertas, solo está en una
    premio <- sample(x = puertas, size = 1)
    
    # De la no elegida, el presentador te abre una no premiada
    puerta_abierta <-
      puertas[puertas != eleccion_inicial & puertas != premio]
    
    # si solo hay una opción (es decir que tu eleccion inicial
    # y el premio son puertas distintas) no haces nada
    
    # Si hubiese 2 opciones (si tu eleccion y el premio coinciden)
    # te abrirá una al azar
    if (length(puerta_abierta) > 1) {
      
      puerta_abierta <- sample(x = puerta_abierta, size = 1)
    }
      
    
    # si mantienes es como antes
    exitos_mantengo[i] <-
      if_else(eleccion == premio, exitos_mantengo[i] + 1, exitos_mantengo[i])
    
    # si cambias es a una puerta distinta de la inicial y de la abierta, la que quede
    cambio <- puertas[puertas != eleccion_inicial & puertas != puerta_abierta]
    exitos_cambio[i] <-
      if_else(cambio == premio, exitos_cambio[i] + 1, exitos_cambio[i])
    
  }
  # Tras jugar, lo dividimos entre el número de veces que has jugado
  # para tener una proporción
  exitos_mantengo[i] <- exitos_mantengo[i] / intentos[i]
  exitos_cambio[i] <- exitos_cambio[i] / intentos[i]
}
exitos_mantengo
exitos_cambio

¿Qué sucede?

6 🐣 Caso práctico III: tidy data

En el paquete {tidyr} contamos con el dataset who2 (dataset de la Organización Mundial de la Salud - OMS)

library(tidyr)
who2
# A tibble: 7,240 × 58
   country      year sp_m_014 sp_m_1524 sp_m_2534 sp_m_3544 sp_m_4554 sp_m_5564
   <chr>       <dbl>    <dbl>     <dbl>     <dbl>     <dbl>     <dbl>     <dbl>
 1 Afghanistan  1980       NA        NA        NA        NA        NA        NA
 2 Afghanistan  1981       NA        NA        NA        NA        NA        NA
 3 Afghanistan  1982       NA        NA        NA        NA        NA        NA
 4 Afghanistan  1983       NA        NA        NA        NA        NA        NA
 5 Afghanistan  1984       NA        NA        NA        NA        NA        NA
 6 Afghanistan  1985       NA        NA        NA        NA        NA        NA
 7 Afghanistan  1986       NA        NA        NA        NA        NA        NA
 8 Afghanistan  1987       NA        NA        NA        NA        NA        NA
 9 Afghanistan  1988       NA        NA        NA        NA        NA        NA
10 Afghanistan  1989       NA        NA        NA        NA        NA        NA
# ℹ 7,230 more rows
# ℹ 50 more variables: sp_m_65 <dbl>, sp_f_014 <dbl>, sp_f_1524 <dbl>,
#   sp_f_2534 <dbl>, sp_f_3544 <dbl>, sp_f_4554 <dbl>, sp_f_5564 <dbl>,
#   sp_f_65 <dbl>, sn_m_014 <dbl>, sn_m_1524 <dbl>, sn_m_2534 <dbl>,
#   sn_m_3544 <dbl>, sn_m_4554 <dbl>, sn_m_5564 <dbl>, sn_m_65 <dbl>,
#   sn_f_014 <dbl>, sn_f_1524 <dbl>, sn_f_2534 <dbl>, sn_f_3544 <dbl>,
#   sn_f_4554 <dbl>, sn_f_5564 <dbl>, sn_f_65 <dbl>, ep_m_014 <dbl>, …

Échale un vistazo y piensa qué cosas te podría pedir para convertirlo a tidydata.

6.1 Pregunta 1

¿Qué significan los datos? ¿Cuántas variables y observaciones tenemos?

Code
library(tidyr)
? who2

6.2 Pregunta 2

¿Es tidy data? ¿Por qué

Code
# No es tidy data porque en realidad todas las variable
# a partir de year es lo mismo: casos de tuberculosis
# Son todo casos, solo que en distintas edades o sexo o tipos
# de diagnosis, pero la variable es "cases"

6.3 Pregunta 3

Primer paso para tidy data: pivota la tabla (consejo: usa papel y boli para bocetar como debería quedar la base de datos) tal que exista una columna llamada cases

Code
who_tidy <-
  who2 |> 
  # fíjate que en lugar de elegir las 56 columnas
  # le decimos las que NO queremos pivotar
  # los nombres de columnas los mandamos a "type" (tipo de caso)
  pivot_longer(cols = -(country:year),
               names_to = "type",
               values_to = "cases")
who_tidy 

6.4 Pregunta 4

Si te fijas hay muchísimas filas que no tiene sentido mantener ya que ¡no tenemos casos! Investiga las opciones de pivot_longer() para ver como podemos directamente eliminarlas en el pivotaje

Code
who_tidy <-
  who2 |> 
  # con values_drop_na = TRUE eliminamos los NA
  pivot_longer(cols = -(country:year),
               names_to = "type",
               values_to = "cases",
               values_drop_na = TRUE)
who_tidy 

6.5 Pregunta 5

Si te fijas ahora en type tenemos codificada la información como diagnosis_sexo_edad. ¿Cómo separarlo en 3 columnas? Investiga tanto separate() como las opciones de pivot_longer()

Code
# con separate
who_tidy <-
  who_tidy |> 
  separate(col = "type", into = c("diagnosis", "sex", "age"))

# con pivot_longer
who_tidy <-
  who2 |> 
  pivot_longer(cols = -(country:year),
               names_to = c("diagnosis", "sex", "age"),
               values_to = "cases",
               values_drop_na = TRUE,
               names_sep = "_")
who_tidy 

6.6 Pregunta 6

Por último, separa en dos (age_inf, age_sup) el tramo etario (que sean números). Piensa cómo hacerlo ya que no siempre son 4 números

Code
# Usamos separate y le indicamos las posiciones, pero desde atrás 
# ya que siempre el límite superior es un número de 2 cifras
# y usamos convert = TRUE para convertir a números
who_tidy <-
  who_tidy |> 
  separate(col = "age", into = c("age_inf", "age_sup"),
           sep = -2, convert = TRUE)