Équipe ONCOSTAT
Modules de cours de R

Programmation fonctionnelle

Dan Chaltiel

Les boucles et les fonctions




Plan du module

  1. Les boucles
  2. Les fonctions

Note

Cette présentation est fortement inspirée de l’excellent site juba.github.io, en particulier les chapitres 14 (fonctions) et 18 (purrr). Allez voir ce site pour en apprendre encore plus!

Les boucles

C’est quoi, une boucle ?

  • Répéter le même code sur plusieurs éléments (DRY)
  • Remplacer du copier-coller par une structure générique
  • Améliore la robustesse (toujours)
  • Améliore la lisibilité (parfois)

Trois grandes façons en R

  • for
    • Style impératif, proche des autres langages
    • Très flexible mais peu lisible, parfois lent
  • sapply()
    • Style fonctionnel en base R
    • Simplification automatique des sorties, parfois piégeuse
  • purrr::map()
    • Style fonctionnel tidyverse
    • Type de sortie explicite
    • S’intègre naturellement dans les pipes et les workflows modernes

Exemple simple (vecteur)

On a un vecteur character, et on veut séparer chaque item en Nom/Prénom

library(tidyverse)
x = c("Jeanne Alyse", "Bertrand Domise", "XXXX", NA)

On va donc utiliser str_split() (package stringr) dans une boucle.

res_for = list()
for (i in seq_along(x)) {
  res_for[[i]] = stringr::str_split(x[i], " ")
}
res_for
#> [[1]]
#> [[1]][[1]]
#> [1] "Jeanne" "Alyse" 
#> 
#> 
#> [[2]]
#> [[2]][[1]]
#> [1] "Bertrand" "Domise"  
#> 
#> 
#> [[3]]
#> [[3]][[1]]
#> [1] "XXXX"
#> 
#> 
#> [[4]]
#> [[4]][[1]]
#> [1] NA
res_sapply = sapply(x, function(s) stringr::str_split(s, " "))
res_sapply
#> $`Jeanne Alyse`
#> [1] "Jeanne" "Alyse" 
#> 
#> $`Bertrand Domise`
#> [1] "Bertrand" "Domise"  
#> 
#> $XXXX
#> [1] "XXXX"
#> 
#> $<NA>
#> [1] NA
res_purrr = x %>% purrr::map(~ stringr::str_split(.x, " "))
res_purrr
#> [[1]]
#> [[1]][[1]]
#> [1] "Jeanne" "Alyse" 
#> 
#> 
#> [[2]]
#> [[2]][[1]]
#> [1] "Bertrand" "Domise"  
#> 
#> 
#> [[3]]
#> [[3]][[1]]
#> [1] "XXXX"
#> 
#> 
#> [[4]]
#> [[4]][[1]]
#> [1] NA

Pour ou contre: les boucles for

👎 Contre :

  • Beaucoup de code pour pas grand chose (5 lignes vs 2)
  • Plus ou moins lisible selon la complexité
  • Ne s’applique pas en pipeline
  • Déclaration manuelle de l’objet final (souvent mal optimisé)
  • Programmation non-fonctionnelle → Effets de bords !

👍 Pour :

  • Très flexible
  • Effets de bords ?

Exemple de boucle for maléfique

res_for = list()
for (i in seq_along(x)) {
  x = stringr::str_split(x[i], " ") #variable temporaire
  res_for[i] = x
}
res_for
#> [[1]]
#> [1] "Jeanne" "Alyse" 
#> 
#> [[2]]
#> [1] "NULL"
#> 
#> [[3]]
#> [1] "NULL"
#> 
#> [[4]]
#> [1] "NULL"

La ligne x = str_split(...) a eu un effet de bord, c’est à dire un effet hors de sa zone, et ça a écrasé notre vecteur d’origine 😱😭

Dans une très grande boucle for, ce genre de conflit de noms peut arriver !

Pour ou contre: sapply()

Maintenant, on veut appliquer table() à chaque colonne d’une dataframe :

library(crosstable)
x = mtcars2 %>% select(am, vs)
sapply(x, table)
#>        am vs
#> auto   19 14
#> manual 13 18
map(x, table)
#> $am
#> 
#>   auto manual 
#>     19     13 
#> 
#> $vs
#> 
#> straight  vshaped 
#>       14       18

Par défaut, sapply() simplifie implicitement l’output en vecteur ou en matrice. Comme ici, ça peut donner lieu à des étrangetés.

purrr à la rescousse!

À l’opposé, les fonctions de purrr assurent à 100% le type de l’output:

  • map() retourne une liste
  • map_chr() retourne un vecteur character
  • map_dbl() retourne un vecteur numeric (double)
  • map_lgl() retourne un vecteur logical

Ça donne aussi accès à la syntaxe “lambda-function” : argument par défaut = .x

mon_vecteur = 1:5
map_dbl(mon_vecteur, ~.x+1)
#> [1] 2 3 4 5 6
map_dbl(mon_vecteur, function(x) x+1)
#> [1] 2 3 4 5 6
map_dbl(mon_vecteur, \(x) x+1)  #depuis R 4.1
#> [1] 2 3 4 5 6
mon_vecteur+1
#> [1] 2 3 4 5 6

Quand utiliser des boucles ?

Sur des fonctions vectorisées, on n’utilise pas de boucles :

textes = c("fantastique", "effectivement", "igloo")
str_count(textes, "f")
#> [1] 1 2 0
map_int(textes, ~ str_count(.x, "f")) #moins lisible + moins rapide
#> [1] 1 2 0

On l’utilise sur les fonctions non vectorisées :

fichiers = c("fichier1.csv", "fichier2.csv")
l = fichiers %>% map(read.csv2)

Dans str_count(), l’argument pattern n’est pas vectorisé :

c(a = "a", e = "e") %>% 
  map(~ str_count(textes, pattern = .x))
#> $a
#> [1] 2 0 0
#> 
#> $e
#> [1] 1 4 0

Oui mais, et les noms ?

Un des “avantages” de sapply(), c’est que les noms sont automatiques :

x = c("Jeanne Alyse", "Bertrand Domise", "XXXX", NA)
get_prenom = function(s) stringr::str_split(s, " ")[[1]][1]
sapply(x, get_prenom) 
#>    Jeanne Alyse Bertrand Domise            XXXX            <NA> 
#>        "Jeanne"      "Bertrand"          "XXXX"              NA

Ca peut poser un problème en cas de noms très longs, donc map() ne prend que les noms explicites :

x2 = c("JA"="Jeanne Alyse", BD="Bertrand Domise", XX="Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua", KO=NA)
map_chr(x2, get_prenom) 
#>         JA         BD         XX         KO 
#>   "Jeanne" "Bertrand"    "Lorem"         NA

Pro tip:

La fonction set_names() assigne à un objet des noms basés sur ses propres valeurs.

x %>% set_names()
#>      Jeanne Alyse   Bertrand Domise              XXXX              <NA> 
#>    "Jeanne Alyse" "Bertrand Domise"            "XXXX"                NA

D’ailleurs, c’est cool les noms, non ?

La fonction imap() permet de boucler sur un vecteur (.x) et ses noms (.y) en même temps.

mtcars %>% 
  select(1:3) %>% 
  imap_chr(~{
    paste(.y, ":", n_distinct(.x), "valeurs distinctes")
  })
#>                            mpg                            cyl                           disp 
#>  "mpg : 25 valeurs distinctes"   "cyl : 3 valeurs distinctes" "disp : 27 valeurs distinctes"

Note

Notez que boucler sur une dataframe permet d’itérer sur les colonnes en conservant les noms.

Parallélisation

A partir de purrr v1.1.0, on peut utiliser in_parallel() pour paralléliser nos fonctions.

library(tictoc)
very_long_addition = function(x, y){
  Sys.sleep(0.1) #on attend 0.1 sec
  x+y
}

tic("Sans parallelisation")
a = map_dbl(1:20, ~ very_long_addition(.x, 1))
toc()
#> Sans parallelisation: 2.19 sec elapsed
tic("Avec parallelisation")
mirai::daemons(4) #on utilise 4 coeurs

crate = in_parallel(
  ~ very_long_addition(.x, 1), 
  very_long_addition=very_long_addition
)

b = map_dbl(1:20, crate)

mirai::daemons(0)
toc()
#> Avec parallelisation: 1.25 sec elapsed

Ici, on double presque notre vitesse, et ce serait encore plus impressionnant sur plus d’itérations ou une fonction plus longue !

Les fonctions

Principe

“En R, tout ce qui existe est un objet, et tout ce qui se passe est une fonction.”
John Chambers

Une fonction prend des arguments, les transforme, et retourne un résultat :

n_events_schoenfeld = function(za, zb, hr){
  x1 = 4*(za+zb)^2
  x2 = log(hr)^2
  rslt = x1/x2
  return(rslt)
}
n_events_schoenfeld(za=1.96, zb=1.28, hr=0.7)
#> [1] 330.0691
n_events_schoenfeld(za=1.96, zb=1.28, hr=0.6)
#> [1] 160.918
n_events_schoenfeld = function(za, zb, hr){
  x1 = 4*(za+zb)^2
  x2 = log(hr)^2
  x1/x2
}
n_events_schoenfeld(za=1.96, zb=1.28, hr=0.7)
#> [1] 330.0691
n_events_schoenfeld(za=1.96, zb=1.28, hr=0.6)
#> [1] 160.918

Par défaut, la dernière ligne vaut return() donc on peut simplifier.

n_events_schoenfeld = function(za, zb, hr){
  4*(za+zb)^2 / (log(hr)^2)
}
n_events_schoenfeld(za=1.96, zb=1.28, hr=0.7)
#> [1] 330.0691
n_events_schoenfeld(za=1.96, zb=1.28, hr=0.6)
#> [1] 160.918

On peut toujours faire encore plus court, mais essayez de privilégier la lisibilité !

Note

Une fonction dont le résultat ne dépend que de ses arguments est appelée une “fonction pure”.

DRY: Don’t repeat yourself

  • Règle générale de programmation

  • Au 3ème copié-collé d’un même bout de code, on se force à faire une fonction !

  • On lui donne un nom explicite.

  • Ça améliore la lisibilité, la maintenabilité, et la robustesse du code.

Bonnes pratiques

Comparons deux fonctions :

a = 1
very_bad_function = function(x){
  x + a + rnorm(1)
}
very_bad_function(x=1)
#> [1] 2.235631
a=100
very_bad_function(x=1)
#> [1] 100.8557
very_bad_function(x=1)
#> [1] 98.88467
add_with_noise = function(x, y, seed=1234){
  set.seed(seed)
  x + y + rnorm(1)
}
add_with_noise(x=1, y=1)
#> [1] 0.7929343
add_with_noise(x=1, y=100)
#> [1] 99.79293
add_with_noise(x=1, y=100)
#> [1] 99.79293

Pourquoi la fonction de gauche est fort vilaine ?

  1. Elle utilise une variable globale, a, qui peut changer / ne pas exister
  2. Elle n’est pas déterministe / reproductible
  3. Elle a un nom qui ne donne aucune information

Wrappers

Un wrapper est une fonction courte qui vise à “réécrire” une fonction connue pour l’appliquer à un cas particulier.

Exemple de EDCimport::fct_yesno() :

fct_yesno = function(x){
  x = as.character(x)
  fct_recode(x, "Yes" = "1", "No" = "0")
}
input = c(0,1,0,0)
fct_yesno(input)
#> [1] No  Yes No  No 
#> Levels: No Yes

Tip

La vraie fonction dans EDCimport est bien plus utile, elle gère les 0/1, les TRUE/FALSE, mais aussi les Yes/No encodés, par exemple 1-Yes/2-No ou encore les Oui/Non.

L’ellipsis ...

L’ellipsis (...) s’utilise quand on ne sait pas de combien d’arguments on a besoin.
C’est le cas par exemple avec dplyr::select().

C’est surtout utile avec les wrappers pour passer des arguments :

moyenne = function(v, ...){
  mean(v, ...) #on passe ... à mean()
}
a = c(1:10, NA)
moyenne(a, na.rm=TRUE)
#> [1] 5.5
moyenne(a, na.rm=FALSE)
#> [1] NA

On en a rarement besoin directement, mais on l’utilise avec c(...) ou list(...) :

f = function(...){
  x = c(...)
  paste(names(x), x, sep=":", collapse=" -- ")
}
f(a="A", b="B", c="C")
#> [1] "a:A -- b:B -- c:C"
f(alpha=0.95, power=0.90, na.rm=TRUE, type="Schoenfeld")
#> [1] "alpha:0.95 -- power:0.9 -- na.rm:TRUE -- type:Schoenfeld"

Appel par ordre vs Appel par nom

moyenne = function(x, remove_missing=TRUE){
  mean(x, na.rm=remove_missing)
}
a = c(1:10, NA)
moyenne(a, TRUE)
#> [1] 5.5
moyenne(x=a, remove_missing=TRUE)
#> [1] 5.5
moyenne(x=a, remove_missing=FALSE)
#> [1] NA

En général, on n’appelle pas le 1er argument par nom car peu informatif.

Certaines fonctions obligent à appeler par nom certains arguments.

Notez la valeur par défaut de remove_missing qui permet d’écrire simplement:

moyenne(a)
#> [1] 5.5
mean(a)
#> [1] NA

Erreurs et warnings

Contrairement à SAS, les erreurs R sont bloquantes et arrêtent le script : c’est le signe d’un problème majeur qu’il faut corriger.
Les warnings ne sont pas bloquants.

On les émet avec stop() et warning() (base R) ou cli_abort()et cli_warn() (package cli).

n_events_schoenfeld  = function(za, zb, hr){
  if(za<0) stop("`za` doit être positif!")
  if(zb<0) stop("`zb` doit être positif!")
  if(hr<0 || hr>1) stop("`hr` doit être compris entre 0 et 1!")
  if(za>10) warning("`za` est vraiment très grand, t'es sûr de toi ?")
  4*(za+zb)^2 / (log(hr)^2)
}
n_events_schoenfeld(za=-1.96, zb=1.28, hr=0.7)
#> Error in n_events_schoenfeld(za = -1.96, zb = 1.28, hr = 0.7): `za` doit être positif!
n_events_schoenfeld(za=19.6, zb=1.28, hr=0.7)
#> Warning in n_events_schoenfeld(za = 19.6, zb = 1.28, hr = 0.7): `za` est vraiment très grand, t'es
#> sûr de toi ?
#> [1] 13708.05
library(cli)
n_events_schoenfeld  = function(za, zb, hr){
  if(za<0) cli_abort("{.arg za} doit être positif, il vaut {.val {za}}!")
  if(za>10) cli_warn("{.arg za} est vraiment très grand ({.val {za}}), t'es sûr de toi ?")
  4*(za+zb)^2 / (log(hr)^2)
}
n_events_schoenfeld(za=-1.96, zb=1.28, hr=0.7)
#> Error in `n_events_schoenfeld()`:
#> ! `za` doit être positif, il vaut -1.96!
n_events_schoenfeld(za=19.6, zb=1.28, hr=0.7)
#> Warning: `za` est vraiment très grand (19.6), t'es sûr de toi ?
#> [1] 13708.05

Side-effects

Certaines fonctions ont une action concrete, par exemple print() ou write.csv().

On appelle ça des “effets de bord”.

On a rarement besoin d’en écrire nous-même.

Exemple final: simulation

On va prendre un exemple un peu compliqué pour revoir tout ce qu’on a appris !

simulate_lm = function(n, beta=0, sd=1, seed=NULL){
  set.seed(seed)
  data = tibble(
    x = rnorm(n),
    y = beta * x + rnorm(n, sd=sd)
  )
  
  lm(y ~ x, data = data) %>% 
    broom::tidy() %>% 
    filter(term=="x") %>% #on enlève l'intercept
    mutate(args_beta=beta, args_sd=sd)
}
simulate_lm(n=10000, beta=1) %>% as.data.frame()
#>   term estimate   std.error statistic p.value args_beta args_sd
#> 1    x 1.002794 0.009959544  100.6867       0         1       1

Warning

Attention set.seed(NULL) dans une fonction est parfois déconseillé, je ne suis pas expert !

run_simulations = function(N, n, beta = 0, sd = 1){
  tibble(seed = seq_len(N)) %>% 
    mutate(res = map(seed, ~simulate_lm(n=n, beta=beta, sd=sd, seed=.x)))
}


run_simulations(N=10, n=10, beta=0) %>% head(2)
#> # A tibble: 2 × 2
#>    seed res             
#>   <int> <list>          
#> 1     1 <tibble [1 × 7]>
#> 2     2 <tibble [1 × 7]>

run_simulations(N=10, n=10, beta=0) %>% unnest(res) %>% head(2)
#> # A tibble: 2 × 8
#>    seed term  estimate std.error statistic p.value args_beta args_sd
#>   <int> <chr>    <dbl>     <dbl>     <dbl>   <dbl>     <dbl>   <dbl>
#> 1     1 x       -0.516     0.449    -1.15    0.283         0       1
#> 2     2 x        0.240     0.415     0.577   0.580         0       1
show_results = function(sim){
  sim %>% 
    unnest(res) %>% 
    summarise(avg_beta = mean(estimate)) %>% 
    pull(avg_beta)
}


run_simulations(N=10, n=10, beta=0) %>% show_results()   #peu précis avec n=10
#> [1] -0.1713649

run_simulations(N=100, n=100, beta=0) %>% show_results() #mieux
#> [1] 0.007437185

run_simulations(N=100, n=100, beta=5) %>% show_results() #mieux aussi
#> [1] 5.007437

Débuguer avec browser()

Imaginons une fonction qui fait un calcul métier un peu obscur:

library(purrr)
f = function(bias, penalty, values){
  ref = mean(values) + bias - penalty
  values %>% 
    map(~{
      data.frame(
        input = .x,
        score = (.x + bias) / (.x - ref)
      )
    }) %>% 
    list_rbind()
}
f(bias=2, penalty=1, values=c(3, 4, 5, 6))
#>   input score
#> 1     3    -2
#> 2     4    -4
#> 3     5   -14
#> 4     6    16
f(bias=1, penalty=1, values=c(2, 4, 6))
#>   input score
#> 1     2  -1.5
#> 2     4   Inf
#> 3     6   3.5

Essayons d’ajouter if(.x - ref == 0) browser() dans la boucle.

Conclusion

Les boucles c’est génial

  • Ca réduit la duplication du code (DRY)
  • Ca améliore la maintenabilité du code
  • Ca améliore (parfois) la lisibilité du code

Les fonctions c’est génial

  • Ca réduit la duplication du code (DRY)
  • Ca améliore la maintenabilité du code
  • Ca améliore (parfois) la lisibilité du code