Sesión 7 Funciones e iteracion

“To understand computations in R, two slogans are helpful:
Everything that exists is an object.
Everything that happens is a function call.”
— John Chambers

Funciones

En R todas las operaciones son producto de la llamada a una función, esto incluye operaciones como +, operadores que controlan flujo como for, if y while, e incluso operadores para obtener subconjuntos como [ ] y $.

a <- 3
b <- 4
a + b
#> [1] 7
`+`(a, b)
#> [1] 7

for (i in 1:2) print(i)
#> [1] 1
#> [1] 2
`for`(i, 1:2, print(i))
#> [1] 1
#> [1] 2

Para escribir código eficiente y fácil de leer es importante saber escribir funciones, se dice que si hiciste copy-paste de una sección de tu código 3 o más veces es momento de escribir una función.

Escribimos una función para calcular un promedio ponderado:

wtd_mean <- function(x, wt = rep(1, length(x))) {
  sum(x * wt) / sum(wt)
}

Notemos que esta función recibe hasta dos argumentos:

x: el vector a partir del cual calcularemos el promedio y
wt: un vector de ponderadores para cada componente del vector x.

Notemos además que al segundo argumento le asignamos un valor predeterminado, esto implica que si no especificamos los ponderadores la función usará el valor predeterminado y promediara con mismo peso a todas las componentes.

wtd_mean(c(1:10))
#> [1] 5.5
wtd_mean(1:10, 10:1)
#> [1] 4

Veamos como escribir una función que reciba un vector y devuelva el mismo vector centrado en cero.

Comenzamos escribiendo el código para un caso particular, por ejemplo, reescalando el vector $(0, 5, 10)$.

vec <- c(0, 5, 10)
vec - mean(vec)
#> [1] -5  0  5

Una vez que lo probamos lo convertimos en función:

center_vector <- function(vec) {
  vec - mean(vec)
}
center_vector(c(0, 5, 10))
#> [1] -5  0  5

Ejercicio

Escribe una función que reciba un vector y devuelva el mismo vector reescalado al rango 0 a 1.

Comienza escribiendo el código para un caso particular, por ejemplo, empieza reescalando el vector . Tip: la función range() devuelve el rango de un vector.

Estructura de una función

Las funciones de R tienen tres partes:

El cuerpo: el código dentro de la función

body(wtd_mean)
#> {
#>     sum(x * wt)/sum(wt)
#> }

Los formales: la lista de argumentos que controlan como puedes llamar a la función,

formals(wtd_mean)
#> $x
#> 
#> 
#> $wt
#> rep(1, length(x))

El ambiente: el mapeo de la ubicación de las variables de la función.

library(ggplot2)
environment(wtd_mean)
#> <environment: R_GlobalEnv>
environment(ggplot)
#> <environment: namespace:ggplot2>

Veamos mas ejemplos, ¿qué regresan las siguientes funciones?

# 1
x <- 5
f <- function(){
  y <- 10
  c(x = x, y = y) 
}
rm(x, f)

# 2
x <- 5
g <- function(){
  x <- 20
  y <- 10
  c(x = x, y = y)
}
rm(x, g)

# 3
x <- 5
h <- function(){
  y <- 10
  i <- function(){
    z <- 20
    c(x = x, y = y, z = z)
  }
  i() 
}

# 4 ¿qué ocurre si la corremos por segunda vez?
j <- function(){
  if (!exists("a")){
    a <- 5
  } else{
    a <- a + 1 
}
  print(a) 
}
x <- 0
y <- 10

# 5 ¿qué regresa k()? ¿y k()()?
k <- function(){
  x <- 1
  function(){
    y <- 2
    x + y 
  }
}

Las reglas de búsqueda determinan como se busca el valor de una variable libre en una función. A nivel lenguaje R usa lexical scoping, esto implica que en R los valores de los símbolos se basan en como se anidan las funciones cuando fueron creadas y no en como son llamadas.

f <- function(x) {
  x + y
}

f(2)
#> Error in f(2): object 'y' not found

Si creamos el objeto y.

y <- 1
f(2)
#> [1] 3

Como consecuencia de las reglas de búsqueda de R, todos los objetos deben ser guardados en memoria y, si uno no es cuidadoso se pueden cometer errores fácilmente.

y <- 100
f(2)
#> [1] 102

Observaciones del uso de funciones

Cuando llamamos a una función podemos especificar los argumentos en base a posición, nombre completo o nombre parcial:

f <- function(abcdef, bcde1, bcde2) {
  c(a = abcdef, b1 = bcde1, b2 = bcde2)
}
# Posición
f(1, 2, 3)
#>  a b1 b2 
#>  1  2  3
f(2, 3, abcdef = 1)
#>  a b1 b2 
#>  1  2  3
# Podemos abreviar el nombre de los argumentos
f(2, 3, a = 1)
#>  a b1 b2 
#>  1  2  3
# Siempre y cuando la abreviación no sea ambigua
f(1, 3, b = 1)
#> Error in f(1, 3, b = 1): argument 3 matches multiple formal arguments

Los argumentos de las funciones en R se evalúan conforme se necesitan (lazy evaluation),

f <- function(a, b){
  a ^ 2
}
f(2)
#> [1] 4

La función anterior nunca utiliza el argumento b, de tal manera que f(2) no produce ningún error.

Funciones con el mismo nombre en distintos paquetes:

La función filter() (incluida en R base) aplica un filtro lineal a una serie de tiempo de una variable.

x <- 1:100
filter(x, rep(1, 3))
#> Error in UseMethod("filter"): no applicable method for 'filter' applied to an object of class "c('integer', 'numeric')"

Ahora cargamos dplyr.

library(dplyr)
filter(x, rep(1, 3))
#> Error in UseMethod("filter"): no applicable method for 'filter' applied to an object of class "c('integer', 'numeric')"

R tiene un conflicto en la función a llamar, nosotros requerimos usar filter de stats y no la función filter de dplyr. R utiliza por default la función que pertenece al último paquete que se cargó.

search()
#>  [1] ".GlobalEnv"        "package:lubridate" "package:forcats"  
#>  [4] "package:stringr"   "package:dplyr"     "package:purrr"    
#>  [7] "package:readr"     "package:tidyr"     "package:tibble"   
#> [10] "package:ggplot2"   "package:tidyverse" "package:stats"    
#> [13] "package:graphics"  "package:grDevices" "package:utils"    
#> [16] "package:datasets"  "package:methods"   "Autoloads"        
#> [19] "package:base"

Una opción es especificar el paquete en la llamada de la función:

stats::filter(x, rep(1, 3))
#> Time Series:
#> Start = 1 
#> End = 100 
#> Frequency = 1 
#>   [1]  NA   6   9  12  15  18  21  24  27  30  33  36  39  42  45  48  51  54
#>  [19]  57  60  63  66  69  72  75  78  81  84  87  90  93  96  99 102 105 108
#>  [37] 111 114 117 120 123 126 129 132 135 138 141 144 147 150 153 156 159 162
#>  [55] 165 168 171 174 177 180 183 186 189 192 195 198 201 204 207 210 213 216
#>  [73] 219 222 225 228 231 234 237 240 243 246 249 252 255 258 261 264 267 270
#>  [91] 273 276 279 282 285 288 291 294 297  NA

Como alternativa surge el paquete conflicted que alerta cuando hay conflictos y tiene funciones para especificar a que paquete se desea dar preferencia en una sesión de R.

Iteración

En análisis de datos es común implementar rutinas iteraivas, esto es, cuando debemos aplicar los mismos pasos a distintas entradas. Veremos que hay dos paradigmas de iteración:

Programación imperativa: ciclos for y ciclos while.
Programación funcional: los ciclos se implmentan mediante funciones,

La ventaja de la programación imperativa es que hacen la iteración de manera clara, sin embargo, veremos que una vez que nos familiarizamos con el paradigma de programación funcional, resulta en código más fácil de mantener y menos propenso a errores.

Ciclos for

Supongamos que tenemos una base de datos y queremos calcular la media de sus columnas numéricas.

df <- data.frame(id = 1:10, a = rnorm(10), b = rnorm(10, 2), c = rnorm(10, 3), 
    d = rnorm(10, 4))
df
#>    id          a          b        c        d
#> 1   1  1.0974355  1.7549684 3.768092 4.474497
#> 2   2  0.2383376  0.0877771 2.789389 2.620226
#> 3   3 -0.2600241  3.7346894 2.432871 5.710107
#> 4   4  0.8884876  2.5403336 3.889400 4.877224
#> 5   5  0.9268678  2.7673062 2.907595 3.984483
#> 6   6  0.9859076 -0.6774733 2.809649 2.957303
#> 7   7 -0.7168732  2.0166886 3.368445 3.758661
#> 8   8 -0.6908190  2.2058280 3.132857 2.808530
#> 9   9 -0.9696831  2.1189732 1.890708 2.721669
#> 10 10  3.0073544  2.5412250 1.062003 4.681626

Podemos crear el código para cada columna pero esto involucra copy-paste y no será muy práctico si aumenta el número de columnas:

mean(df$a)
#> [1] 0.4506991
mean(df$b)
#> [1] 1.909032
mean(df$c)
#> [1] 2.805101
mean(df$d)
#> [1] 3.859433

Con un ciclo for sería:

salida <- vector("double", 4)  
for (i in 1:4) {            
  salida[[i]] <- mean(df[[i+1]])      
}
salida
#> [1] 0.4506991 1.9090316 2.8051009 3.8594326

Los ciclos for tienen 3 componentes:

La salida: salida <- vector("double", 4). Es importante especificar el tamaño de la salida antes de iniciar el ciclo for, de lo contrario el código puede ser muy lento.
La secuencia: determina sobre que será la iteración, la función seq_along puede resultar útil.

salida <- vector("double", 5)  
for (i in seq_along(df)) {            
  salida[[i]] <- mean(df[[i]])      
}
seq_along(df)
#> [1] 1 2 3 4 5

El cuerpo: salida[[i]] <- mean(df[[i]]), el código que calcula lo que nos interesa sobre cada objeto en la iteración.

Ejercicio

Calcula el número de valores únicos en cada columna de los datos iris.

head(iris)
#>   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
#> 1          5.1         3.5          1.4         0.2  setosa
#> 2          4.9         3.0          1.4         0.2  setosa
#> 3          4.7         3.2          1.3         0.2  setosa
#> 4          4.6         3.1          1.5         0.2  setosa
#> 5          5.0         3.6          1.4         0.2  setosa
#> 6          5.4         3.9          1.7         0.4  setosa

Recordando la limpieza de datos de la sección anterior en uno de los últimos ejercicios leíamos archivos de manera iteativa, el vector paths contenía la ruta a distintos archivos csv. Crea la tabla de datos final usando un ciclo for.

paths <- list.files("data/specdata", pattern = ".csv", full.names = TRUE)

Programación funcional

Ahora veremos como abordar iteración usando programación funcional.

Regresando al ejemplo inicial de calcular la media de las columnas de una tabla de datos:

salida <- vector("double", 4)  
for (i in 1:4) {            
  salida[[i]] <- mean(df[[i+1]])      
}
salida
#> [1] 0.4506991 1.9090316 2.8051009 3.8594326

Podemos crear una función que calcula la media de las columnas de un data.frame:

col_media <- function(df) {
  salida <- vector("double", length(df))
  for (i in seq_along(df)) {
    salida[i] <- mean(df[[i]])
  }
  salida
}
col_media(df)
#> [1] 5.5000000 0.4506991 1.9090316 2.8051009 3.8594326
col_media(select(iris, -Species))
#> [1] 5.843333 3.057333 3.758000 1.199333

Y podemos extender a crear más funciones que describan los datos:

col_mediana <- function(df) {
  salida <- vector("double", length(df))
  for (i in seq_along(df)) {
    salida[i] <- median(df[[i]])
  }
  salida
}
col_sd <- function(df) {
  salida <- vector("double", length(df))
  for (i in seq_along(df)) {
    salida[i] <- sd(df[[i]])
  }
  salida
}

Podemos hacer nuestro código más general y compacto escribiendo una función que reciba los datos sobre los que queremos iterar y la función que queremos aplicar:

col_describe <- function(df, fun) {
  salida <- vector("double", length(df))
  for (i in seq_along(df)) {
    salida[i] <- fun(df[[i]])
  }
  salida
}
col_describe(df, median)
#> [1] 5.5000000 0.5634126 2.1624006 2.8586218 3.8715720
col_describe(df, mean)
#> [1] 5.5000000 0.4506991 1.9090316 2.8051009 3.8594326

Ahora utilizaremos esta idea de pasar funciones a funciones para eliminar los ciclos for.

La iteración a través de funciones es muy común en R, hay funciones para hacer esto en R base (lapply(), apply(), sapply()). Nosotros utilizaremos las funciones del paquete purrr,

La familia de funciones del paquete iteran siempre sobre un vector (vector atómico o lista), aplican una función a cada parte y regresan un nuevo vector de la misma longitud que el vector entrada. Cada función especifica en su nombre el tipo de salida:

map() devuelve una lista.
map_lgl() devuelve un vector lógico.
map_int() devuelve un vector entero.
map_dbl() devuelve un vector double.
map_chr() devuelve un vector caracter.

En el ejemplo de las medias, map puede recibir un data.frame (lista de vectores) y aplicará las funciones a las columnas del mismo.

library(purrr)
map_dbl(df, mean)
#>        id         a         b         c         d 
#> 5.5000000 0.4506991 1.9090316 2.8051009 3.8594326
map_dbl(select(iris, -Species), median)
#> Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
#>         5.80         3.00         4.35         1.30

Podemos crear un vector de tipo caracter que nos almacene el tipo de dato almacenado en cada columna de la tabla iris.

map_chr(iris, class)
#> Sepal.Length  Sepal.Width Petal.Length  Petal.Width      Species 
#>    "numeric"    "numeric"    "numeric"    "numeric"     "factor"

Si queremos calcular el número de faltantes en cada columna podemos hacer:

sum(is.na(airquality$Ozone))
#> [1] 37
map_int(airquality, function(x) sum(is.na(x)))
#>   Ozone Solar.R    Wind    Temp   Month     Day 
#>      37       7       0       0       0       0

O utilizar la notación .:

map_int(airquality, ~ sum(is.na(.)))
#>   Ozone Solar.R    Wind    Temp   Month     Day 
#>      37       7       0       0       0       0

Usaremos map para ajustar un modelo lineal a subconjuntos de los datos mtcars determinados por el cilindraje del motor.

models <- mtcars |> 
  group_split(cyl) |> 
  map(function(df) lm(mpg ~ wt, data = df))

Podemos usar la notación . para hacer código más corto:

models <- mtcars |> 
  group_split(cyl) |> 
  map(~lm(mpg ~ wt, data = .))

models
#> [[1]]
#> 
#> Call:
#> lm(formula = mpg ~ wt, data = .)
#> 
#> Coefficients:
#> (Intercept)           wt  
#>      39.571       -5.647  
#> 
#> 
#> [[2]]
#> 
#> Call:
#> lm(formula = mpg ~ wt, data = .)
#> 
#> Coefficients:
#> (Intercept)           wt  
#>       28.41        -2.78  
#> 
#> 
#> [[3]]
#> 
#> Call:
#> lm(formula = mpg ~ wt, data = .)
#> 
#> Coefficients:
#> (Intercept)           wt  
#>      23.868       -2.192

Usemos map_** para unir tablas de datos que están almacenadas en múltiples archivos csv.

names(paths) <- basename(paths)
specdata_us_vec <- map(paths, ~readr::read_csv(., col_types = "Tddi"), 
    .id = "filename")
specdata_us_vec[[2]]
#> # A tibble: 3,652 × 4
#>    Date                sulfate nitrate    ID
#>    <dttm>                <dbl>   <dbl> <int>
#>  1 2001-01-01 00:00:00      NA      NA     2
#>  2 2001-01-02 00:00:00      NA      NA     2
#>  3 2001-01-03 00:00:00      NA      NA     2
#>  4 2001-01-04 00:00:00      NA      NA     2
#>  5 2001-01-05 00:00:00      NA      NA     2
#>  6 2001-01-06 00:00:00      NA      NA     2
#>  7 2001-01-07 00:00:00      NA      NA     2
#>  8 2001-01-08 00:00:00      NA      NA     2
#>  9 2001-01-09 00:00:00      NA      NA     2
#> 10 2001-01-10 00:00:00      NA      NA     2
#> # ℹ 3,642 more rows
class(specdata_us_vec)
#> [1] "list"

En este caso es más apropiado usar map_df

specdata_us <- map_df(paths, ~readr::read_csv(., col_types = "Tddi"), 
    .id = "filename")

Ejercicio

Usa la función map_** para calcular el número de valores únicos en las columnas de iris.
Usa la función map_** para extraer el coeficiete de la variable wt para cada modelo:

models[[1]]$coefficients[2]
#>        wt 
#> -5.647025

Utiliza map_*** para crear un vector con el valor mínimo de cada columna de los datos airquality. ¿Cómo eliminas los faltantes?
Con el siguiente código agregamos una variable categórica de velocidad del aire a los datos airqualuty, para cada categoría ajusta un modelo lineal (con la función lm) para explorar la relación entre ozono y radiación solar.

airquality$Wind.cat <- cut(airquality$Wind, breaks = c(1.5, 9, 12, 21))

Puedes descargar un acordeón de purrr aquí.